Formal power series: Difference between revisions
Content added Content deleted
(→{{header|D}}: fix & remove incorrect tag) |
m (Fixed typo in task description) |
||
| Line 7: | Line 7: | ||
The ''a<sub>i</sub>'' 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. |
The ''a<sub>i</sub>'' 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 |
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. |
'''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. |
||
Revision as of 09:55, 6 April 2008
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
module FPS:
| <d>
module fps; import std.format : sformat = format; import std.math ; interface Term(U) { U coef(int n) ; } struct FPS(U) { alias Term!(U) UT ;
alias FPS!(U) UF ;
static int EvalTerm = 8 ;
static int DispTerm = 8 ;
static string XVar = "x" ;
UT term ;
static UF opCall(UT t) { UF f ; f.term = t ; return f ;} // from a coef. generator
static UF opCall(U[] polynomial) { // from a finite 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 opCall(U x, int evTerm = -1) { // evaluate the series at x
U y = term.coef(0) ;
if( evTerm <= 0) evTerm = EvalTerm ;
for(int i = 1; i < evTerm ; i++, x *= x) {
y = y + x*term.coef(i) ;
}
return y ;
}
bool Equal(UF rhs, int upto = -1) {
if (upto <= 0) upto = EvalTerm ;
bool res = true ;
for(int i = 0 ; i < upto ; i++)
// if (term.coef(i) != rhs.term.coef(i)) /* not work due to _real_ representation is not exact, * _approxEqual_ should break this template struct's generics.
* 1e-10 below is relative difference
*/
if (!approxEqual(term.coef(i), rhs.term.coef(i), 1e-10))
res = false ;
return res ;
}
U inverseCoef(int n) { // the coef of the inverse, not all FPS has an inverse
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 opNeg() { return UF(new class() UT { U coef(int n) {
return -term.coef(n) ;
};});}
UF opPos() { return UF(term); }
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++) // Cauchy Product with rhs
res += term.coef(i)*rhs.term.coef(n -i) ;
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++) // Cauchy Product with inverse of rhs
res += term.coef(i)*rhs.inverseCoef(n -i) ;
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 ;
};});}
U IntgEval(U a, U c = cast(U)0, int evTerm = -1)
{ UF intg = Intg ; return intg(a, evTerm) - c ; }
U DiffEval(U a, int evTerm = -1)
{ UF diff = Diff ; return diff(a, evTerm) ; }
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 ;
}
} class ternaryDIVEval(U) { alias FPS!(U) UF ;
alias ternaryDIVEval!(U) DV ;
U evalValue, tempValue ;
this(U value) { evalValue = value ; }
DV opDiv_r(UF lhs) { tempValue = lhs(evalValue) ; return this ; }
U opDiv(UF rhs) { return tempValue / rhs(evalValue) ; }
} ternaryDIVEval!(U) DiV(U)(U value) { return new ternaryDIVEval!(U)(value) ;} </d> |
Test program:
| <d>
module fpx; import fps ; import std.stdio, std.math ; alias FPS!(real) RF ; alias Term!(real) RT ; void main() { RF COS ; RF SIN = COS.Intg ; COS = RF([1]) - SIN.Intg ; RF DCOS = COS.Diff ; RF DSIN = SIN.Diff ; RF ICOS = COS.Intg ; RF ISIN = SIN.Intg ; writefln("Is COS == 1 - \∫SIN ? %s", COS.Equal(RF([1]) - ISIN, 128)) ;
writefln() ;
writefln("SIN'(x) = ",DSIN) ;
writefln("COS(x) = ",COS) ;
writefln("1-\∫SIN(x) = ",RF([1]) - ISIN) ;
writefln() ;
writefln("-COS'(x) = ", -DCOS) ;
writefln("SIN(x) = ",SIN) ;
writefln("\∫COS(x) = ",+ICOS) ;
writefln() ;
real NR = 128.0 ;
real RY = 1/NR*PI ; // smaller value for faster convergence
writefln("y = %.8f", RY) ;
writefln() ;
writefln("sin(y) = %.19f <- std.math function",sin(RY)) ;
writefln("SIN(y) = %.19f",SIN(RY)) ;
writefln("\∫COS(y) = %.19f",ICOS(RY)) ;
writefln("SIN'(y) = %.19f",DSIN(RY)) ;
writefln("COS(y) = %.19f",COS(RY)) ;
writefln("cos(y) = %.19f <- std.math function",cos(RY)) ;
writefln() ;
RF SxC = SIN * COS ;
writefln("(SIN*COS)(x) = ",SxC) ;
writefln("sin(2y) = %.19f <- std.math function",sin(2*RY)) ;
writefln("SIN(2y) = %.19f",SIN(2*RY)) ;
writefln("2*SIN(y)*COS(y) = %.19f",2*SIN(RY)*COS(RY)) ;
writefln("2*(SIN*COS)(y) = %.19f",2*SxC(RY)) ;
writefln() ;
RF SdC = SIN / COS ;
writefln("(SIN/COS)(x) = ",SdC) ;
writefln("sin(y)/cos(y) = %.19f <- std.math function",sin(RY)/cos(RY)) ;
writefln("SIN/DiV(y)/COS = %.19f",SIN/DiV(RY)/COS) ;
writefln("SIN(y)/COS(y) = %.19f",SIN(RY)/COS(RY)) ;
writefln("(SIN/COS)(y) = %.19f",SdC(RY)) ;
writefln() ;
RF UNI = COS*COS + SIN*SIN ;
writefln("compare COS(y)*COS(y)+SIN(y)*SIN(y) with (COS*COS+SIN*SIN)(y)") ;
for(int i=-4 ; i <= 4 ; i += 2) {
real Y = i/NR*PI ;
real a = COS(Y)*COS(Y) + SIN(Y)*SIN(Y) ;
real b = UNI(Y) ;
writefln("y = %+.8f : %.19f - %.8f",Y, a, b) ;
}
}</d> |
Output:
Is COS == 1 - ∫SIN ? true SIN'(x) = 1-0.5x2+0.0416666x4-0.00138888x6+... COS(x) = 1-0.5x2+0.0416666x4-0.00138888x6+... 1-∫SIN(x) = 1-0.5x2+0.0416666x4-0.00138888x6+... -COS'(x) = x-0.166667x3+0.00833333x5-0.000198412x7+... SIN(x) = x-0.166667x3+0.00833333x5-0.000198412x7+... ∫COS(x) = x-0.166667x3+0.00833333x5-0.000198412x7+... y = 0.02454369 sin(y) = 0.0245412285229122880 <- std.math function SIN(y) = 0.0245436321266466245 ∫COS(y) = 0.0245436321266466245 SIN'(y) = 0.9996988035766323983 COS(y) = 0.9996988035766323983 cos(y) = 0.9996988186962042201 <- std.math function (SIN*COS)(x) = x-0.666667x3+0.133333x5-0.0126984x7+... sin(2y) = 0.0490676743274180142 <- std.math function SIN(2y) = 0.0490864175399623567 2*SIN(y)*COS(y) = 0.0490724793448672567 2*(SIN*COS)(y) = 0.0490869013761514380 (SIN/COS)(x) = x+0.333333x3+0.133333x5+0.0539682x7+... sin(y)/cos(y) = 0.0245486221089254441 <- std.math function SIN/DiV(y)/COS = 0.0245510268081112297 SIN(y)/COS(y) = 0.0245510268081112297 (SIN/COS)(y) = 0.0245438135652175300 compare COS(y)*COS(y)+SIN(y)*SIN(y) with (COS*COS+SIN*SIN)(y) y = -0.09817477 : 1.0000262651249077323 - 1.00000000 y = -0.04908738 : 1.0000015465133229911 - 1.00000000 y = +0.00000000 : 1.0000000000000000000 - 1.00000000 y = +0.04908738 : 1.0000013565112958463 - 1.00000000 y = +0.09817477 : 1.0000201850600390979 - 1.00000000
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.