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
a0 + a1 * x + a2 * x2 + a3 * x3 + ...
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
sinx = ∫ cosx cosx = 1 - ∫ sinx
Goals: Demonstrate how the language handles new numeric types and delayed (or lazy) evaluation.
D
| <d>
module fps; import q ; // a module for Rational Number import std.format : sformat = format; import std.math, std.stdio ; interface Term(U) { U coef(int n) ; } struct FPS(U) { alias Term!(U) UT ;
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 class() UT { 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.coef(0) ;
for(int i = 1 ; i <= n ; i++) {
res[i] = cast(U)0 ;
for(int j = 0 ; j < i ; j++) {
res[i] += term.coef(i - j) * res[j] ;
}
res[i] = -res[0] * res[i] ;
}
return res[n] ;
}
UF opAdd(UF rhs) { return UF(new class() UT { U coef(int n) {
return term.coef(n) + rhs.term.coef(n) ;
};});}
UF opSub(UF rhs) { return UF(new class() UT { U coef(int n) {
return term.coef(n) - rhs.term.coef(n) ;
};});}
UF opMul(UF rhs) { return UF(new class() UT { U coef(int n) {
U res = cast(U) 0 ;
for(int i = 0 ; i <=n ; i++)
res += term.coef(i)*rhs.term.coef(n -i) ; // Cauchy product with rhs
return res ;
};});}
UF opDiv(UF rhs) { return UF(new class() UT { U coef(int n) {
U res = cast(U) 0 ;
for(int i = 0 ; i <=n ; i++)
res += term.coef(i)*rhs.inverseCoef(n -i) ; // Cauchy product with inverse of rhs
return res ;
};});}
UF Diff() { return UF(new class() UT { U coef(int n) {
return term.coef(n + 1) * (n + 1) ;
};});}
UF Intg() { return UF(new class() UT { U coef(int n) {
if (n == 0) return cast(U)0 ;
return term.coef(n-1)/n ;
};});}
string toString(int dpTerm = -1) {
string s ;
if(dpTerm <= 0) dpTerm = DispTerm ;
U c = term.coef(0) ;
if(c != cast(U)0) s ~= sformat("%s", c) ;
for(int i = 1 ; i < dpTerm ; i++)
if((c = term.coef(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() { RF COS ;
RF SIN = COS.Intg ;
COS = RF([QT(1,1)]) - SIN.Intg ;
writefln("SIN(x) = ",SIN) ;
writefln("COS(x) = ",COS) ;
}</d> |
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 limited precision of long, the denominator will be quickly overflow when multiplied, but sufficient to show first few coefficients in this task.
| <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) ;
}
}</d> |
Haskell
newtype Series a = S [a] deriving (Eq, Show) 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 fs - S gs = S $ zipWith (-) fs gs S (f:ft) * S gs@(g:gt) = S $ f*g : ft*gs + map (f*) gt instance Fractional a => Fractional (Series a) where S (f:ft) / S (g:gt) = S qs where qs = f/g : map ((recip g) *) (ft-qs*gt) 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
Output (with manual interruption):
*Main> sinx S [0 % 1,1 % 1,0 % 1,(-1) % 6,0 % 1,1 % 120,0 % 1,(-1) % 5040,0 % 1,1 % 362880,Interrupted. *Main> cosx S [1 % 1,0 % 1,(-1) % 2,0 % 1,1 % 24,0 % 1,(-1) % 720,0 % 1,1 % 40320,0 % 1,(-1) % 3628800,Interrupted.