Formal power series: Difference between revisions
Rename Perl 6 -> Raku, alphabetize, minor clean-up
SqrtNegInf (talk | contribs) m (→{{header|Perl}}: typo) |
Thundergnat (talk | contribs) (Rename Perl 6 -> Raku, alphabetize, minor clean-up) |
||
Line 207:
Sin = + 1 X ** 1 - 1 / 6 X ** 3 + 1 / 120 X ** 5 - 1 / 5040 X ** 7 + 1 / 362880X ** 9 - 1 / 39916800 X ** 11 ...
</pre>
=={{header|C}}==
Following is a simple implementation of formal power series in C. It's not "new datatype for the language" per se, but does demonstrate how lazy evaluation and infinite list generation can be done for this task. Note that, to be of real use, one should also cache terms looked up and free up memory. Both are trivially done (I actually had them, but removed them for simplicity).
<lang C>#include <stdio.h>
#include <stdlib.h>
#include <math.h> /* for NaN */
enum fps_type {
FPS_CONST = 0,
FPS_ADD,
FPS_SUB,
FPS_MUL,
FPS_DIV,
FPS_DERIV,
FPS_INT,
};
typedef struct fps_t *fps;
typedef struct fps_t {
int type;
fps s1, s2;
double a0;
} fps_t;
fps fps_new()
{
fps x = malloc(sizeof(fps_t));
x->a0 = 0;
x->s1 = x->s2 = 0;
x->type = 0;
return x;
}
/* language limit of C; when self or mutual recursive definition is needed,
* one has to be defined, then defined again after it's used. See how
* sin and cos are defined this way below
*/
void fps_redefine(fps x, int op, fps y, fps z)
{
x->type = op;
x->s1 = y;
x->s2 = z;
}
fps _binary(fps x, fps y, int op)
{
fps s = fps_new();
s->s1 = x;
s->s2 = y;
s->type = op;
return s;
}
fps _unary(fps x, int op)
{
fps s = fps_new();
s->s1 = x;
s->type = op;
return s;
}
/* Taking the n-th term of series. This is where actual work is done. */
double term(fps x, int n)
{
double ret = 0;
int i;
switch (x->type) {
case FPS_CONST: return n > 0 ? 0 : x->a0;
case FPS_ADD:
ret = term(x->s1, n) + term(x->s2, n); break;
case FPS_SUB:
ret = term(x->s1, n) - term(x->s2, n); break;
case FPS_MUL:
for (i = 0; i <= n; i++)
ret += term(x->s1, i) * term(x->s2, n - i);
break;
case FPS_DIV:
if (! term(x->s2, 0)) return NAN;
ret = term(x->s1, n);
for (i = 1; i <= n; i++)
ret -= term(x->s2, i) * term(x, n - i) / term(x->s2, 0);
break;
case FPS_DERIV:
ret = n * term(x->s1, n + 1);
break;
case FPS_INT:
if (!n) return x->a0;
ret = term(x->s1, n - 1) / n;
break;
default:
fprintf(stderr, "Unknown operator %d\n", x->type);
exit(1);
}
return ret;
}
#define _add(x, y) _binary(x, y, FPS_ADD)
#define _sub(x, y) _binary(x, y, FPS_SUB)
#define _mul(x, y) _binary(x, y, FPS_MUL)
#define _div(x, y) _binary(x, y, FPS_DIV)
#define _integ(x) _unary(x, FPS_INT)
#define _deriv(x) _unary(x, FPS_DERIV)
fps fps_const(double a0)
{
fps x = fps_new();
x->type = FPS_CONST;
x->a0 = a0;
return x;
}
int main()
{
int i;
fps one = fps_const(1);
fps fcos = fps_new(); /* cosine */
fps fsin = _integ(fcos); /* sine */
fps ftan = _div(fsin, fcos); /* tangent */
/* redefine cos to complete the mutual recursion; maybe it looks
* better if I said
* *fcos = *( _sub(one, _integ(fsin)) );
*/
fps_redefine(fcos, FPS_SUB, one, _integ(fsin));
fps fexp = fps_const(1); /* exponential */
/* make exp recurse on self */
fps_redefine(fexp, FPS_INT, fexp, 0);
printf("Sin:"); for (i = 0; i < 10; i++) printf(" %g", term(fsin, i));
printf("\nCos:"); for (i = 0; i < 10; i++) printf(" %g", term(fcos, i));
printf("\nTan:"); for (i = 0; i < 10; i++) printf(" %g", term(ftan, i));
printf("\nExp:"); for (i = 0; i < 10; i++) printf(" %g", term(fexp, i));
return 0;
}</lang>Output:<pre>Sin: 0 1 0 -0.166667 0 0.00833333 0 -0.000198413 0 2.75573e-06
Cos: 1 0 -0.5 0 0.0416667 0 -0.00138889 0 2.48016e-05 0
Tan: 0 1 0 0.333333 0 0.133333 0 0.0539683 0 0.0218695
Exp: 1 1 0.5 0.166667 0.0416667 0.00833333 0.00138889 0.000198413 2.48016e-05 2.75573e-06</pre>
=={{header|Clojure}}==
Line 440 ⟶ 588:
(* 3 (force *sinx*))))))
(3 -3 -1/2 1/2 1/24 -1/40 -1/720 1/1680 1/40320 -1/120960)</pre>
=={{header|D}}==
Line 670:
</lang>
=={{header|Elisa}}==
The generic component FormalPowerSeries is instantiated with the Rational type as provided by the task [[Arithmetic/Rational]]
<lang Elisa>component FormalPowerSeries(Number);
type PowerSeries;
PowerSeries(Size = integer) -> PowerSeries;
+ PowerSeries -> PowerSeries;
- PowerSeries -> PowerSeries;
PowerSeries + PowerSeries -> PowerSeries;
PowerSeries - PowerSeries -> PowerSeries;
PowerSeries * PowerSeries -> PowerSeries;
Integral(PowerSeries) -> PowerSeries;
Differential(PowerSeries) -> PowerSeries;
Zero -> PowerSeries;
One -> PowerSeries;
Array(PowerSeries) -> array(Number);
begin
PowerSeries(Size) = PowerSeries:[T = array(Number, Size); Size];
+ A = A;
- A = [ C = PowerSeries(A.Size);
[ i = 1 .. A.Size; C.T[i] := - A.T[i] ];
C];
A + B = [ if A.Size > B.Size then return(B + A);
C = PowerSeries(B.Size);
[ i = 1 .. A.Size; C.T[i] := A.T[i] + B.T[i] ];
[ i = (A.Size +1) .. B.Size; C.T[i] := B.T[i] ];
C];
A - B = A + (- B );
A * B = [ C = PowerSeries(A.Size + B.Size - 1);
[ i = 1 .. A.Size;
[j = 1.. B.Size;
C.T[i + j - 1] := C.T[i + j - 1] + A.T[i] * B.T[j] ] ];
C];
Integral(A) = [ if A.Size == 0 then return (A);
C = PowerSeries(A.Size + 1);
[ i = 1 .. A.Size; C.T[i +1] := A.T[i] / Number( i )];
C.T[1]:= Number(0);
C ];
Differential(A) = [ if A.Size == 1 then return (A);
C = PowerSeries(A.Size - 1);
[ i = 1 .. C.Size; C.T[i] := A.T[i + 1] * Number( i )];
C ];
Zero = [ C = PowerSeries (1); C.T[1]:= Number(0); C];
One = [ C = PowerSeries (1); C.T[1]:= Number(1); C];
Array(PowerSeries) -> array(Number);
Array(TS) = TS.T;
end component FormalPowerSeries;
</lang>
Tests
<lang Elisa>use RationalNumbers;
use FormalPowerSeries(Rational);
X => symbol;
term + term => term;
term / term => term;
term * term => term;
symbol ** integer => term;
Output(text,PowerSeries) -> term;
Output(Name,PS) = [ E1 := term:symbol(Name); E2:= null(term);
[ i = 1..size(Array(PS));
Num = Numerator(Array(PS)[i]);
if Num <> 0 then
[ E2:= term: Num / term: Denominator(Array(PS)[i]) * X ** (i-1);
E1:= E1 + E2 ];
];
E1];
Cos(integer) -> PowerSeries;
Cos(Limit) = [ if Limit == 1 then return(One);
( One - Integral(Integral(Cos (Limit - 1)))) ];
Sin(integer) -> PowerSeries;
Sin(Limit) = Integral(Cos (Limit));
Output("cos = ",Cos(5))?
Output("sin = ",Sin(5))?
</lang>
Output
<lang Elisa>cos = + 1 / 1 * X ** 0 + -1 / 2 * X ** 2 + 1 / 24 * X ** 4 + -1 / 720 * X ** 6 + 1 / 40320 * X ** 8
sin = + 1 / 1 * X ** 1 + -1 / 6 * X ** 3 + 1 / 120 * X ** 5 + -1 / 5040 * X ** 7 + 1 / 362880 * X ** 9</lang>
=={{header|Go}}==
Line 909 ⟶ 1,005:
cos: 1.00000 0.00000 -0.50000 0.00000 0.04167 0.00000
</pre>
=={{header|Haskell}}==
Line 1,091:
=={{header|Java}}==
See [[Formal power series/Java]]
=={{header|jq}}==
{{works with|jq|1.4}}
====Introduction and Examples====
Since a formal power series can be viewed as a function from the non-negative integers onto a suitable range,
we shall identify a jq filter that maps integers to the appropriate range as a power series. For example, the jq function
<lang jq>1/(1+.)</lang>
represents the power series 1 + x/2 + x/3 + ... because 1/(1+.) maps i to 1/(i+1).
Similarly, the jq filter 1 (i.e. the filter that always returns 1) represents the power series Σ x^i.
The exponential power series, Σ (x^i)/i!, can be represented in jq by the filter:
<lang jq>1/factorial</lang>
assuming "factorial" is defined in the usual way:
<lang jq>def factorial:
reduce range(1; . + 1) as $i
(1; . * $i);</lang>
For ease of reference, we shall also define ps_exp as 1/factorial:
<lang jq>def ps_exp: 1/factorial;</lang>
In a later subsection of this article, we will define another function, ps_evaluate(p), for evaluating the power series, p, at the value specified by the input, so for example:
<lang jq>1 | ps_evaluate(ps_exp)</lang>
should evaluate to the number e approximately; using the version of ps_evaluate defined below, we find:
<lang jq>1 | ps_evaluate(1/factorial)</lang>
evaluates to 2.7182818284590455.
The following function definitions are useful for other power series:
<lang jq>def pow(n):
. as $x | n as $n
| reduce range(0;$n) as $i (1; . * $x);</lang>
For example, the power series 1 + Σ ( (x/i)^n ) where the summation is over i>0 can be written:
<lang jq>1/pow(.)</lang>
The power series representation of ln(1 + x) is as follows:
<lang jq># ln(1+x) = x - x^2 / 2 + ...
def ln_1px:
def c: if . % 2 == 0 then -1 else 1 end;
. as $i | if $i == 0 then 0 else ($i|c) / $i end;
</lang>
jq numbers are currently implemented using IEEE 754 64-bit arithmetic, and therefore this article will focus on power series that can be adequately represented using jq numbers. However, the approach used here can be used for power series defined on other domains, e.g. rationals, complex numbers, and so on.
====Finite power series====
To make it easy to represent finite power series, we define poly(ary) as follows:
<lang jq>def poly(ary): ary[.] // 0;</lang>
For example, poly( [1,2,3] ) represents the finite power series: 1 + 2x + 3x^2.
(The "// 0" ensures that the result is 0 for integers that are out-of-range with respect to the array.)
====Addition and Subtraction====
jq's "+" operator can be used to add two power series with the intended semantics;
for example:
<lang jq>(poly([1,2,3]) + poly([-1,-2,-3]))</lang>
is equal to poly([]), i.e. 0.
This is simply because in jq, (i | (poly([1,2,3]) + poly([-1,-2,-3]))) evaluates to (i | (poly([1,2,3])) + (i|poly([-1,-2,-3]))).
Subtraction works in the same way and for the same reason. The product of two power series, however, must be handled specially.
====Multiplication, Differentiation and Integration====
<lang jq># Multiply two power series, s and t:
def M(s;t):
. as $i | reduce range(0; 1+$i) as $k
(0; . + ($k|s) * (($i - $k)|t));
# Derivative of the power series, s:
def D(s): (. + 1) as $i | $i * ($i|s);
# Integral of the power series, s,
# with an integration constant equal to 0:
def I(s):
. as $i
| if $i == 0 then 0 else (($i-1)|s) /$i end;
</lang>
====Equality and Evaluation====
The following function, ps_equal(s;t;k;eps) will check whether the first k coefficients of the two power series agree to within eps:
<lang jq>def ps_equal(s; t; k; eps):
def abs: if . < 0 then -. else . end;
reduce range(0;k) as $i
(true;
if . then ((($i|s) - ($i|t))|abs) <= eps
else .
end);</lang>
To evaluate a power series, P(x), at a particular point, say y, we
can define a function, ps_evaluate(p), so that (y|ps_evaluate(p))
evaluates to P(y), assuming that P(x) converges sufficiently rapidly
to a value that can be represented using IEEE 754 64-bit arithmetic.
<lang jq># evaluate p(x) based on the first k terms of polynomial p, where x is the input
def ps_eval(p; k):
. as $x
| reduce range(0;k) as $i
# state: [sum, x^i]
([0, 1];
.[1] as $xn
| ($i|p) as $coeff
| [ .[0] + $coeff * $xn, $x * $xn])
| .[0];
# If |x| < 1 then ps_evaluate(x) will evaluate to p(x) with high precision
# if the coefficients of the polynomial are eventually bounded.
#
# WARNING: ps_evaluate(p) will not detect divergence and is not intended to
# produce accurate results unless the terms of p(x) are reasonably well-behaved.
# For |x| > 1, the result will be null if x^n overflows before convergence is achieved.
#
def ps_evaluate(p):
def abs: if . < 0 then -. else . end;
def eval(p;x):
# state: [i, x^i, sum of i terms, delta, prevdelta]
recurse(
.[0] as $i
| .[1] as $xi
| .[2] as $sum
| .[3] as $delta
| .[4] as $prevdelta
| if $delta < 1e-17 and $prevdelta < 1e-17
and ( $xi < 1e-100
or ( $sum != 0 and
(($delta/$sum) | abs) < 1e-10 and
(($prevdelta/$sum) | abs) < 1e-10) )
then empty
else
($xi * ($i|p)) as $newdelta
| [ $i + 1,
x*$xi,
$sum+$newdelta,
($newdelta|abs), $delta]
end ) ;
. as $x
| [0, 1, 0, 1, 1]
| reduce eval(p; $x) as $vector (0; $vector[2]);</lang>
====Examples====
<lang jq># Utility functions:
def abs: if . < 0 then -. else . end;
# The power series whose only non-zero coefficient is 1 at x^i:
def ps_at(i): if . == i then 1 else 0 end;
# Create an array consisting of the first . coefficients of the power series, p:
def ps_to_array(p): . as $in | reduce range(0;$in) as $i ([]; . + [$i|p]);
def pi: 4 * (1|atan);
</lang>
''''cos == I(sin)''''
<lang jq>
# Verify that the first 100 terms of I(cos) and of sin are the same:
ps_equal( I(ps_cos); ps_sin; 100; 1e-15)
# => true
# Verify that the two power series agree when evaluated at pi:
((pi | ps_evaluate(I(ps_cos))) - (pi | ps_evaluate(ps_sin))) | abs < 1e-15
# => true</lang>
''''cos == 1 - I(sin)''''
<lang jq># Verify that the first 100 terms of cos and (1 - I(sin)) are the same:
ps_equal( ps_cos; ps_at(0) - I(ps_sin); 100; 1e-5)
# => true
# Verify that the two power series agree at pi:
((pi | ps_evaluate(ps_cos)) - (pi | ps_evaluate(ps_at(0) - I(ps_sin)))) | abs < 1e-15
# => true</lang>
=={{header|Julia}}==
Line 1,516 ⟶ 1,684:
COS(x) = 1 - 1/2x^2 + 1/24x^4 - 1/720x^6 + 1/40320x^8 - 1/3628800x^10 + ...
</pre>
=={{header|Lua}}==
Line 1,768 ⟶ 1,752:
{{out}}
<pre>O[x]^9</pre>
=={{header|Maxima}}==
<lang maxima>deftaylor(f(x), sum(n! * x^n, n, 0, inf))$
taylor(f(x), x, 0, 10);
/ * 1 + x + 2 * x^2 + 6 * x^3 + 24 * x^4 + 120 * x^5 + 720 * x^6 + 5040 * x^7 + 40320 * x^8 + 362880 * x^9 + 3628800 * x^10 + ... * /
taylor(f(x)^2, x, 0, 10);
/ * 1 + 2 * x + 5 * x^2 + 16 * x^3 + 64 * x^4 + 312 * x^5 + 1812 * x^6 + 12288 * x^7 + 95616 * x^8 + 840960 * x^9 + 8254080 * x^10 + ... * /
deftaylor(fcos(x), sum((-1)^n * x^(2 * n) / (2 * n)!, n, 0, inf))$
deftaylor(fsin(x), sum((-1)^n * x^(2 * n + 1) / (2 * n + 1)!, n, 0, inf))$
taylor(fcos(x)^2 + fsin(x)^2, x, 0, 20);
/ * 1 + ... * /</lang>
=={{header|PARI/GP}}==
Line 1,909:
For a version which *does* use proper lazy lists, see [[Formal power series/Perl]]
=={{header|Phix}}==
Line 2,553 ⟶ 2,462:
'(0 1 0 1/3 0 2/15 0 17/315 0 62/2835)
</lang>
=={{header|Raku}}==
(formerly Perl 6)
<lang perl6>class DerFPS { ... }
class IntFPS { ... }
role FPS {
method coeffs { ... }
method differentiate { DerFPS.new(:x(self)) }
method integrate { IntFPS.new(:x(self)) }
method pretty($n) {
sub super($i) { $i.trans('0123456789' => '⁰¹²³⁴⁵⁶⁷⁸⁹') }
my $str = $.coeffs[0];
for flat 1..$n Z $.coeffs[1..$n] -> $p, $c {
when $c > 0 { $str ~= " + { $c .nude.join: '/'}∙x{super($p)}" }
when $c < 0 { $str ~= " - {-$c .nude.join: '/'}∙x{super($p)}" }
}
$str;
}
}
class ExplicitFPS does FPS { has @.coeffs }
class SumFPS does FPS {
has FPS ($.x, $.y);
method coeffs { $.x.coeffs Z+ $.y.coeffs }
}
class DifFPS does FPS {
has FPS ($.x, $.y);
method coeffs { $.x.coeffs Z- $.y.coeffs }
}
class ProFPS does FPS {
has FPS ($.x, $.y);
method coeffs { (0..*).map: { [+] ($.x.coeffs[0..$_] Z* $.y.coeffs[$_...0]) } }
}
class InvFPS does FPS {
has FPS $.x;
method coeffs {
# see http://en.wikipedia.org/wiki/Formal_power_series#Inverting_series
flat gather {
my @a = $.x.coeffs;
@a[0] != 0 or fail "Cannot invert power series with zero constant term.";
take my @b = (1 / @a[0]);
take @b[$_] = -@b[0] * [+] (@a[1..$_] Z* @b[$_-1...0]) for 1..*;
}
}
}
class DerFPS does FPS {
has FPS $.x;
method coeffs { (1..*).map: { $_ * $.x.coeffs[$_] } }
}
class IntFPS does FPS {
has FPS $.x;
method coeffs { 0, |(0..*).map: { $.x.coeffs[$_] / ($_+1) } }
}
class DeferredFPS does FPS {
has FPS $.realized is rw;
method coeffs { $.realized.coeffs }
}
# some arithmetic operations for formal power series
multi infix:<+>(FPS $x, FPS $y) { SumFPS.new(:$x, :$y) }
multi infix:<->(FPS $x, FPS $y) { DifFPS.new(:$x, :$y) }
multi infix:<*>(FPS $x, FPS $y) { ProFPS.new(:$x, :$y) }
multi infix:</>(FPS $x, FPS $y) { $x * InvFPS.new(:x($y)) }
# an example of a mixed-type operator:
multi infix:<->(Numeric $x, FPS $y) { ExplicitFPS.new(:coeffs(lazy flat $x, 0 xx *)) - $y }
# define sine and cosine in terms of each other
my $sin = DeferredFPS.new;
my $cos = 1 - $sin.integrate;
$sin.realized = $cos.integrate;
# define tangent in terms of sine and cosine
my $tan = $sin / $cos;
say 'sin(x) ≈ ' ~ $sin.pretty(10);
say 'cos(x) ≈ ' ~ $cos.pretty(10);
say 'tan(x) ≈ ' ~ $tan.pretty(10);</lang>
{{out}}
<pre>sin(x) ≈ 0 + 1/1∙x¹ - 1/6∙x³ + 1/120∙x⁵ - 1/5040∙x⁷ + 1/362880∙x⁹
cos(x) ≈ 1 - 1/2∙x² + 1/24∙x⁴ - 1/720∙x⁶ + 1/40320∙x⁸ - 1/3628800∙x¹⁰
tan(x) ≈ 0 + 1/1∙x¹ + 1/3∙x³ + 2/15∙x⁵ + 17/315∙x⁷ + 62/2835∙x⁹</pre>
=={{header|Scheme}}==
| |||