/* -*- gp-script -*- */ \\% Pade approximant via EGCD \\ Author: Joerg Arndt \\ License: GPL version 3 or later \\ online at http://www.jjj.de/pari/ \\ version: 2014-October-16 (18:29) pade(s, d)= /* Compute Pade approximant A/B (of the power series S) such that deg(A)=deg(S)-d. Must have: d < deg(S); S must have a nonzero constant term. */ { my(n, t, q); s = truncate(s); \\ remove O(x^(n+1)) term if present n = poldegree(s); u = [1, 0, x^(n+1)]; v = [0, 1, s]; t = [0, 0, 0]; \\ print1(" P[",n, ",",0, "] = (", v[3],") / (",v[2],")"); \\ == s / 1 \\ print(); \\ d = min(n, d); for ( j=1, d, q = (u[3] \ v[3]); \\ division without remainder t = u - v*q; u = v; v = t; \\ print1(" P[",n-j, ",",j, "] = (", v[3],") / (",v[2],")"); print(); ); return( v[3]/v[2] ); } /* ----- */ \\ ==== end of file ====