/* -*- gp-script -*- */ \\% Conversion between power series and Lambert series \\ Author: Joerg Arndt \\ License: GPL version 3 or later \\ online at http://www.jjj.de/pari/ \\ version: 2014-October-16 (18:30) ser2lambert(t)= { /* Let t=[a1,a2,a3, ...], n=length(v), where t(x)=sum_{k=1}^{n}{a_k*x^k}; * Return L=[l1,l2,l3,...] so that (up to order n) * t(x)=\sum_{j=1}^{n}{l_j*x^j/(1-x^j)} */ my(n, L); n = length(t); L = vector(n); for (k=1, n, fordiv(k, d, L[k]+=moebius(k/d)*t[d]); ); \\ for (k=1, n, \\ greedy algorithm \\ tk = t[k]; \\ L[k] = tk; \\ \\ subtract tk * x^k/(1-x^k): \\ forstep(j=k, n, k, t[j] -= tk); \\ ); return( L ); } /* ----- */ lambert2ser(L)= { /* inverse of ser2lambert() */ my(n, t); n = length(L); t = sum(k=1, length(L), O('x^(n+1))+L[k]*'x^k/(1-'x^k) ); t = Vec(t); return( t ); } /* ----- */ ser2lambertplus(t)= { /* Let t=[a1,a2,a3, ...], n=length(v), where t(x)=sum_{k=1}^{n}{a_k*x^k}; * Return L=[l1,l2,l3,...] so that (up to order n) * t(x)=\sum_{j=1}^{n}{l_j*x^j/(1+x^j)} */ my(n, L, k4); n = length(t); L = vector(n); for (k=1, n, \\ greedy algorithm tk = t[k]; L[k] = tk; \\ subtract tk * x^k/(1+x^k): forstep(j=k, n, 2*k, t[j] -= tk); forstep(j=k+k, n, 2*k, t[j] += tk); ); return( L ); } /* ----- */ lambertplus2ser(L)= { /* inverse of ser2lambertplus() */ my(n, t); n = length(L); t = sum(k=1, length(L), O('x^(n+1))+L[k]*'x^k/(1+'x^k) ); t = Vec(t); return( t ); } /* ----- */ \\ ==== end of file ====