/* -*- gp-script -*- */ \\% Conversion of products to continued fractions. \\ Author: Joerg Arndt \\ License: GPL version 3 or later \\ online at http://www.jjj.de/pari/ \\ version: 2014-October-16 (18:32) \\ Let P := \prod(k=0, n, (1+Y_k) \\ then CF(P) = CF(A, B) where \\ A and B are vectors returned by \\ cfproda(V,n) and cfprodb(V, n), respectively. \\ V is the vector of the values Y_k. cfproda(yv, n)= { my( y2, y3, d ); if ( n<=1, return(1) ); y3 = yv[n]; y2 = yv[n-1]; return( (y2+y3*(1+y2)) ); } /* ----- */ cfprodb(yv, n)= { my( y1, y2, y3 ); if (0==n, return(1) ); \\ unused if (1==n, return(+yv[1+0]) ); y3 = yv[n]; y2 = yv[n-1]; y1 = if ( n==2, 1, yv[n-2] ); return( -y1*y3*(1+y2) ); } /* ----- */ cfprod(yv)= { /* Example: Y(k)=eval(Str("Y"k)); yv=vector(n,j,Y(j-1)) ==> [Y0, Y1, Y2] t=cfprod(yv) ==> [[1, 1, (Y1 + 1)*Y0 + Y1, (Y2 + 1)*Y1 + Y2], [1, Y0, -Y1*Y0 - Y1, (-Y2*Y1 - Y2)*Y0]] */ my(n, a, b); n = length(yv); n += 1; \\ n+1 terms in continued fraction a = vector(n); b = vector(n); for (k=0, n-1, a[k+1] = cfproda(yv, k); b[k+1] = cfprodb(yv, k); ); return( [a, b] ); } /* ----- */ \\ ==== end of file ====