| 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.
.-. .-.
/ \ .-. .-. / \
/ \ / \ .-. _ .-. / \ / \
/ \ / \ / \ / \ / \ / \ / \
/ \ / \ / `-' `-' \ / \ / \
\ / `-' `-' \ /
`-' `-'