|
|
> No.3881[元記事へ]
GAIさんへのお返事です。
k=p[1]^a[1] * p[2]^a[2] * … * p[n]^a[n] と素因数分解されるとすると、
k' = k*(a[1]/p[1] + a[2]/p[2] + … + a[n]/p[n]) と表される。
なので、
素因数分解のプログラムで処理することができます。
配列（10^7程度）が確保できれば、篩い法で処理するのが速いと思います。
どちらにしても、5*10^15は天文学的な数です。
LET t0=TIME
LET K=10^6 !2からkまで
DIM D(K) !k'
MAT D=CON
LET S=0 !Σgcd(k,k')
FOR N=2 TO K !エラトステネスの篩いによる
!!!PRINT N;D(N) !k,k' debug
IF D(N)=1 THEN !素数なら
LET S=S+1 !p'=1
FOR J=N*N TO K STEP N !倍数を篩う
LET dK=J/N+N*D(J/N) !k'=a'b+ab'
LET D(J)=dK
NEXT J
ELSE !合成数なら
LET S=S+GCD(N,D(N))
END IF
NEXT N
PRINT S !結果を表示する
PRINT TIME-t0
END
EXTERNAL FUNCTION gcd(a,b) !最大公約数 a,b≧0
DO UNTIL b=0
LET t=b
LET b=MOD(a,b)
LET a=t
LOOP
LET gcd=a
END FUNCTION
|
|