Skip to content

Instantly share code, notes, and snippets.

@shlomif

shlomif/pari.gp

Created July 11, 2019 23:06
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save shlomif/e7a5214f608d0e69a14a244187ec0cee to your computer and use it in GitHub Desktop.
Save shlomif/e7a5214f608d0e69a14a244187ec0cee to your computer and use it in GitHub Desktop.
Download ZIP
Raw
pari.gp
\\ Proposed PARI/GP slution for https://projecteuler.net/problem=483
\\ Based on https://oeis.org/A060014
\\ Contact: https://www.shlomifish.org/
T(n, k)={n!*polcoeff(sumdiv(k, i, moebius(k/i)*exp(sumdiv(i, j, x^j/j) + O(x*x^n))), n)}
Tlog(n, k)={log(0.000000000001+n!*polcoeff(sumdiv(k, i, moebius(k/i)*exp(sumdiv(i, j, x^j/j) + O(x*x^n))), n))}
a(n)=sum(k=1,n*n,k*T(n,k))
b(n)=sum(k=1,n*n,k*k*T(n,k))
for(n=1, 22,print(b(n)))
d(n,k,s)={t=k*k*exp(Tlog(n,k));s+=t;printp(" ",n," ",n*n," ",k," ",t," ",s);s}
c(n)={s=0;for(k=1,n*n,s=d(n,k,s));s}
\\ (* for(n=1, 22,if(b(n)!=c(n),pari_err_BUG(5))) *)
e()=print(c(350))
e()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment