Formal power series

Formal power series
You are encouraged to solve this task according to the task description, using any language you may know.

A power series is an infinite sum of the form

$a_{0}+a_{1}\cdot x+a_{2}\cdot x^{2}+a_{3}\cdot x^{3}+\cdots$

The a_i are called the coefficients of the series. Such sums can be added, multiplied etc., where the new coefficients of the powers of x are calculated according to the usual rules.

If one is not interested in evaluating such a series for particular values of x, or in other words, if convergence doesn't play a role, then such a collection of coefficients is called formal power series. It can be treated like a new kind of number.

Task: Implement formal power series as a numeric type. Operations should at least include addition, multiplication, division and additionally non-numeric operations like differentiation and integration (with an integration constant of zero). Take care that your implementation deals with the potentially infinite number of coefficients.

As an example, define the power series of sine and cosine in terms of each other using integration, as in

$\sin x=\int _{0}^{x}\cos t\,dt$

$\cos x=1-\int _{0}^{x}\sin t\,dt$

Goals: Demonstrate how the language handles new numeric types and delayed (or lazy) evaluation.

Ada

The Taylor series package is generic to be instantiated with any rational type implementation provided by the task Rational Arithmetic: <lang ada>with Generic_Rational;

generic

  with package Rational_Numbers is new Generic_Rational (<>);

package Generic_Taylor_Series is

  use Rational_Numbers;
  type Taylor_Series is array (Natural range <>) of Rational;
  
  function "+" (A : Taylor_Series) return Taylor_Series;
  function "-" (A : Taylor_Series) return Taylor_Series;

  function "+" (A, B : Taylor_Series) return Taylor_Series;
  function "-" (A, B : Taylor_Series) return Taylor_Series;
  function "*" (A, B : Taylor_Series) return Taylor_Series;

  function Integral (A : Taylor_Series) return Taylor_Series;
  function Differential (A : Taylor_Series) return Taylor_Series;

  function Value (A : Taylor_Series; X : Rational) return Rational;

  Zero : constant Taylor_Series := (0 => Rational_Numbers.Zero);
  One  : constant Taylor_Series := (0 => Rational_Numbers.One);

end Generic_Taylor_Series;</lang> The package implementation: <lang ada>package body Generic_Taylor_Series is

  function Normalize (A : Taylor_Series) return Taylor_Series is
  begin
     for Power in reverse A'Range loop
        if A (Power) /= 0 then
           return A (0..Power);
        end if;
     end loop;
     return Zero;
  end Normalize;

  function "+" (A : Taylor_Series) return Taylor_Series is
  begin
     return A;
  end "+";

  function "-" (A : Taylor_Series) return Taylor_Series is
     Result : Taylor_Series (A'Range);
  begin
     for Power in A'Range loop
        Result (Power) := -A (Power);
     end loop;
     return Result;
  end "-";

  function "+" (A, B : Taylor_Series) return Taylor_Series is
  begin
     if A'Last > B'Last then
        return B + A;
     else
        declare
           Result : Taylor_Series (0..B'Last);
        begin
           for Power in A'Range loop
              Result (Power) := A (Power) + B (Power);
           end loop;
           for Power in A'Last + 1..B'Last loop
              Result (Power) := B (Power);
           end loop;
           return Normalize (Result);
        end;
     end if;
  end "+";

  function "-" (A, B : Taylor_Series) return Taylor_Series is
  begin
     return A + (-B);
  end "-";

  function "*" (A, B : Taylor_Series) return Taylor_Series is
     Result : Taylor_Series (0..A'Last + B'Last);
  begin
     for I in A'Range loop
        for J in B'Range loop
           Result (I + J) := A (I) * B (J);
        end loop;
     end loop;
     return Normalize (Result);
  end "*";

  function Integral (A : Taylor_Series) return Taylor_Series is
  begin
     if A = Zero then
        return Zero;
     else
        declare
           Result : Taylor_Series (0..A'Last + 1);
        begin
           for Power in A'Range loop
              Result (Power + 1) := A (Power) / Number (Power + 1);
           end loop;
           Result (0) := Rational_Numbers.Zero;
           return Result;
        end;
     end if;
  end Integral;

  function Differential (A : Taylor_Series) return Taylor_Series is
  begin
     if A'Length = 1 then
        return Zero;
     else
        declare
           Result : Taylor_Series (0..A'Last - 1);
        begin
           for Power in Result'Range loop
              Result (Power) := A (Power + 1) * Number (Power);
           end loop;
           return Result;
        end;
     end if;         
  end Differential;

  function Value (A : Taylor_Series; X : Rational) return Rational is
     Sum : Rational := A (A'Last);
  begin
     for Power in reverse 0..A'Last - 1 loop
        Sum := Sum * X + A (Power);
     end loop;
     return Sum;
  end Value;

end Generic_Taylor_Series;</lang> The procedure Normalize is used to truncate the series when the coefficients are zero. The summation of a series (function Value) uses Horner scheme.

Test task

<lang ada>with Ada.Text_IO; use Ada.Text_IO;

with Generic_Taylor_Series; with Generic_Rational;

procedure Test_Taylor_Series is

  package Integer_Rationals is new Generic_Rational (Integer);
  package Integer_Taylor_Series is new Generic_Taylor_Series (Integer_Rationals);
  use Integer_Taylor_Series;
     -- Procedure to print a series
  procedure Put (A : Taylor_Series) is
     use Integer_Rationals;
     procedure Put (A : Rational) is
     begin
        if Numerator (A) = 1 then
           Put (" 1");
        else
           Put (Integer'Image (Numerator (A)));
        end if;
        if Denominator (A) /= 1 then
           Put (" /");
           Put (Integer'Image (Denominator (A)));
        end if;
     end Put;
  begin
     if A (0) /= 0 then
        Put (A (0));
     end if;
     for Power in 1..A'Last loop
        if A (Power) > 0 then
           Put (" +");
           Put (A (Power));
           Put (" X **" & Integer'Image (Power));
        elsif A (Power) < 0 then
           Put (" -");
           Put (abs A (Power));
           Put (" X **" & Integer'Image (Power));
        end if;
     end loop;
  end Put;
     -- Cosine generator
  function Cos (N : Natural) return Taylor_Series is
  begin
     if N = 0 then
        return One;
     else
        return One - Integral (Integral (Cos (N - 1)));
     end if;
  end Cos;

begin

  Put ("Cos ="); Put (Cos (5)); Put_Line (" ...");
  Put ("Sin ="); Put (Integral (Cos (5))); Put_Line (" ...");

end Test_Taylor_Series;</lang> Sample output:

Cos = 1 - 1 / 2 X ** 2 + 1 / 24 X ** 4 - 1 / 720 X ** 6 + 1 / 40320 X ** 8 - 1 / 3628800 X ** 10 ...
Sin = + 1 X ** 1 - 1 / 6 X ** 3 + 1 / 120 X ** 5 - 1 / 5040 X ** 7 + 1 / 362880X ** 9 - 1 / 39916800 X ** 11 ...

Clojure

This version takes advantage of the laziness of most of Clojure's sequence functions, including map, for, take-while, concat, and drop. A formal power series (FPS) is represented as a sequence of coefficients; for example, [1 2 1] represents 1 + 2*x + 1*x*x.

First addition and subtraction. Note that most of the complication arises in allowing for finite and infinite FPSs; if only infinite power series were at issue, the function (defn ips+ [ips0 ips1] (map + ips0 ips1)) would suffice. <lang lisp> (defn ps+ [ps0 ps1]

 (letfn [(+zs   [ps] (concat ps (repeat :z)))
         (notz? [a] (not= :z a))
         (nval  [a] (if (notz? a) a 0))
         (z+    [a0 a1] (if (= :z a0 a1) :z (+ (nval a0) (nval a1))))]
   (take-while notz? (map z+ (+zs ps0) (+zs ps1)))))

(defn ps- [ps0 ps1] (ps+ ps0 (map - ps1))) </lang> Multiplication next; again most of the complication is dealing with both finite and infinite FPS. This function explicitly uses the standard function lazy-seq to define the product sequence. <lang lisp> (defn indexed [ps] (map vector (iterate inc 0) ps))

(defn ps*

 ([ps0 ps1] (ps* [0] ps0 ps1))
 ([[a0 & resta] [p0 & rest0] [p1 & rest1 :as ps1]]
   (lazy-seq
     (cons
       (+ a0 (* p0 p1))
       (let [mrest1 (if (or (nil? rest1) (zero? p0)) nil, (map #(* p0 %) rest1))
             accum  (cond (nil? resta) mrest1, (nil? mrest1) resta, :else (ps+ resta mrest1))]
         (if (nil? rest0) accum, (ps* (or accum [0]) rest0 ps1)))))))

</lang> As with most of the other examples on this page, there's no definition for division. Mathematically, FPS is a commutative ring (in fact a Euclidean domain), but not a field: the set of FPSs is not closed under division.

Now we can define integration and differentiation: <lang> (defn differentiate [ps]

 (drop 1 (for [[n a] (indexed ps)] (* n a))))

(defn integrate [ps]

 (cons 0 (for [[n a] (indexed ps)] (/ a (inc n)))))

</lang> Some examples of using these functions; in each case a println call (which forces the lazy sequence) is followed by a comment showing the output. <lang lisp> (println (ps+ [1 2] [3 4 5]))

(4 6 5)

(println (ps* [1 2] [3 4 5]))

(3 10 13 10)

</lang> And some examples using infinite FPSs. First define the sequence of factorials (facts), then define sin and cos Taylor series. <lang lisp> (def nfacts (iterate (fn f n [(* f n) (inc n)]) [1 1])) (def facts (map first nfacts))

(def sin (map / (cycle [0 1 0 -1]) facts)) (def cos (map / (cycle [1 0 -1 0]) facts))

(println (take 10 sin))

(0 1 0 -1/6 0 1/120 0 -1/5040 0 1/362880)

(println (take 10 (integrate cos)))

(0 1 0 -1/6 0 1/120 0 -1/5040 0 1/362880)

(println (take 20 (ps+ (ps* sin sin) (ps* cos cos))))

(1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0)

</lang> By using letfn, which supports defining mutually recursive functions, we can define the sin and cos power series directly in terms of integrals of the other series: <lang lisp> (letfn [(fsin [] (lazy-seq (integrate (fcos))))

       (fcos [] (ps- [1] (integrate (fsin))))]
 (def sinx (fsin))
 (def cosx (fcos)))

(println (take 10 sinx))

(0 1 0 -1/6 0 1/120 0 -1/5040 0 1/362880)

</lang>

Common Lisp

Common Lisp isn't lazy, and doesn't define the arithmetic operators as generic functions. As such, this implementation defines lazy primitives (delay and force), and a lazy list built on top of them (lons, lar, ldr). This implementation also defines a package #:formal-power-series which uses all the symbols of the "COMMON-LISP" package except for +, -, *, and /, which are shadowed. Shadowing these symbols allows for definitions of generic functions which can be specialized for the power series, can default to the normal CL arithmetic operations for other kinds of objects, and can do "the right thing" for mixes thereof.

<lang lisp>(defpackage #:formal-power-series

 (:nicknames #:fps)
 (:use "COMMON-LISP")
 (:shadow
  #:+ #:- #:* #:/))

(in-package #:formal-power-series)</lang>

Lazy Primitives

<lang lisp>(defstruct promise

 thunk value)

(defmacro delay (form)

 `(make-promise :thunk #'(lambda () ,form)))

(defun force (object)

 (cond
  ((not (promise-p object))
   object)
  ((null (promise-thunk object))
   (promise-value object))
  (t (let ((val (funcall (promise-thunk object))))
       (setf (promise-thunk object) nil
             (promise-value object) val)))))</lang>

Lazy Lists

<lang lisp>(defstruct lons

 lar
 ldr)

(defun lar (lons)

 (lons-lar lons))

(defun ldr (lons)

 (if (not (promise-p (lons-ldr lons)))
   (lons-ldr lons)
   (setf (lons-ldr lons)
         (force (lons-ldr lons)))))

(defmacro lons (lar ldr)

 `(make-lons :lar ,lar :ldr (delay ,ldr)))</lang>

A few utilities to make working with lazy lists easier.

<lang lisp>(defun invoke-with-lons (function lons)

 (funcall function (lar lons) (ldr lons)))

(defmacro with-lons ((lar ldr) lons &body body)

 `(invoke-with-lons #'(lambda (,lar ,ldr) ,@body) ,lons))

(defun maplar (function llist &rest llists)

 (let ((llists (list* llist llists)))
   (if (some 'null llists) nil
     (lons (apply function (cl:mapcar 'lar llists))
           (apply 'maplar function (cl:mapcar 'ldr llists))))))

(defun take (n llist)

 (if (zerop n) '()
   (lons (lar llist)
         (take (1- n) (ldr llist)))))

(defun force-list (llist)

 (do ((fl '() (cons (lar l) fl))
      (l llist (ldr l)))
     ((null l) (nreverse fl))))

(defun repeat (x)

 (lons x (repeat x)))

(defun up-from (n)

 (lons n (up-from (1+ n))))</lang>

Formal Power Series

The mathematical operations here are translations of the Haskell code, but we specialize the operations in various ways so that behavior for normal numeric operations is preserved.

<lang lisp>(defstruct (series (:constructor series (coeffs)) (:conc-name))

 coeffs)

(defgeneric negate (f)

 (:method (f)
  (cl:- f))
 (:method ((f series))
  (series (maplar 'negate (coeffs f)))))

(defgeneric + (f g)

 (:method (f g)
  (cl:+ f g))
 (:method (f (g series))
  (series (lons (+ f (lar (coeffs g))) (ldr (coeffs g)))))
 (:method ((f series) g)
  (+ g f))
 (:method ((f series) (g series))
  (series (maplar '+ (coeffs f) (coeffs g)))))

(defun - (f g)

 (+ f (negate g)))

(defun series-* (f g)

 (with-lons (f ft) (coeffs f)
   (with-lons (g gt) (coeffs g)
     (series (lons (* f g)
                   (coeffs (+ (* (series ft)
                                 (series gt))
                              (* f (series gt)))))))))

(defgeneric * (f g)

 (:method (f g)
  (cl:* f g))
 (:method ((f series) g)
  (series (maplar #'(lambda (x) (* x g)) (coeffs f))))
 (:method (f (g series))
  (* g f))
 (:method ((f series) (g series))
  (series-* f g)))

(defun series-/ (f g)

 (with-lons (f ft) (coeffs f)
   (with-lons (g gt) (coeffs g)
     (let ((qs nil))
       (setf qs (lons (/ f g)
                      (maplar #'(lambda (x) (/ x g))
                              (coeffs (- (series ft)
                                         (* (series qs)
                                            (series gt)))))))))))

(defgeneric / (f g)

 (:method (f g)
  (cl:/ f g))
 (:method ((f series) g)
  (series (maplar #'(lambda (x) (/ x g)) (coeffs f))))
 (:method (f (g series))
  (/ (series (lons f (repeat 0))) g))
 (:method ((f series) (g series))
  (series-/ f g)))

(defun int (f)

 (series (lons 0 (maplar '/ (coeffs (force f)) (up-from 1)))))

(defun diff (f)

 (series (maplar '* (ldr (coeffs f)) (up-from 1))))</lang>

Example

<lang lisp>(defparameter *sinx*

 (locally (declare (special *cosx*))
   (delay (int (force *cosx*)))))

(defparameter *cosx*

 (delay (- 1 (int *sinx*))))</lang>

FPS > (force-list (take 10 (coeffs (force *sinx*))))
(0 1 0 -1/6 0 1/120 0 -1/5040 0 1/362880)

FPS > (force-list (take 10 (coeffs (force *cosx*))))
(1 0 -1/2 0 1/24 0 -1/720 0 1/40320 0)

FPS > (force-list (take 10 (coeffs (+ (force *cosx*) (force *sinx*)))))
(1 1 -1/2 -1/6 1/24 1/120 -1/720 -1/5040 1/40320 1/362880)

FPS > (force-list (take 10 (coeffs (- (+ 2 (force *cosx*))
                                      (* 3 (force *sinx*))))))
(3 -3 -1/2 1/2 1/24 -1/40 -1/720 1/1680 1/40320 -1/120960)

D

Works with: D version 2.007

<lang d>module fps; import q ; // a module for Rational Number import std.format : sformat = format; import std.math, std.stdio ; interface Gene(U) { U coef(int n) ; } class Term(U) { U[] cache ; Gene!(U) gene ; this(Gene!(U) g) { gene = g ; } U opIndex(int n) { if(n < 0) return cast(U)0 ; if(n >= cache.length) for(int i = cache.length ; i <= n ; i++) cache ~= gene.coef(i) ; return cache[n] ; } } struct FPS(U) { alias Term!(U) UT ; alias Gene!(U) UG ; alias FPS!(U) UF ; static int DispTerm = 12 ; static string XVar = "x" ; UT term ; static UF opCall(UT t) { UF f ; f.term = t ; return f ;} static UF opCall(U[] polynomial) { UF f ; f.term = new UT( new class() UG { U coef(int n) { if (n < 0 || n >= polynomial.length ) return cast(U)0 ; return polynomial[n] ; };}) ; return f ; } U inverseCoef(int n) { U[] res = new U[n + 1] ; res[0] = 1/term[0] ; for(int i = 1 ; i <= n ; i++) { res[i] = cast(U)0 ; for(int j = 0 ; j < i ; j++) { res[i] += term[i - j] * res[j] ; } res[i] = -res[0] * res[i] ; } return res[n] ; } UF opAdd(UF rhs) { return UF(new UT(new class() UG { U coef(int n) { return term[n] + rhs.term[n] ; }}));} UF opSub(UF rhs) { return UF(new UT(new class() UG { U coef(int n) { return term[n] - rhs.term[n] ; }}));} UF opMul(UF rhs) { return UF(new UT(new class() UG { U coef(int n) { U res = cast(U) 0 ; for(int i = 0 ; i <=n ; i++) res += term[i]*rhs.term[n -i] ; // Cauchy product with rhs return res ; }}));} UF opDiv(UF rhs) { return UF(new UT(new class() UG { U coef(int n) { U res = cast(U) 0 ; for(int i = 0 ; i <=n ; i++) res += term[i]*rhs.inverseCoef(n -i) ; // Cauchy product with inverse of rhs return res ; }}));} UF Diff() { return UF(new UT(new class() UG { U coef(int n) { return term[n + 1] * (n + 1) ; }}));} UF Intg() { return UF(new UT(new class() UG { U coef(int n) { if (n == 0) return cast(U)0 ; return term[n-1]/n ; }}));} string toString(int dpTerm = -1) { string s ; if(dpTerm <= 0) dpTerm = DispTerm ; U c = term[0] ; if(c != cast(U)0) s ~= sformat("%s", c) ; for(int i = 1 ; i < dpTerm ; i++) if((c = term[i]) != cast(U)0) { string t ; if(c > 0 && s.length > 0) t = "+" ; if (c == cast(U)1) t ~= XVar ; else if ( c == cast(U)-1) t ~= sformat("-%s", XVar) ; else t ~= sformat("%s%s", c, XVar) ; if(i > 1) t ~= sformat("%s", i) ; s ~= t ; } if(s.length == 0) s = "0" ; s ~= "+..." ; return s ; } } alias Q!(long) QT ; alias FPS!(QT) RF ; void main() { QT.DenoLimit = 0 ; RF COS  ; RF SIN = COS.Intg ; COS = RF([QT(1,1)]) - SIN.Intg ; writefln("SIN(x) = ",SIN) ; writefln("COS(x) = ",COS) ; }</lang>

Output:

SIN(x) = x-1/6x3+1/120x5-1/5040x7+1/362880x9-1/39916800x11+...
COS(x) = 1-1/2x2+1/24x4-1/720x6+1/40320x8-1/3628800x10+...

The q module. Because of the limited precision of long, the denominator will quickly overflow when multiplied, but it is sufficient to show first few coefficients in this task.

<lang d>module q; import std.format : sformat = format; T gcd(T)(T u, T v) { // binary gcd algorithm if(u < 0) u = -u ; if(v < 0) v = -v ; if(u == 0 || v == 0) return u | v ; int shift = 0 ; for(; 0 == ((u|v) & 1) ; shift++) { u >>>= 1 ; v >>>= 1 ; } while((u & 1) == 0) u >>>= 1 ; while (v != 0) { while((v & 1) == 0) v >>>= 1 ; if (u <= v) v -= u ; else { v = u - v ; u -= v ;} v >>>= 1 ; } return u << shift ; } struct Q(T) { static T DenoLimit = int.max ; alias Q!(T) QX ; T num, den ; static QX opCall(long n, long d) { QX q ; with(q) { num = n ; den = d; reduce ; } return q ; } static QX opCall(real r) { QX q ; T multiple = 1_000_000 ; with(q) { num = cast(T)(r*multiple)  ; den = multiple; reduce ; } return q ; } void reduce() { T g = gcd(num, den) ; if(g == 0 || den == 0) throw new Exception(sformat("Zero denominator:%s/%s gcd:%s", num, den, g )) ; if(num == 0) { den = 1 ; return ; } num /= g ; den /= g ; if(den < 0) { num = -num ; den = -den ;} // keep den positive if(DenoLimit != 0 && den > DenoLimit) { // prevent denominator overflow, buggy? den /= 2; num /= 2; reduce ; } } T toBaseType() {T result = num /den ; return result ; } real toReal() { real result = num /cast(real)den ; return result ; } real opCast() { return toReal ; } T opCmp(QX rhs) { return num*rhs.den - den*rhs.num ; } T opCmp(int rhs) { return num - den*rhs ; } T opCmp(real rhs) { return cast(T)(num - den*rhs) ; } bool opEquals(QX rhs) { return opCmp(rhs) == 0 ; } bool opEquals(int rhs) { return opCmp(rhs) == 0 ; } bool opEquals(real rhs) { return opCmp(rhs) == 0 ; } QX opPostInc() { num += den ; reduce ; return *this ; } QX opPostDec() { num -= den ; reduce ; return *this ; } QX opAddAssign(QX rhs) { num = num*rhs.den + rhs.num*den ; den *= rhs.den ; reduce ; return *this ;} QX opSubAssign(QX rhs) { num = num*rhs.den - rhs.num*den ; den *= rhs.den ; reduce ; return *this ;} QX opMulAssign(QX rhs) { num *= rhs.num ; den *= rhs.den ; reduce ; return *this ;} QX opDivAssign(QX rhs) { num *= rhs.den ; den *= rhs.num ; reduce ; return *this ;} QX opAddAssign(int rhs) { num += rhs*den ; reduce ; return *this ; } QX opSubAssign(int rhs) { num -= rhs*den ; reduce ; return *this ; } QX opMulAssign(int rhs) { num *= rhs ; reduce ; return *this ; } QX opDivAssign(int rhs) { den *= rhs ; reduce ; return *this ; } QX opNeg() { return QX(-num, den) ;} QX opPos() { return QX(num, den) ;} QX opAdd(QX rhs) { return QX(num*rhs.den + rhs.num*den, den*rhs.den);} QX opSub(QX rhs) { return QX(num*rhs.den - rhs.num*den, den*rhs.den);} QX opMul(QX rhs) { return QX(num*rhs.num, den*rhs.den);} QX opDiv(QX rhs) { return QX(num*rhs.den, den*rhs.num);} QX opAdd(int rhs) { return QX(num + rhs*den, den);} QX opSub(int rhs) { return QX(num - rhs*den, den);} QX opSub_r(int lhs) { return QX(lhs*den - num, den);} QX opMul(int rhs) { return QX(num*rhs, den);} QX opDiv(int rhs) { return QX(num, den*rhs);} QX opDiv_r(int lhs) { return QX(lhs*den, num);} string toString() { if(den == 1) return sformat("%s", num) ; return sformat("%s/%s", num, den) ; } }</lang>

Haskell

It's simpler to assume we are always dealing with an infinite list of coefficients. Mathematically, a finite power series can be generalized to an infinite power series with trailing zeros. <lang haskell>newtype Series a = S { coeffs :: [a] } deriving (Eq, Show) -- Invariant: coeffs must be an infinite list

instance Num a => Num (Series a) where

 fromInteger n = S $ fromInteger n : repeat 0
 negate (S fs) = S $ map negate fs
 S fs + S gs   = S $ zipWith (+) fs gs
 S (f:ft) * S gs@(g:gt) = S $ f*g : coeffs (S ft * S gs + S (map (f*) gt))

instance Fractional a => Fractional (Series a) where

 fromRational n = S $ fromRational n : repeat 0
 S (f:ft) / S (g:gt) = S qs where qs = f/g : map (/g) (coeffs (S ft - S qs * S gt))

-- utility function to convert from a finite polynomial fromFiniteList xs = S (xs ++ repeat 0)

int (S fs) = S $ 0 : zipWith (/) fs [1..]

diff (S (_:ft)) = S $ zipWith (*) ft [1..]

sinx,cosx :: Series Rational sinx = int cosx cosx = 1 - int sinx

fiboS = 1 / fromFiniteList [1,-1,-1]</lang>

Output:

*Main> take 11 $ coeffs sinx
[0 % 1,1 % 1,0 % 1,(-1) % 6,0 % 1,1 % 120,0 % 1,(-1) % 5040,0 % 1,1 % 362880,0 % 1]

*Main> take 11 $ coeffs cosx
[1 % 1,0 % 1,(-1) % 2,0 % 1,1 % 24,0 % 1,(-1) % 720,0 % 1,1 % 40320,0 % 1,(-1) % 3628800]

*Main> take 11 $ coeffs $ sinx / cosx -- tangent
[0 % 1,1 % 1,0 % 1,1 % 3,0 % 1,2 % 15,0 % 1,17 % 315,0 % 1,62 % 2835,0 % 1]

*Main> take 11 $ map truncate $ coeffs $ fiboS  -- some fibonaccis
[1,1,2,3,5,8,13,21,34,55,89]

*Main> take 11 $ coeffs $ fromFiniteList [1,5] * fromFiniteList [2,3,7] -- multiplying polynomials
[2,13,22,35,0,0,0,0,0,0,0]

Java

Works with: Java version 1.5+

Copied from the D example. Java has no generic numeric interface, and has no templates, so we cannot make it work for all numeric types at the same time. Because the Java library does not come with a Fraction class (and I am too lazy to implement one, although there is one in Apache Commons Math), here I will just hard-code it to use doubles (64-bit floating-point numbers) instead of fractions. It won't be as pretty, but it can be changed to Fraction or any other type trivially by substituting the types. <lang java5>import java.util.*;

interface Gene {

   double coef(int n);

}

class Term {

   private final List<Double> cache = new ArrayList<Double>();
   private final Gene gene;

   public Term(Gene g) { gene = g; }

   public double get(int n) {
       if (n < 0)
           return 0;
       else if (n >= cache.size())
           for (int i = cache.size(); i <= n; i++)
               cache.add(gene.coef(i));
       return cache.get(n);
   }

}

public class FormalPS {

   private static final int DISP_TERM = 12;
   private static final String X_VAR = "x";
   private Term term;

   public FormalPS() { }
   public void copyFrom(FormalPS foo) {
       term = foo.term;
   }

   public FormalPS(Term t) {
       term = t;
   }

   public FormalPS(final double[] polynomial) {
       this(new Term(new Gene() {
               public double coef(int n) {
                   if (n < 0 || n >= polynomial.length)
                       return 0;
                   else
                       return polynomial[n];
               }
           }));
   }

   public double inverseCoef(int n) {
       double[] res = new double[n + 1];
       res[0] = 1 / term.get(0);
       for (int i = 1; i <= n; i++) {
           res[i] = 0;
           for (int j = 0; j < i; j++)
               res[i] += term.get(i-j) * res[j];
           res[i] *= -res[0];
       }
       return res[n];
   }

   public FormalPS add(final FormalPS rhs) {
       return new FormalPS(new Term(new Gene() {
               public double coef(int n) {
                   return term.get(n) + rhs.term.get(n);
               }
           }));
   }

   public FormalPS sub(final FormalPS rhs) {
       return new FormalPS(new Term(new Gene() {
               public double coef(int n) {
                   return term.get(n) - rhs.term.get(n);
               }
           }));
   }

   public FormalPS mul(final FormalPS rhs) {
       return new FormalPS(new Term(new Gene() {
               public double coef(int n) {
                   double res = 0;
                   for (int i = 0; i <= n; i++)
                       res += term.get(i) * rhs.term.get(n-i);
                   return res;
               }
           }));
   }

   public FormalPS div(final FormalPS rhs) {
       return new FormalPS(new Term(new Gene() {
               public double coef(int n) {
                   double res = 0;
                   for (int i = 0; i <= n; i++)
                       res += term.get(i) * rhs.term.get(n-i);
                   return res;
               }
           }));
   }

   public FormalPS diff() {
       return new FormalPS(new Term(new Gene() {
               public double coef(int n) {
                   return term.get(n+1) * (n+1);
               }
           }));
   }

   public FormalPS intg() {
       return new FormalPS(new Term(new Gene() {
               public double coef(int n) {
                   if (n == 0)
                       return 0;
                   else
                       return term.get(n-1) / n;
               }
           }));
   }

   public String toString() {
       return toString(DISP_TERM);
   }

   public String toString(int dpTerm) {
       StringBuffer s = new StringBuffer();
       {
           double c = term.get(0);
           if (c != 0)
               s.append(c);
       }
       for (int i = 1; i < dpTerm; i++) {
           double c = term.get(i);
           if (c != 0) {
               if (c > 0 && s.length() > 0)
                   s.append("+");
               if (c == 1)
                   s.append(X_VAR);
               else if (c == -1)
                   s.append("-" + X_VAR);
               else
                   s.append(c + X_VAR);
               if (i > 1)
                   s.append(i);
           }
       }
       if (s.length() == 0)
           s.append("0");
       s.append("+...");
       return s.toString();
   }

   public static void main(String[] args) {
       FormalPS cos = new FormalPS();
       FormalPS sin = cos.intg();
       cos.copyFrom(new FormalPS(new double[]{1}).sub(sin.intg()));
       System.out.println("SIN(x) = " + sin);
       System.out.println("COS(x) = " + cos);
   }

}</lang> Output:

SIN(x) = x-0.16666666666666666x3+0.008333333333333333x5-1.984126984126984E-4x7+2.7557319223985893E-6x9-2.505210838544172E-8x11+...
COS(x) = 1.0-0.5x2+0.041666666666666664x4-0.001388888888888889x6+2.48015873015873E-5x8-2.7557319223985894E-7x10+...

Python

<lang python> \ For a discussion on pipe() and head() see

 http://paddy3118.blogspot.com/2009/05/pipe-fitting-with-python-generators.html

from itertools import izip from fractions import Fraction import sys

def head(n):

    return a generator that passes through at most n items
   
   def head(x):
       i = n
       for val in x:
           if i>0:
               i -= 1
               yield val
           else:
               raise StopIteration
   return head

def pipe(*cmds):

    pipe(a,b,c,d, ...) -> yield from ...d(c(b(a())))
   
   gen = cmds[0]
   for cmd in cmds[1:]:
       gen = cmd(gen)
   for x in gen:
       yield x

def sinepower():

   n = 0
   fac = 1
   sign = +1
   zero = 0
   yield zero
   while True:
       n +=1
       fac *= n
       yield Fraction(1, fac*sign)
       sign = -sign
       n +=1
       fac *= n
       yield zero

def cosinepower():

   n = 0
   fac = 1
   sign = +1
   yield Fraction(1,fac)
   zero = 0
   while True:
       n +=1
       fac *= n
       yield zero
       sign = -sign
       n +=1
       fac *= n
       yield Fraction(1, fac*sign)

def pluspower(*powergenerators):

   for elements in izip(*powergenerators):
       yield sum(elements)

def minuspower(*powergenerators):

   for elements in izip(*powergenerators):
       yield elements[0] - sum(elements[1:])

def mulpower(fgen,ggen):

   'From: http://en.wikipedia.org/wiki/Power_series#Multiplication_and_division'
   a,b = [],[]
   for f,g in izip(fgen, ggen):
       a.append(f)
       b.append(g)
       yield sum(f*g for f,g in izip(a, reversed(b)))

def constpower(n):

   yield n
   while True:
       yield 0

def diffpower(gen):

   'differentiatiate power series'
   n = 0
   gen.next()
   for an1 in gen:
       n +=1
       yield an1*n

def intgpower(k=0):

   'integrate power series with constant k'
   def _intgpower(gen):
       n = 0
       yield k
       an1 = gen.next()
       for an in gen:
           n +=1
           yield an1 * Fraction(1,n)
           an1 = an
   return _intgpower

print "cosine" c = [x for x in pipe(cosinepower(), head(10))] print c print "sine" s = [x for x in pipe(sinepower(), head(10))] print s

1. integrate cosine

integc = [x for x in pipe(cosinepower(),intgpower(0), head(10))]

1. 1 - (integrate sine)

integs1 = [x for x in minuspower(pipe(constpower(1), head(10)),

                                pipe(sinepower(),intgpower(0), head(10)))]

assert s == integc, "The integral of cos should be sin" assert c == integs1, "1 minus the integral of sin should be cos"</lang>

Sample output

cosine
[Fraction(1, 1), 0, Fraction(-1, 2), 0, Fraction(1, 24), 0, Fraction(-1, 720), 0, Fraction(1, 40320), 0]
sine
[0, Fraction(1, 1), 0, Fraction(-1, 6), 0, Fraction(1, 120), 0, Fraction(-1, 5040), 0, Fraction(1, 362880)]

Scheme

Definitions of operations on lazy lists: <lang scheme>(define-syntax lons

 (syntax-rules ()
   ((_ lar ldr) (delay (cons lar (delay ldr))))))

(define (lar lons)

 (car (force lons)))

(define (ldr lons)

 (force (cdr (force lons))))

(define (lap proc . llists)

 (lons (apply proc (map lar llists)) (apply lap proc (map ldr llists))))

(define (take n llist)

 (if (zero? n)
     (list)
     (cons (lar llist) (take (- n 1) (ldr llist)))))

(define (iota n)

 (lons n (iota (+ n 1))))

(define (repeat n)

 (lons n (repeat n)))</lang>

Definitions of operations on formal power series: <lang scheme>(define (fps+ . llists)

 (apply lap + llists))

(define (fps- . llists)

 (apply lap - llists))

(define (fps* . llists)

 (define (*fps* p q)
   (let ((larp (lar p)) (larq (lar q)) (ldrp (ldr p)) (ldrq (ldr q)))
     (lons (* larp larq)
           (fps+ (lap (lambda (p) (* p larp)) ldrq)
                 (lap (lambda (p) (* p larq)) ldrp)
                 (lons 0 (*fps* ldrp ldrq))))))
 (cond ((null? llists) (lons 1 (repeat 0)))
       ((null? (cdr llists)) (car llists))
       (else
        (apply fps* (cons (*fps* (car llists) (cadr llists)) (cddr llists))))))

(define (fps/ n . llists)

 (define (*fps/ n d)
   (let ((q (/ (lar n) (lar d))))
     (lons q (*fps/ (fps- (ldr n) (lap (lambda (p) (* p q)) (ldr d))) d))))
 (if (null? llists)
     (*fps/ (lons 1 (repeat 0)) n)
     (*fps/ n (apply fps* llists))))

(define (fpsint llist)

 (lons 0 (lap * llist (lap / (iota 1)))))

(define (fpsdif llist)

 (lap * (iota 1) (ldr llist)))</lang>

Now the sine and cosine functions can be defined in terms of eachother using integrals: <lang scheme>(define fpscos

 (fps- (lons 1 (repeat 0)) (fpsint (delay (force fpssin)))))

(define fpssin

 (fpsint (delay (force fpscos))))

(display (take 10 fpssin)) (newline)

(display (take 10 fpscos)) (newline)</lang> Output:

(0 1 0 -1/6 0 1/120 0 -1/5040 0 1/362880)
(1 0 -1/2 0 1/24 0 -1/720 0 1/40320 0)

Now we can do some calculations with these, e.g. show that $\sin ^{2}x+\cos ^{2}x=1$ or define $\tan x={\frac {\sin x}{\cos x}}$: <lang scheme>(display (take 10 (fps+ (fps* fpssin fpssin) (fps* fpscos fpscos)))) (newline)

(define fpstan

 (fps/ fpssin fpscos))

(display (take 10 fpstan)) (newline)</lang> Output:

(1 0 0 0 0 0 0 0 0 0)
(0 1 0 1/3 0 2/15 0 17/315 0 62/2835)

Tcl

Works with: Tcl version 8.6

Tcl doesn't arbitrary definitions of numbers without extension packages, so we'll represent these formal power series as objects, which are really just wrappers around a pair of functions: one determining how many terms there are in the series (possibly including "infinitely many") and the other producing the factor for a particular term.

This code makes extensive use of the fact that objects can have methods and variables independent of their class. This greatly reduces the requirement for singleton classes. <lang tcl>package require Tcl 8.6

oo::class create PowerSeries {

   variable name
   constructor {{body {}} args} {
       # Use the body to adapt the methods of the _object_

oo::objdefine [self] $body

       # Use the rest to configure variables in the object

foreach {var val} $args { set [my varname $var] $val }

       # Guess the name if not already set

if {![info exists [my varname name]]} { set name [namespace tail [self]] }

   }
   method name {} {

return $name

   }
   method term i {

return 0

   }
   method limit {} {

return inf

   }

   # A pretty-printer, that prints the first $terms non-zero terms
   method print {terms} {

set result "${name}(x) == " set limit [my limit] if {$limit == 0} { # Special case return $result[my term 0] } set tCount 0 for {set i 0} {$tCount<$terms && $i<=$limit} {incr i} { set t [my term $i] if {$t == 0} continue incr tCount set t [format %.4g $t]

           if {$t eq "1" && $i != 0} {set t ""}

if {$i == 0} { append result "$t + " } elseif {$i == 1} { append result "${t}x + " } else { set p [string map { 0 \u2070 1 \u00b9 2 \u00b2 3 \u00b3 4 \u2074 5 \u2075 6 \u2076 7 \u2077 8 \u2078 9 \u2079 } $i] append result "${t}x$p + " } } return [string trimright $result "+ "]

   }

   # Evaluate (a prefix of) the series at a particular x
   # The terms parameter gives the number; 5 is enough for show
   method evaluate {x {terms 5}} {

set result 0 set limit [my limit] set tCount 0 for {set i 0} {$tCount<$terms && $i<=$limit} {incr i} { set t [my term $i] if {$t == 0} continue incr tCount set result [expr {$result + $t * ($x ** $i)}] } return $result

   }

   # Operations to build new sequences from old ones
   method add {s} {

PowerSeries new { variable S1 S2 method limit {} {expr {max([$S1 limit],[$S2 limit])}} method term i { set t1 [expr {$i>[$S1 limit] ? 0 : [$S1 term $i]}] set t2 [expr {$i>[$S2 limit] ? 0 : [$S2 term $i]}] expr {$t1 + $t2} } } S1 [self] S2 $s name "$name+[$s name]"

   }
   method subtract {s} {

PowerSeries new { variable S1 S2 method limit {} {expr {max([$S1 limit],[$S2 limit])}} method term i { set t1 [expr {$i>[$S1 limit] ? 0 : [$S1 term $i]}] set t2 [expr {$i>[$S2 limit] ? 0 : [$S2 term $i]}] expr {$t1 - $t2} } } S1 [self] S2 $s name "$name-[$s name]"

   }
   method integrate Template:Name "" {

if {$Name eq ""} {set Name "Integrate\[[my name]\]"} PowerSeries new { variable S limit method limit {} { if {[info exists limit]} {return $limit} try { return [expr {[$S limit] + 1}] } on error {} { # If the limit spirals out of control, it's infinite! return [set limit inf] } } method term i { if {$i == 0} {return 0} set t [$S term [expr {$i-1}]] expr {$t / double($i)} } } S [self] name $Name

   }
   method differentiate Template:Name "" {

if {$Name eq ""} {set Name "Differentiate\[[my name]\]"} PowerSeries new { variable S method limit {} {expr {[$S limit] ? [$S limit] - 1 : 0}} method term i {expr {[incr i] * [$S term $i]}} } S [self] name $Name

   }
   # Special constructor for making constants
   self method constant n {

PowerSeries new { variable n method limit {} {return 0} method term i {return $n} } n $n name $n

   }

}

1. Define the two power series in terms of each other

PowerSeries create cos ;# temporary dummy object... rename [cos integrate "sin"] sin cos destroy  ;# remove the dummy to make way for the real one... rename [[PowerSeries constant 1] subtract [sin integrate]] cos</lang> Demonstrating: <lang tcl>% sin print 7 sin(x) == x + -0.1667x³ + 0.008333x⁵ + -0.0001984x⁷ + 2.756e-06x⁹ + -2.505e-08x¹¹ + 1.606e-10x¹³ % cos print 7 1-Integrate[sin](x) == 1 + -0.5x² + 0.04167x⁴ + -0.001389x⁶ + 2.48e-05x⁸ + -2.756e-07x¹⁰ + 2.088e-09x¹² % sin evaluate [expr acos(0)] 1.0000035425842861 % cos evaluate [expr acos(0)] 2.473727636463901e-5</lang>