CodeIQ:「LCM・パレード」問題の解答

Kawazoe(@riverplus)さんのCodeIqの問題です。
問題とテストケースはこちらにお世話になりました。解答もあります。
->https://blog.goo.ne.jp/r-de-r/e/cc4c0751d6d5658a45dbe573488e2e8f
こちらにも情報があります。
->https://togetter.com/li/1198588

aの10^8を何とかしないと計算量が多すぎますので、aをbで表すことを考えます。
そこでまず
lcm(x*j+i,x)=(x*j+i)*x/gcd(x*j+i,x)=(x*j+i)*x/gcd(i,x)
そして
y=d*x+rとして1<=k<=d*xとするとv=Σlcm(k,x)/xは
for i in 1:x , j in 0:d-1
　 v+=lcm(x*j+i,x)/x

ここでlcm(x*j+i,x)/x=( (x*j+i)*x/gcd(x,i) )/x=(x*j+i)/gcd(x,i)

またΣx*j(j=0..d-1)=x*d*(d-1)/2、Σi(j=0..d-1)=d*i となりますので上の式は
for i in 1:x
　v+=(d*(d-1)/2*x+d*i)/gcd(x,i)
剰余分は別に
for i in 1:r
　v+=(d*x+i)/gcd(x,i)
となります。
ところで投稿のため見直していたら、出題者の「 O(k)とO(√k)の解法を紹介したよ」
というtweetを見つけました。 O(√k)というのは見当がつきません。
まずJuliaから。

#julia ver 1.41,1.53,1.73

function solve(y,x,a)
    ans = 0
    d,r = divrem(y,x)
    m = div(d*(d-1),2)
    for i in 1:x
        ans += div((m*x+d*i),gcd(x,i))
        if i <= r
           ans += div(d*x+i,gcd(x,i))
        end
    end
    s = a == ans ? "ok " : "no "
   println(s,ans)
end

solve(10,3,43)
solve(20,9,162)
solve(9, 12, 22)
solve(40, 50, 502)
solve(312, 41, 47708)
solve(7122, 360, 10778909)
solve(1234567, 17200, 414701577895)
solve(98777, 98999, 4878497253)
solve(100000000, 98765, 4199787396686968)

Prologも同じ方針で解いています。

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

loop(X,X,_,_,_,R,R):-!.
loop(I0,X,D,M,N,R1,R):-
    I is I0+1,Y is gcd(X,I),R2 is R1+(N*X+D*I) div Y,
    (I =< M->R3 is R2+(D*X+I) div Y;R3=R2),loop(I,X,D,M,N,R3,R).

solve(Y,X,A):-
    divmod(Y,X,D,M),N is D*(D-1) div 2,
    (loop(0,X,D,M,N,0,R)->(R==A->Str=" ok  ";Str=" no  ");
    write(fail)),write(" "),write(Str),writeln(R),!.

start:-f(Y,X,A),solve(Y,X,A),fail.
start.

f(10,3,43).
f(20,9,162).
f(9, 12, 22).
f(40, 50, 502).
f(312, 41, 47708).
f(7122, 360, 10778909).
f(1234567, 17200, 414701577895).
f(98777, 98999, 4878497253).
f(100000000, 98765, 4199787396686968).