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
The ai 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
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) > Integer_Rationals.Zero then Put (" +"); Put (A (Power)); Put (" X **" & Integer'Image (Power)); elsif A (Power) < Integer_Rationals.Zero 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 ...
D
<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
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 java> 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+...