Conference rusure::math

	PARTITIONS.FOR
C	Author: Lynn Yarbrough 
C	10 July 1989
C
C	For an explanation of the Algorithm, see
C	"A Course in  Number Theory" by H. E. Rose
C	Published by Oxford Science Publications, 1988
C	Chapter 11, Page 205
C
	REAL*16 p, partitions(0:200)
	INTEGER m, m1, m2, n, s
C======================================================================
	partitions(0) = 1
	DO n = 1, 200
	    s = +1
	    p = 0.0
	    m = 1
	    DO WHILE (n .GE. m*(3*m-1)/2)
		m1 = n - m*(3*m-1)/2
		m2 = n - m*(3*m+1)/2
		IF (m1 .GE. 0) p = p + s*partitions(m1)
		IF (m2 .GE. 0) p = p + s*partitions(m2)
		s = -s
		m = m + 1
	    END DO
	    partitions(n) = p
	END DO
	TYPE *, 200, partitions(200)
	END

     Does anyone have the series "approximation" of Ramanujan
     (and Hardy?) that is accurate enough that it rounds off to
     the correct answer [the number of partitions is always an
     integer, the approximation is off by less than 1/2]?
     
     Dan

    The best (most accurate) "closed form" expression I know of is due to H.
    Radamacher.

    It's in AMS55 in the section on combinatorics (and is fairly hairy, but
    would not be hard to program.)  Also, there are numerous references
    on the problem in that chapter.

    There are simpler expressions, but I can't quote them off the top
    of my head.  Radamacher's paper may make interesting reading if you
    can get the library to obtain a copy...

    - Jim

I have a copy of Hardy's monograph/lecture notes that contains the analysis.

I had originally thought of writing up the method in MAPLE until I found 
the algorithm in .1, which is much simpler (I never could quite wade through
all of H & R's method) and very fast. There is a curious footnote in
Hardy's notes to the effect that Lehmer has shown that Ramunajan's series
did not in fact converge (although it certainly appears to in the first few
terms)! 

This, by the way, is a problem that has quietly bugged me since my graduate 
school days in the mid-fifties! It felt really good to be able to display 
this algorithm.

Lynn 

     Asymptotic series often fail to converge.  For instance,
     Stirling's series for ln n!, the series in powers of 1/n.
     (Stirling's well known approximation is just the first
     term). If you take a fixed number of terms of the series, it
     is a more and more accurate relative approximation (i.e.,
     percentage-wise) as n -> oo.  But for any fixed n, the
     series diverges as the number of terms -> oo.  I suppose
     that for a particular n you keep adding up terms until the
     terms stop decreasing in absolute value and start
     increasing, to get the best approximation that the series
     can deliver.
     
     Dan

OK, now for an obscure problem: what's the smallest N for which 
partitions(n) contains all 10 digits?


Answer after the spoiler ...

   n		p(n)
  247    182973889854026.	then, shortly after that comes 
  250    230793554364681.

I checked these visually, and I might have missed an earlier one. - Lynn 

The FORTRAN (REAL*16) and MAPLE results agree up to 
	p(1142) = 5523864341942100491068450472029219.
After that, REAL*16 arithmetic is inadequate to represent the results 
correctly.

Lynn Yarbrough 

        What does the algorithm compute for the example in the
        700's in the paper you sent me?  I have it at home I'll
        post the result of the approximation.
        
        Dan

        The paper gives an approximation for p(721), using the
        first 21 terms of the series, of
        
        	161,061,755,750,279,477,635,534,762.0041
        
        Since p(721) is an integer, just truncate that.  Does
        that give the same result as Maple does?  I predict yes. :-)
        
        Dan

Both the MAPLE and FORTRAN codes duplicate Lehmer's result, 

	 p(721)=   161061755750279477635534762.

This particular number was calculated to disprove a conjecture about the 
divisibility of partition numbers by certain powers of divisors of the
arguments. Much of Ramanujan's original work - done long before computers 
existed - had to do with certain patterns in the divisors of p(n); e.g. the
above number is also 

                      2     3                            
               (2) (7)  (11)  (97) (12729652951228673767)

and some of the smaller divisors can be determined from the argument, 721, 
directly.

Conjecture:
	for all finite [integer] k,
	lim(N->inf) partitions(N+k)/partitions(N) = 1

This is not at all obvious, especially for large k, as the convergence even 
for k=1 is very slow. If it can be proved for k=1, the rest will follow, I
think, since if 
	partitions(M+1)/partitions(M) = 1+e(M), (e=epsilon)
then	partitions(M+k)/partitions(M) ~= (1+e)^k ~= 1+ke, etc.

Lynn Yarbrough 

    Re .11:
    
    I don't think so, Lynn.  It seems pretty easy to prove that:
    
    	Partitions(N+1) >= 2 * Partitions(N)
    
    Right?
    
    --  Barry

        re .11    >> Right?
        
        Wrong.  See reply .6 for a counterexample.
        
        Dan

>     <<< Note 1096.11 by CIVAGE::LYNN "Lynn Yarbrough @WNP DTN 427-5663" >>>
>                      -< ratio of partitions: conjecture >-
>
>Conjecture:
>	for all finite [integer] k,
>	lim(N->inf) partitions(N+k)/partitions(N) = 1

   This follows at once because the radius of convergence of the
   generating function for p(n) is 1, and the ratio test works for the
   series.

   - Jim

   Speaking of generating functions, here is a program that shows that
   you can approximate the coefficients of a power series generating
   function by numerically Fourier analyzing the function values around
   a circular contour with an FFT.

   Cauchy's residue theorem shows that the coefficients of a Laurent
   series can be picked off one by one with contour integrals; if the
   contour is circular, the integrals basically amount to calculating
   the Fourier coefficients of the "periodic" function values around
   the circle.

   Point sampling around the contour amounts to use of the trapezoidal
   rule; since the function is analytic on the circle, if it doesn't
   pass through any singularity the "aliasing" errors from the sampling
   involved are nicely bounded, and the error terms at the "endpoints"
   of the interval vanish since the function and all derivatives are
   periodic.  There is a tradeoff of roundoff error in the FFT, and
   aliasing error.  Bringing the radius nearer singularities lets in
   more higher harmonics, but too small a radius and roundoff error will
   contaminate the answer.  Using more terms in the FFT will of course
   be slower, and also adds a little roundoff noise as well.

   In VAX D float, this program gives the first 256 partitions correctly.

   By merely putting in other "generating functions", you can get other
   combinatorial functions; for example, the coefficients of Klein's
   modular J invariant from elliptic curve theory are really hairy to
   calculate, but the program got the first dozen or so correctly,
   as if by magic :-)

   By the way, if you have a multi-sheeted Riemann surface and evaluate
   the (multivalued) function around a circular contour that winds
   around the sheets and returns to the starting point after k turns,
   the method will give a Laurent series in powers of z^(1/k) - the
   fractional power disambiguates which sheet the function falls on.
   It is not necessary for the center of the series to be at a branch
   point for this to work - it can be done around any circular
   contour on a covering of the Riemann sphere.

   - Jim

	PROGRAM PN
C
C Compute unrestricted partition function p(n) via residue theorem and FFT
C
	IMPLICIT NONE

	INTEGER*4 LOGN, NPTS, NPN, IPASS, I, J, K, L
	COMPLEX*16 T, U, V, W, Z, Z0, A(1024), B(1024)
	REAL*8 PI, X, R

	DATA PI/3.14159265358979323846264338D0/
C
C Generating function for p(n)
C
	COMPLEX*16 F
C
C Number of points around contour, number of series coefficients
C
	LOGN = 10
	NPTS = 2**LOGN
	NPN = 2**(LOGN-2)

	Z0 = 0.0
C
C Adhoc contour radius
C
	R = 0.919
C
C Approximate contour integrals with trapezoidal rule using FFT
C
	A(1) = F(Z0+R)
	A(2) = F(Z0-R)
	J = 0
	DO I = 1,NPTS/2-1
	  V = R*DCMPLX(DCOS((2.0*PI*I)/NPTS), DSIN((2.0*PI*I)/NPTS))
	  K = NPTS
100	  CONTINUE
	  K = K/2
	  J = J.XOR.K
	  IF ((J.AND.K).EQ.0) GOTO 100
	  A(J+1) = F(Z0+V)
	  A(J+2) = F(Z0-V)
	ENDDO

	L = 1
	DO IPASS = 1,LOGN
	  U = 1.0
	  DO J = 1,L
	    DO I = J,NPTS,2*L
	      K = I+L
	      T = A(K)*U
	      A(K) = A(I)-T
	      A(I) = A(I)+T
	    ENDDO
	    U =  DCMPLX(DCOS((PI*J)/L), -DSIN((PI*J)/L))
	  ENDDO
	  L = L+L
	ENDDO
C
C Denormalize series coefficients
C
	X = NPTS
	DO I = 1,NPTS
	  B(I) = A(I)/X
	  X = X*R
	ENDDO
C
C Show the answer
C
	TYPE 101, (I-1, DREAL(B(I)), DREAL(A(I)), I=1,NPN)
101	FORMAT (1X, I10, F36.6, G28.17)

	STOP

	END

	COMPLEX*16 FUNCTION F(Z)
C
C Generating function for p(n) = 1/PROD(k>0) (1-z^k)
C
	COMPLEX*16 Z, ZN, T

	ZN = Z
	F = 1.0
100	CONTINUE
	T = F		
	F = F*(1.0-ZN)
	IF (CDABS(F-T).EQ.0.0) GOTO 200
	ZN = ZN*Z
	GOTO 100
200	CONTINUE

	F = 1.0/F

	RETURN

	END