Dr. Robert Harley on Tue, 6 May 2003 13:05:12 +0200 (CEST) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: Rademacher formula for p(n) with PARI |
Damien Wyart wrote: >I [...] plan to implement Rademacher's formula for p(n) using PARI. [...] >Any hints or code snippets would then be very welcome. Here is a simple implementation for GP: ------------------------------------------------------------------------------ partfac1(q, n) = real(sum(p = 0, q-1, if(gcd(p, q)>1, 0, exp(Pi*I*(sum(m = 1, q-1, m*(m*p/q-floor(m*p/q)-1/2))-2*n*p)/q)))) partfac2(q, n) = local(K, l); K = Pi*sqrt(2/3); l = sqrt(n-1/24); (cosh(K*l/q)*K/q-sinh(K*l/q)/l)*sqrt(q/2)/2/Pi/(n-1/24) partbla(d) = local(K); K = Pi*sqrt(2/3); sqrt(27)*K^6/d^2*(sinh(K*d)/K^3/d^3+1/6)^3 partition(n) = { local(prec, oprec, delta, mq, d, t); if ( n == 0 , 1 , oprec = precision(1.0) ; prec = ceil(log(exp(Pi*sqrt(2*n/3))/4/sqrt(3)/n)/log(10))+10 ; if ( n < 200 , delta = 2 , default(realprecision, 20) ; d = 2 ; while(partbla(d) < n, d++) ; delta = solve ( delta = d-1, d, partbla(delta)-n ) ) ; mq = ceil(sqrt(n)/delta) ; default(realprecision, prec) ; t = round(sum(q = 1, mq, partfac1(q, n)*partfac2(q, n))) ; default(realprecision, oprec) ; t ) } ------------------------------------------------------------------------------ For instance: ------------------------------------------------------------------------------ > # timer = 1 (on) > partition(12345) time = 152 ms. %1 = 69420357953926116819562977205209384460667673094671463620270321700806074195845953959951425306140971942519870679768681736 ------------------------------------------------------------------------------ Regards, Rob. .-. .-. / \ .-. .-. / \ / \ / \ .-. _ .-. / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / `-' `-' \ / \ / \ \ / `-' `-' \ / `-' `-'