• API
• FAQ
• Tools
• Trends
• Archive
SHARE
TWEET

# Paritiotions into at most 60 parts

a guest Jul 27th, 2011 168 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1.
2. \\% Basic manipulations of vectors.
3. \\ Author: Ralf Stephan, (reformatted by Joerg Arndt)
4. \\ License: GPL version 3 or later
5. \\ online at http://www.jjj.de/pari/
6. \\ version: 2011-April-08 (19:18)
7.
8.
9. \\ gp/pari code for Doron's  challenge at http://mathoverflow.net/questions/71092/how-many-integer-partitions-of-a-googol-10100-into-at-most-60-parts/71132#71132
10. \\Needs http://www.jjj.de/pari/fastrec.gpi
11.
12. ggf(v)=
13. { /* Guess generating function for supplied vector of integers */
14. /* Example: ggf(vector(10,j,fibonacci(j)^2)) ==> (-x + 1)/(x^3 - 2*x^2 - 2*x + 1) */
15.     local(l,m,p,q,B);
16.     l = length(v);
17.     B = floor(l/2);
18.     if ( B<4,  return(0) );
19.     m = matrix(B, B, x, y, v[x-y+B+1]);
20.     q = qflll(m, 4)[1];
21.     if ( length(q)==0,  return(0) );
22. \\    p = sum(k=1, B, x^(k-1)*q[k,1]);
23.     p = sum(k=1, B, 'x^(k-1)*q[k,1]);
24. \\    q = Pol( Pol( vector(l, n, v[l-n+1]) )*p + O(x^(B+1)) );
25.     q = Pol( Pol( vector(l, n, v[l-n+1]) )*p + O('x^(B+1)) );
26.     if ( polcoeff(p,0)<0,  q=-q; p=-p);
27.     q = q/p;
28. \\    p = Ser(q+O(x^(l+1)));
29.     p = Ser(q+O('x^(l+1)));
30.     for (m=1, l, if( polcoeff(p,m-1)!=v[m],  return(0)) );
31.     return(q);
32. } /* ------- */
33.
34.
35. \\ aux:
36. strsgn(a)=if(a>=0, Str("+" a), Str(a));
37. \\prtpol(c)=
38. \\{ /* print coefficients with ascending degrees */
39. \\    print1("(");
40. \\    for (d=(poldegree(c)-...)
41. \\    print1(")");
42. \\} /* ----- */
43. \\
44. \\prtgf(v)=
45. \\{ /* Print Generating function */
46. \\    local(gf, t);
47. \\    gf = ggf(v);
48. \\    if ( gf==0, return(0) );
49. \\    t = numerator(gf);
50. \\    print1
51. \\    return(1);
52. \\} /* ----- */
53.
54.
55. prtrec(v)=
56. { /* Print recurrence relation */
57.     local(c,d,c0,t);
58.     c=ggf(v);
59.     if ( c==0, return(0) );
60.     c = denominator(c);
61. \\    print(c);
62. \\    c = polrecip(c);
63. \\    print(c);
64.     d = poldegree(c);
65.     c0 = polcoeff(c,0);
66.     print1("a(n) =");
67.     for (k=1, d,
68.         t=polcoeff(c,k)/c0;
69.         if (t, print1(" "); print1( strsgn( -t ), "*a(n-",k,")" ) );
70.     );
71.     print();
72.     return(1);
73. } /* ----- */
74.
75. prtrec2(c)=
76. { /* Print recurrence relation */
77.     local(d,c0,t,r,t2);
78.     \\c=ggf(v);
79.     if ( c==0, return(0) );
80.     c = denominator(c);
81. \\    print(c);
82. \\    c = polrecip(c);
83. \\    print(c);
84.     d = poldegree(c);
85.     r=[];\\listcreate(d);
86.     c0 = polcoeff(c,0);
87.     print1("a(n) =");
88.     forstep (k=1, d,1,
89.         t=polcoeff(c,k)/c0;
90.         t2=-polcoeff(c,k)/c0;
91.         r=concat(r,t2);
92.         if (t, print1(" "); print1( strsgn( -t ), "*aa(n-",k,")" ) );
93.     );
94.     print();
95.     return(r);
96. } /* ----- */
97.
98. k4=60;
99. d=prod(n=1,k4,(1-x^n));
100. {aa(n)=return(polcoeff(1/d+O(x^(n+1)),n ));};
101. zz=prtrec2(1/d);
102. a0=[];\\listcreate(poldegree(d));
103. for(n=0,poldegree(d)-1,a0=concat(a0,aa(n)); );
104. {getd()=d;}
105. {geta0()=a0;}
106. {getzz()=zz;}
107.
108. X=frec(a0,zz,10^100)[1];
109.
110. \\ ==== end of file ====
RAW Paste Data
Top