CodeIQ:「セカンド・イクエイション」の解答

Kawazoe(@riverplus)さんのCodeIQの問題です。
こちらに問題と解答例があります。テストケースも頂いた気がするのですが。
->https://sangaku0418.hatenablog.com/entry/2015/11/17/131915
こちらにも情報があります。
->https://posfie.com/@riverplus/p/gQvM74G

問題の要約
　自然数a,b,cからなるa*x*x+b*x+c=0という方程式が実数根を持つとき
　b<=mの範囲でとれるa,b､cの組み合わせの数F(m)を求めよ。

F(n)=F(n-1)+(b==nの時の(a,c)の組み合わせの数)　です。
b==nの時
実根を持つ条件は、判別式よりn*n-4*a*c>=0、よってa*c<=n*n/4
なので、この条件を満たす自然数 a,c の組が何個あるかという問題
(a,c)は対称なのでまずa<cの場合を数えて2倍します。
この時a^2<a*c<=n^2/4、したがって1<=a<=n div 2
この範囲で c<=(n^2 div 4*a)、a<cなので　cの数は(n^2 div 4*a)-a
これを二倍して足していき最後に c=a の時の1を足します。

Prologで解いていました。
solve1/3で n を、fn/4で a を動かしています。
これでもなんとかテストケースはパスしますが、１桁多くなると解けません。

%SWI-Prolog ver 7.42,7.64,9.29
%:-initialization(start).   %ideone,Online Prolog Compiler
%start.

solve1(1,R,R).
solve1(M,R1,R):-
   X is M^2 div 4,Y is M div 2,fn(Y,X,0,R2),
   R3 is R1+R2,M1 is M-1,solve1(M1,R3,R),!.

fn(0,_,R,R):-!.    %X div Y は　ｎ^2 div (a*4) と同じ
fn(Y,X,R1,R):-
   V is (X div Y-Y)*2+1,R2 is R1 + V,Y1 is Y-1,fn(Y1,X,R2,R).

solve(M:C):-
   (solve1(M,0,R)->(R=:=C->Str=" ok  ";Str=" no  ");
   write(fail)),write(" "),write(Str),writeln(R),!.

start:-f(A:B),solve(A:B),fail.
start.

f(3:4).
f(4:12).
f(10:287). 
f(20:3032).
f(67:167672).
f(100:620174). 
f(125:1280414).
f(1879:7291768620).
f(2615:20635666540).
f(3000:31773763314).

他に方法がないか問題を眺めていたら、
b=n の時の (a,c) の数は 1 から x=div(n^2,4) までの各整数について、
その整数の約数の個数を全部足したものと同じ事に気が付きました。
これはΣ(1<=j<=x) div(x,j) と同じです。ただこのままでは速くありません。
そこでx以下の数の約数の総和をO(√x)で求めるアルゴリズムを利用します。
->https://telvia.jp/mathB/suretsu/suretsu-z05/
->https://en.wikipedia.org/wiki/Dirichlet_hyperbola_method
約数を足し合わせることは、グラフで不等式 j*k<=x　(j,k>0) を満たす
格子点 (j,k) の個数を数えることと同じです。
(j,k) と (k,j) は対で現れますので、y=foor(√x) としますと、
1<=j<=y の k の数をdiv(x,j) で求めて2倍し足していきます。
(1,1) から (y,y) の部分は重複しますので重複分を後から引きます。　
Juliaで書きました。

#julia ver 1.0,1.41,1.53,1.85,1.12.1

function solve(n0)
    sum=0
    for i in 1 : n0
        x = div(i^2,4)
        m = Int(floor(sqrt(x)))
        for j in 1 : m
           sum += div(x,j)*2
        end
        sum -= m*m
    end
    return sum
end
function f(n,a)
　　ret = solve(n)
   s = ret == a ? "ok " : "no "
   println(s,n,"->",ret)
end

f(3,4)
f(4,12)
f(10,287)
f(20,3032)
f(67,167672)
f(100,620174)
f(125,1280414)
f(1879,7291768620)
f(2615,20635666540)
f(3000,31773763314)