Kawazoe(@riverplus)さんのCodeIqの問題です。
問題はこちらの御世話になりました。解答例もあります。
->https://sw-logic.com/Portfolio/Programing/CodeIQ/2903.php
こちらにも情報があります。
->https://togetter.com/li/1049668
テストケースはあちこちから集めた気がします。
例えばz=x^2*y^3(x,y素因数)とすると、zの約数は、
x^0*y^0,x^0*y^1,x^0*y^2,x^0*y^3
x^1*y^0,x^1*y^1,x^1*y^2,x^1*y^3
x^2*y^0,x^2*y^1,x^2*y^2,x^2*y^3
となり合計は
x^0*Σ(0<=k<=3)(y^k)+x^1*Σ(0<=k<=3)(y^k)+x^2*Σ(0<=k<=3)(y^k)
=Σ(0<=h<=2)(x^h)*Σ(0<=k<=3)(y^k)
素因数やその指数が増えても同様です。
更にΣ(0<=h<=m)(x^h)=(x^(m+1)-1)/(x-1)
そしてn!^nのある素因数の指数はn!のその素因数の指数のn倍です。
Juliaではfactorが使えますがパッケージをインストールしないといけませんので、
自前で(n!)^nの素因数を調べるfactor1を実装しました。
n!の素因数分解はルジャンドル関数を用いて求めることができます。
つまりn!を素因数分解すると素数pの数はΣ(p<=p^i<=n)div(n,p^i)
factor1では同時に素数も求めています。n以下のすべての素数が素因数です。
その後約数の和の計算をします。結果は大きな数になりますのでBigIntを使います。
(BigIntを使わないとinvmodという関数が必要のよう)
2024.9.7追加->カッコ内はinvmodを使った場合のみ解けたという意味です。
#import Pkg;Pkg.add("Primes")
#import Primes: isprime, primes, primesmask, factor
function solve46(n,a)
function factor1(n)
#ルジャンドル関数を利用
for i in 2 : n
x = 0
if ar2[i] != 0 # iが素数なら
k = i
for j in i : i : n
ar2[j] = 0 # iの倍数は素数ではない
if j == k # jはiの累乗
x += div(n,j)
k *= i # kはiの累乗
end
end
ar1[i] = x*n
end
end
end
ar1 = zeros(Int,n)
ar2 = Array(1:1:n)
m = 1000003
ret = big(1)
factor1(n)
for i in 2 : n
if ar1[i] != 0
x = big(i)
ret = mod(ret * div(x^(ar1[x]+1)-1,x-1),m)
end
end
s = a == ret ? "ok " : "no "
println(s,ret)
end
solve46(2, 7)
solve46(3, 600)
solve46(4, 991111)
solve46(12, 845575)
solve46(630, 951018)
solve46(997, 772380)
solve46(1000, 985788)
Prologではコードが複雑にならないように先に素数を求めています。
それでfactor/4はルジャンドル関数だけを実装していますので
NをK^iで割るのではなく、
NをKで割った商をさらにKで割りつつ商を合計しています。
その後calc44で約数の和を計算します。
sumfac/2では素数のリストと、素因数の数を入れるリストを作ってfac/4を呼び、
fac4ではすべての素数をfactor/4に送り、結果をN倍してリストに入れています。
SWI-Prologは多倍長整数ですので、大きな数でも問題ありません。
（投稿時私の環境ではideoneでJavaしか表示されず実行できません）
%SWI-Prolog ver 7.42,7.64,9.25
%:-initialization(start). %(Online Prolog Compiler.)
%start.
factor(N,X,R,R):-X>N,!.
factor(N,X,R1,R):-N1 is N div X,R2 is R1+N1,factor(N1,X,R2,R).
fac([],_,R,R):-!.
fac([H|P],N,L,R):-
factor(N,H,0,V),V1 is V*N,
nth1(H,L,_,L1),nth1(H,R1,V1,L1),fac(P,N,R1,R),!.
sumfac(N,R):-length(L,N),maplist(=(0),L),prime(N,P),fac(P,N,L,R).
prime(2,[2]).
prime(N,R):-
N1 is N-1,prime(N1,R1),
(forall(member(X,R1),(N mod X=\=0))->R=[N|R1];R=R1).
calc(1,_,1).
calc(N,L,R):-
N1 is N-1,calc(N1,L,R1),nth1(N,L,X),
(X=:=0->Y=1;Y is (N^(X+1)-1) div (N-1)),R is (R1*Y) mod 1000003.
solve(N,R):-sumfac(N,L),calc(N,L,R1),R is R1,!.% mod 1000003,!.
start:-
f(N,A),(solve(N,R)->(R==A->Str=" ok ";Str=" no ");
write(fail)),write(" "),write(Str),writeln(R),fail.
start.
f(2, 7).
f(3, 600).
f(4, 991111).
f(12, 845575).
f(630, 951018).
f(997, 772380).
f(1000, 985788).