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