Formal power series: Difference between revisions
jq
(Simplified 'cache' argument semantics) |
(jq) |
||
Line 1,012:
=={{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|Maxima}}==
| |||