Pastebin launched a little side project called VERYVIRAL.com, check it out ;-) Want more features on Pastebin? Sign Up, it's FREE!
Guest

Paritiotions into at most 60 parts

By: a guest on Jul 27th, 2011  |  syntax: None  |  size: 2.79 KB  |  views: 140  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  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 ====