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 ====