Formal power series: Difference between revisions
Content deleted Content added
m →{{header|Perl 6}}: Clean up some unneeded semicolons. |
Added FreeBASIC |
||
(39 intermediate revisions by 20 users not shown) | |||
Line 20:
=={{header|Ada}}==
The Taylor series package is generic to be instantiated with any rational type implementation provided by the task [[Rational Arithmetic]]:
<
generic
Line 42:
Zero : constant Taylor_Series := (0 => Rational_Numbers.Zero);
One : constant Taylor_Series := (0 => Rational_Numbers.One);
end Generic_Taylor_Series;</
The package implementation:
<
function Normalize (A : Taylor_Series) return Taylor_Series is
begin
Line 146:
end Value;
end Generic_Taylor_Series;</
The procedure Normalize is used to truncate the series when the coefficients are zero. The summation of a series (function Value) uses [[wp:Horner_scheme|Horner scheme]].
===Test task===
<
with Generic_Taylor_Series;
Line 201:
Put ("Cos ="); Put (Cos (5)); Put_Line (" ...");
Put ("Sin ="); Put (Integral (Cos (5))); Put_Line (" ...");
end Test_Taylor_Series;</
Sample output:
<pre>
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).
<syntaxhighlight 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;
}</syntaxhighlight>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 213 ⟶ 361:
First addition and subtraction. Note that most of the complication arises in allowing for finite and infinite FPSs; if only infinite power series were at issue, the function ''(defn ips+ [ips0 ips1] (map + ips0 ips1))'' would suffice.
<
(letfn [(+zs [ps] (concat ps (repeat :z)))
(notz? [a] (not= :z a))
Line 220 ⟶ 368:
(take-while notz? (map z+ (+zs ps0) (+zs ps1)))))
(defn ps- [ps0 ps1] (ps+ ps0 (map - ps1)))</
Multiplication next; again most of the complication is dealing with both finite and infinite FPS. This function explicitly uses the standard function ''lazy-seq'' to define the product sequence.
<
([ps0 ps1] (ps* [0] ps0 ps1))
([[a0 & resta] [p0 & rest0] [p1 & rest1 :as ps1]]
Line 230 ⟶ 378:
(let [mrest1 (if (or (nil? rest1) (zero? p0)) nil, (map #(* p0 %) rest1))
accum (cond (nil? resta) mrest1, (nil? mrest1) resta, :else (ps+ resta mrest1))]
(if (nil? rest0) accum, (ps* (or accum [0]) rest0 ps1)))))))</
As with most of the other examples on this page, there's no definition for division. Mathematically, FPS is a commutative ring (in fact a Euclidean domain), but not a field: the set of FPSs is not closed under division.
Now we can define integration and differentiation:
<
(defn differentiate [ps]
Line 240 ⟶ 388:
(defn integrate [ps]
(cons 0 (for [[n a] (indexed ps)] (/ a (inc n)))))</
Some examples of using these functions; in each case a ''println'' call (which forces the lazy sequence) is followed by a comment showing the output.
<
; (4 6 5)
(println (ps* [1 2] [3 4 5]))
; (3 10 13 10)</
And some examples using infinite FPSs. First define the sequence of factorials (''facts''), then define ''sin'' and ''cos'' Taylor series.
<
(def facts (map first nfacts))
Line 261 ⟶ 409:
(println (take 20 (ps+ (ps* sin sin) (ps* cos cos))))
; (1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0)</
By using ''letfn'', which supports defining mutually recursive functions, we can define the ''sin'' and ''cos'' power series directly in terms of integrals of the other series:
<
(fcos [] (ps- [1] (integrate (fsin))))]
(def sinx (fsin))
Line 269 ⟶ 417:
(println (take 10 sinx))
; (0 1 0 -1/6 0 1/120 0 -1/5040 0 1/362880)</
=={{header|Common Lisp}}==
Line 275 ⟶ 423:
Common Lisp isn't lazy, and doesn't define the arithmetic operators as generic functions. As such, this implementation defines lazy primitives (delay and force), and a lazy list built on top of them (lons, lar, ldr). This implementation also defines a package #:formal-power-series which uses all the symbols of the "COMMON-LISP" package except for +, -, *, and /, which are shadowed. Shadowing these symbols allows for definitions of generic functions which can be specialized for the power series, can default to the normal CL arithmetic operations for other kinds of objects, and can do "the right thing" for mixes thereof.
<
(:nicknames #:fps)
(:use "COMMON-LISP")
Line 281 ⟶ 429:
#:+ #:- #:* #:/))
(in-package #:formal-power-series)</
===Lazy Primitives===
<
thunk value)
Line 299 ⟶ 447:
(t (let ((val (funcall (promise-thunk object))))
(setf (promise-thunk object) nil
(promise-value object) val)))))</
===Lazy Lists===
<
lar
ldr)
Line 317 ⟶ 465:
(defmacro lons (lar ldr)
`(make-lons :lar ,lar :ldr (delay ,ldr)))</
A few utilities to make working with lazy lists easier.
<
(funcall function (lar lons) (ldr lons)))
Line 347 ⟶ 495:
(defun up-from (n)
(lons n (up-from (1+ n))))</
===Formal Power Series===
Line 353 ⟶ 501:
The mathematical operations here are translations of the Haskell code, but we specialize the operations in various ways so that behavior for normal numeric operations is preserved.
<
coeffs)
Line 417 ⟶ 565:
(defun diff (f)
(series (maplar '* (ldr (coeffs f)) (up-from 1))))</
===Example===
<
(locally (declare (special *cosx*))
(delay (int (force *cosx*)))))
(defparameter *cosx*
(delay (- 1 (int *sinx*))))</
<pre>FPS > (force-list (take 10 (coeffs (force *sinx*))))
Line 441 ⟶ 589:
(3 -3 -1/2 1/2 1/24 -1/40 -1/720 1/1680 1/40320 -1/120960)</pre>
=={{header|
See [[Formal_power_series/D]].
=={{header|EchoLisp}}==
We implement infinite formal power series (FPS) using '''streams'''. No operator overloading in EchoLisp, so we provide the operators '''s-add, s-mul''' ,.. which implement the needed operations. '''poly->stream''' converts a finite polynomial into an infinite FPS, and '''s-value''' gives the value of a FPS at x.
<syntaxhighlight lang="scheme">
(require 'math)
;; converts a finite polynomial (a_0 a_1 .. a_n) to an infinite serie (a_0 ..a_n 0 0 0 ...)
(define (poly->stream list)
(make-stream (lambda(n) (cons (if (< n (length list)) (list-ref list n) 0) (1+ n))) 0))
;; c = a + b , c_n = a_n + b_n
(define (s-add a b)
(make-stream (lambda (n) (cons (+ (stream-ref a n) (stream-ref b n)) (1+ n))) 0))
;; c = a * b , c_n = ∑ (0 ..n) a_i * b_n-i
(define (s-mul-coeff n a b) (sigma (lambda(i) (* (stream-ref a i)(stream-ref b (- n i)))) 0 n))
(define (s-mul a b)
(make-stream (lambda(n) (cons (s-mul-coeff n a b) (1+ n))) 0))
;; b = 1/a ; b_0 = 1/a_0, b_n = - ∑ (1..n) a_i * b_n-i / a_0
(define (s-inv-coeff n a b)
(if (zero? n) (/ (stream-ref a 0))
(- (/ (sigma (lambda(i) (* (stream-ref a i)(stream-ref b (- n i)))) 1 n)
(stream-ref a 0)))))
;; note the self keyword which refers to b = (s-inv a)
(define (s-inv a)
(make-stream (lambda(n) (cons (s-inv-coeff n a self ) (1+ n))) 0))
;; b = (s-k-add k a) = k + a_0, a_1, a_2, ...
(define (s-k-add k a)
(make-stream (lambda(n) (cons
(if(zero? n) (+ k (stream-ref a 0)) (stream-ref a n)) (1+ n))) 0))
;; b = (s-neg a) = -a_0,-a_1, ....
(define (s-neg a)
(make-stream (lambda(n) (cons (- (stream-ref a n)) (1+ n))) 0))
;; b = (s-int a) = ∫ a ; b_0 = 0 by convention, b_n = a_n-1/n
(make-stream (lambda(n) (cons (if (zero? n) 0 (/ (stream-ref a (1- n)) n)) (1+ n))) 0))
;; value of power serie at x, n
(define (s-value a x (n 20))
(poly x (take a n)))
;; stream-cons allows mutual delayed references
;; sin = ∫ cos
(define sin-x (stream-cons 0 (stream-rest (s-int cos-x))))
;; cos = 1 - ∫ sin
(define cos-x (stream-cons 1 (stream-rest (s-k-add 1 (s-neg (s-int sin-x))))))
</syntaxhighlight>
{{out}}
<syntaxhighlight lang="scheme">
(take cos-x 16)
→ (1 0 -1/2 0 1/24 0 -1/720 0 1/40320 0 -1/3628800 0 1/479001600 0 -1.1470745597729725e-11 0)
(take sin-x 16)
→ (0 1 0 -1/6 0 1/120 0 -1/5040 0 1/362880 0 -1/39916800 0 1.6059043836821613e-10 0 -7.647163731819816e-13)
;; compute (cos PI)
(s-value cos-x PI)
→ -1.0000000035290808
;; check that 1 / (1 - x) = 1 + x + x^1 +
(define fps-1 (poly->stream '( 1 -1)))
(take fps-1 13)
→ (1 -1 0 0 0 0 0 0 0 0 0 0
(define inv-fps-1 (s-inv fps-1))
(take inv-fps-1 13)
→ (1 1 1 1 1 1 1 1 1 1 1
(s-value inv-fps-1 0.5) ;; check that 1 / (1 - 0.5) = 2
→ 1.9999980926513672
(s-value inv-fps-1 0.5 100) ;; 100 terms
→ 2
</syntaxhighlight>
=={{header|Elisa}}==
The generic component FormalPowerSeries is instantiated with the Rational type as provided by the task [[Arithmetic/Rational]]
<syntaxhighlight lang="elisa">component FormalPowerSeries(Number);
type PowerSeries;
PowerSeries(Size = integer) -> 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) =
+
-
[
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
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;
</syntaxhighlight>
Tests
<syntaxhighlight 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))?
</syntaxhighlight>
Output
<syntaxhighlight 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</syntaxhighlight>
=={{header|FreeBASIC}}==
{{trans|C}}
<syntaxhighlight lang="vbnet">Type fps_t
As Integer tipo
As fps_t Ptr s1, s2
As Single a0
End Type
Enum fps_type
FPS_CONST = 0
FPS_ADD
FPS_SUB
FPS_MUL
FPS_DIV
FPS_DERIV
FPS_INT
End Enum
Function fps_new() As fps_t Ptr
Dim As fps_t Ptr x = Callocate(1, Sizeof(fps_t))
x->a0 = 0
x->s1 = 0
x->s2 = 0
x->tipo = 0
Return x
End Function
Sub fps_redefine(x As fps_t Ptr, op As Integer, y As fps_t Ptr, z As fps_t Ptr)
x->tipo = op
x->s1 = y
x->s2 = z
End Sub
Function _binary(x As fps_t Ptr, y As fps_t Ptr, op As Integer) As fps_t Ptr
Dim As fps_t Ptr s = fps_new()
s->s1 = x
s->s2 = y
s->tipo = op
Return s
End Function
Function _unary(x As fps_t Ptr, op As Integer) As fps_t Ptr
Dim As fps_t Ptr s = fps_new()
s->s1 = x
s->tipo = op
Return s
End Function
Function term(x As fps_t Ptr, n As Integer) As Single
Dim As Single ret = 0
Dim As Integer i
Select Case x->tipo
Case FPS_CONST
Return Iif(n > 0, 0, x->a0)
Case FPS_ADD
ret = term(x->s1, n) + term(x->s2, n)
Case FPS_SUB
ret = term(x->s1, n) - term(x->s2, n)
Case FPS_MUL
For i = 0 To n
ret += term(x->s1, i) * term(x->s2, n - i)
Next i
Case FPS_DIV
If term(x->s2, 0) = 0 Then Return 0 / 0 ' NaN
ret = term(x->s1, n)
For i = 1 To n
ret -= term(x->s2, i) * term(x, n - i) / term(x->s2, 0)
Next i
Case FPS_DERIV
ret = n * term(x->s1, n + 1)
Case FPS_INT
If n = 0 Then Return x->a0
ret = term(x->s1, n - 1) / n
Case Else
Print "Unknown operator "; x->tipo
End 1
End Select
Return ret
End Function
Function fps_const(a0 As Single) As fps_t Ptr
Dim As fps_t Ptr x = fps_new()
x->tipo = 0
x->a0 = a0
Return x
End Function
Dim As Integer i
Dim As fps_t Ptr one = fps_const(1)
Dim As fps_t Ptr fcos = fps_new() ' cosine
Dim As fps_t Ptr fsin = _unary(fcos, FPS_INT) ' sine
Dim As fps_t Ptr ftan = _binary(fsin, fcos, FPS_DIV) ' tangent
' redefine cos to complete the mutual recursion
fps_redefine(fcos, FPS_SUB, one, _unary(fsin, FPS_INT))
Dim As fps_t Ptr fexp = fps_const(1) ' exponential
' make exp recurse on self
fps_redefine(fexp, FPS_INT, fexp, 0)
Print "Sin:";
For i = 0 To 9
Print " "; term(fsin, i);
Next i
Print
Print "Cos:";
For i = 0 To 9
Print " "; term(fcos, i);
Next i
Print
Print "Tan:";
For i = 0 To 9
Print " "; term(ftan, i);
Next i
Print
Print "Exp:";
For i = 0 To 9
Print " "; term(fexp, i);
Next i
Print
Sleep</syntaxhighlight>
{{out}}
<pre>Sin: 0 1 0 -0.1666667 0 0.008333334 0 -0.0001984127 0 2.755732e-006
Cos: 1 0 -0.5 0 0.04166667 0 -0.001388889 0 2.480159e-005 0
Tan: 0 1 0 0.3333333 0 0.1333333 0 0.05396825 0 0.02186948
Exp: 1 1 0.5 0.1666667 0.04166667 0.008333334 0.001388889 0.0001984127 2.480159e-005 2.755732e-006</pre>
=={{header|Go}}==
<
import (
Line 824 ⟶ 1,133:
fmt.Println("sin:", partialSeries(sin))
fmt.Println("cos:", partialSeries(cos))
}</
Output:
<pre>
Line 830 ⟶ 1,139:
cos: 1.00000 0.00000 -0.50000 0.00000 0.04167 0.00000
</pre>
=={{header|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.
<
-- Invariant: coeffs must be an infinite list
Line 953 ⟶ 1,166:
cosx = 1 - int sinx
fiboS = 1 / fromFiniteList [1,-1,-1]</
Output:
Line 979 ⟶ 1,192:
J does not allow the definition of types. Also, J requires the programmer be lazy, rather than being lazy itself. With these "limitations" understood, here is an implementation in J:
<
divide=: [ +/ .*~ [:%.&.x: ] +/ .* Ai
diff=: 1 }. ] * i.@#
Line 985 ⟶ 1,198:
mult=: +//.@(*/)
plus=: +/@,:
minus=: -/@,:</
See also section 2 at http://www.jsoftware.com/papers/tot1.htm
Line 1,013 ⟶ 1,226:
See [[Formal power series/Java]]
=={{header|
{{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
<syntaxhighlight lang="jq">1/(1+.)</syntaxhighlight>
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:
<syntaxhighlight lang="jq">1/factorial</syntaxhighlight>
assuming "factorial" is defined in the usual way:
<syntaxhighlight lang="jq">def factorial:
reduce range(1; . + 1) as $i
(1; . * $i);</syntaxhighlight>
For ease of reference, we shall also define ps_exp as 1/factorial:
<syntaxhighlight lang="jq">def ps_exp: 1/factorial;</syntaxhighlight>
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:
<syntaxhighlight lang="jq">1 | ps_evaluate(ps_exp)</syntaxhighlight>
should evaluate to the number e approximately; using the version of ps_evaluate defined below, we find:
<syntaxhighlight lang="jq">1 | ps_evaluate(1/factorial)</syntaxhighlight>
evaluates to 2.7182818284590455.
The following function definitions are useful for other power series:
<syntaxhighlight lang="jq">def pow(n):
. as $x | n as $n
| reduce range(0;$n) as $i (1; . * $x);</syntaxhighlight>
For example, the power series 1 + Σ ( (x/i)^n ) where the summation is over i>0 can be written:
<syntaxhighlight lang="jq">1/pow(.)</syntaxhighlight>
The power series representation of ln(1 + x) is as follows:
<syntaxhighlight 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;
</syntaxhighlight>
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:
<syntaxhighlight lang="jq">def poly(ary): ary[.] // 0;</syntaxhighlight>
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:
<syntaxhighlight lang="jq">(poly([1,2,3]) + poly([-1,-2,-3]))</syntaxhighlight>
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====
<syntaxhighlight 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;
</syntaxhighlight>
====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:
<syntaxhighlight 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);</syntaxhighlight>
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.
<syntaxhighlight 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]);</syntaxhighlight>
====Examples====
<syntaxhighlight 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);
</syntaxhighlight>
''''cos == I(sin)''''
<syntaxhighlight 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</syntaxhighlight>
''''cos == 1 - I(sin)''''
<syntaxhighlight 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</syntaxhighlight>
=={{header|Julia}}==
'''Module''':
<syntaxhighlight lang="julia">module FormalPowerSeries
using Printf
import Base.iterate, Base.eltype, Base.one, Base.show, Base.IteratorSize
import Base.IteratorEltype, Base.length, Base.size, Base.convert
_div(a, b) = a / b
_div(a::Union{Integer,Rational}, b::Union{Integer,Rational}) = a // b
abstract type AbstractFPS{T<:Number} end
Base.IteratorSize(::AbstractFPS) = Base.IsInfinite()
Base.IteratorEltype(::AbstractFPS) = Base.HasEltype()
Base.eltype(::AbstractFPS{T}) where T = T
Base.one(::AbstractFPS{T}) where T = ConstantFPS(one(T))
function Base.show(io::IO, fps::AbstractFPS{T}) where T
itr = Iterators.take(fps, 8)
a, s = iterate(itr)
print(io, a)
a, s = iterate(itr, s)
@printf(io, " %s %s⋅x",
ifelse(sign(a) ≥ 0, '+', '-'), abs(a))
local i = 2
while (it = iterate(itr, s)) != nothing
a, s = it
@printf(io, " %s %s⋅x^%i",
ifelse(sign(a) ≥ 0, '+', '-'), abs(a), i)
i += 1
end
print(io, "...")
end
struct MinusFPS{T,A<:AbstractFPS{T}} <: AbstractFPS{T}
a::A
end
Base.:-(a::AbstractFPS{T}) where T = MinusFPS{T,typeof(a)}(a)
function Base.iterate(fps::MinusFPS)
v, s = iterate(fps.a)
return -v, s
end
function Base.iterate(fps::MinusFPS, st)
v, s = iterate(fps.a, st)
return -v, s
end
struct SumFPS{T,A<:AbstractFPS,B<:AbstractFPS} <: AbstractFPS{T}
a::A
b::B
end
Base.:+(a::AbstractFPS{A}, b::AbstractFPS{B}) where {A,B} =
SumFPS{promote_type(A, B),typeof(a),typeof(b)}(a, b)
Base.:-(a::AbstractFPS, b::AbstractFPS) = a + (-b)
function Base.iterate(fps::SumFPS{T,A,B}) where {T,A,B}
a1, s1 = iterate(fps.a)
a2, s2 = iterate(fps.b)
return T(a1 + a2), (s1, s2)
end
function Base.iterate(fps::SumFPS{T,A,B}, st) where {T,A,B}
stateA, stateB = st
valueA, stateA = iterate(fps.a, stateA)
valueB, stateB = iterate(fps.b, stateB)
return T(valueA + valueB), (stateA, stateB)
end
struct ProductFPS{T,A<:AbstractFPS,B<:AbstractFPS} <: AbstractFPS{T}
a::A
b::B
end
Base.:*(a::AbstractFPS{A}, b::AbstractFPS{B}) where {A,B} =
ProductFPS{promote_type(A, B),typeof(a),typeof(b)}(a, b)
function Base.iterate(fps::ProductFPS{T}) where T
a1, s1 = iterate(fps.a)
a2, s2 = iterate(fps.b)
T(sum(a1 .* a2)), (s1, s2, T[a1], T[a2])
end
function Base.iterate(fps::ProductFPS{T,A,B}, st) where {T,A,B}
stateA, stateB, listA, listB = st
valueA, stateA = iterate(fps.a, stateA)
valueB, stateB = iterate(fps.b, stateB)
push!(listA, valueA)
pushfirst!(listB, valueB)
return T(sum(listA .* listB)), (stateA, stateB, listA, listB)
end
struct DifferentiatedFPS{T,A<:AbstractFPS} <: AbstractFPS{T}
a::A
end
differentiate(fps::AbstractFPS{T}) where T = DifferentiatedFPS{T,typeof(fps)}(fps)
function Base.iterate(fps::DifferentiatedFPS{T,A}) where {T,A}
_, s = iterate(fps.a)
return Base.iterate(fps, (zero(T), s))
end
function Base.iterate(fps::DifferentiatedFPS{T,A}, st) where {T,A}
n, s = st
n += one(n)
v, s = iterate(fps.a, s)
return n * v, (n, s)
end
struct IntegratedFPS{T,A<:AbstractFPS} <: AbstractFPS{T}
a::A
k::T
end
integrate(fps::AbstractFPS{T}, k::T=zero(T)) where T = IntegratedFPS{T,typeof(fps)}(fps, k)
integrate(fps::AbstractFPS{T}, k::T=zero(T)) where T <: Integer =
IntegratedFPS{Rational{T},typeof(fps)}(fps, k)
function Base.iterate(fps::IntegratedFPS{T,A}, st=(0, 0)) where {T,A}
if st == (0, 0)
return fps.k, (one(T), 0)
end
n, s = st
if n == one(T)
v, s = iterate(fps.a)
else
v, s = iterate(fps.a, s)
end
r::T = _div(v, n)
n += one(n)
return r, (n, s)
end
# Examples of FPS: constant
struct FiniteFPS{T} <: AbstractFPS{T}
v::NTuple{N,T} where N
end
Base.iterate(fps::FiniteFPS{T}, st=1) where T =
st > lastindex(fps.v) ? (zero(T), st) : (fps.v[st], st + 1)
Base.convert(::Type{FiniteFPS}, x::Real) = FiniteFPS{typeof(x)}((x,))
FiniteFPS(r) = convert(FiniteFPS, r)
for op in (:+, :-, :*)
@eval Base.$op(x::Number, a::AbstractFPS) = $op(FiniteFPS(x), a)
@eval Base.$op(a::AbstractFPS, x::Number) = $op(a, FiniteFPS(x))
end
struct ConstantFPS{T} <: AbstractFPS{T}
k::T
end
Base.iterate(c::ConstantFPS, ::Any=nothing) = c.k, nothing
struct SineFPS{T} <: AbstractFPS{T} end
SineFPS() = SineFPS{Rational{Int}}()
function Base.iterate(::SineFPS{T}, st=(0, 1, 1)) where T
n, fac, s = st
local r::T
if iseven(n)
r = zero(T)
else
r = _div(one(T), (s * fac))
s = -s
end
n += 1
fac *= n
return r, (n, fac, s)
end
struct CosineFPS{T} <: AbstractFPS{T} end
CosineFPS() = CosineFPS{Rational{Int}}()
function Base.iterate(::CosineFPS{T}, st=(0, 1, 1)) where T
n, fac, s = st
local r::T
if iseven(n)
r = _div(one(T), (s * fac))
else
r = zero(T)
s = -s
end
n += 1
fac *= n
return r, (n, fac, s)
end
end # module FormalPowerSeries</syntaxhighlight>
'''Main''':
<syntaxhighlight lang="julia">using .FormalPowerSeries
@show cosine = FormalPowerSeries.CosineFPS()
@show sine = FormalPowerSeries.SineFPS()
intcosine = FormalPowerSeries.integrate(cosine)
intsine = FormalPowerSeries.integrate(sine)
uminintsine = 1 - FormalPowerSeries.integrate(sine)
# Check coefficients up to the 20th term
coefsine = collect(Iterators.take(sine, 20))
coefintcosine = collect(Iterators.take(intcosine, 20))
coefcosine = collect(Iterators.take(cosine, 20))
coefuminintsine = collect(Iterators.take(uminintsine, 20))
@assert coefsine == coefintcosine "The integral of cos should be sin"
@assert coefcosine == coefuminintsine "1 minus the integral of sin should be cos"
</syntaxhighlight>{{out}}
<pre>
cosine = FormalPowerSeries.CosineFPS() = 1//1 + 0//1⋅x - 1//2⋅x^2 + 0//1⋅x^3 + 1//24⋅x^4 + 0//1⋅x^5 - 1//720⋅x^6 + 0//1⋅x^7...
sine = FormalPowerSeries.SineFPS() = 0//1 + 1//1⋅x + 0//1⋅x^2 - 1//6⋅x^3 + 0//1⋅x^4 + 1//120⋅x^5 + 0//1⋅x^6 - 1//5040⋅x^7...
</pre>
=={{header|Kotlin}}==
This is a translation of the Java entry except that it uses fractions rather than double precision floating point numbers. The Frac class from the [[Arithmetic/Rational]] task has been embedded in the program for this purpose.
<syntaxhighlight lang="scala">// version 1.2.10
fun gcd(a: Long, b: Long): Long = if (b == 0L) a else gcd(b, a % b)
class Frac : Comparable<Frac> {
val num: Long
val denom: Long
companion object {
val ZERO = Frac(0, 1)
val ONE = Frac(1, 1)
}
constructor(n: Long, d: Long) {
require(d != 0L)
var nn = n
var dd = d
if (nn == 0L) {
dd = 1
}
else if (dd < 0) {
nn = -nn
dd = -dd
}
val g = Math.abs(gcd(nn, dd))
if (g > 1) {
nn /= g
dd /= g
}
num = nn
denom = dd
}
constructor(n: Int, d: Int) : this(n.toLong(), d.toLong())
operator fun plus(other: Frac) =
Frac(num * other.denom + denom * other.num, other.denom * denom)
operator fun unaryPlus() = this
operator fun unaryMinus() = Frac(-num, denom)
operator fun minus(other: Frac) = this + (-other)
operator fun times(other: Frac) =
Frac(this.num * other.num, this.denom * other.denom)
operator fun rem(other: Frac) = this - Frac((this / other).toLong(), 1) * other
operator fun inc() = this + ONE
operator fun dec() = this - ONE
fun inverse(): Frac {
require(num != 0L)
return Frac(denom, num)
}
operator fun div(other: Frac) = this * other.inverse()
fun abs() = if (num >= 0) this else -this
override fun compareTo(other: Frac): Int {
val diff = this.toDouble() - other.toDouble()
return when {
diff < 0.0 -> -1
diff > 0.0 -> +1
else -> 0
}
}
override fun equals(other: Any?): Boolean {
if (other == null || other !is Frac) return false
return this.compareTo(other) == 0
}
override fun hashCode() = num.hashCode() xor denom.hashCode()
override fun toString() = if (denom == 1L) "$num" else "$num/$denom"
fun toDouble() = num.toDouble() / denom
fun toLong() = num / denom
}
interface Gene {
fun coef(n: Int): Frac
}
class Term(private val gene: Gene) {
private val cache = mutableListOf<Frac>()
operator fun get(n: Int): Frac {
if (n < 0) return Frac.ZERO
if (n >= cache.size) {
for (i in cache.size..n) cache.add(gene.coef(i))
}
return cache[n]
}
}
class FormalPS {
private lateinit var term: Term
private companion object {
const val DISP_TERM = 12
const val X_VAR = "x"
}
constructor() {}
constructor(term: Term) {
this.term = term
}
constructor(polynomial: List<Frac>) :
this(Term(object : Gene {
override fun coef(n: Int) =
if (n < 0 || n >= polynomial.size)
Frac.ZERO
else
polynomial[n]
}))
fun copyFrom(other: FormalPS) {
term = other.term
}
fun inverseCoef(n: Int): Frac {
val res = Array(n + 1) { Frac.ZERO }
res[0] = term[0].inverse()
for (i in 1..n) {
for (j in 0 until i) res[i] += term[i - j] * res[j]
res[i] *= -res[0]
}
return res[n]
}
operator fun plus(other: FormalPS) =
FormalPS(Term(object : Gene {
override fun coef(n: Int) = term[n] + other.term[n]
}))
operator fun minus(other: FormalPS) =
FormalPS(Term(object : Gene {
override fun coef(n: Int) = term[n] - other.term[n]
}))
operator fun times(other: FormalPS) =
FormalPS(Term(object : Gene {
override fun coef(n: Int): Frac {
var res = Frac.ZERO
for (i in 0..n) res += term[i] * other.term[n - i]
return res
}
}))
operator fun div(other: FormalPS) =
FormalPS(Term(object : Gene {
override fun coef(n: Int): Frac {
var res = Frac.ZERO
for (i in 0..n) res += term[i] * other.inverseCoef(n - i)
return res
}
}))
fun diff() =
FormalPS(Term(object : Gene {
override fun coef(n: Int) = term[n + 1] * Frac(n + 1, 1)
}))
fun intg() =
FormalPS(Term(object : Gene {
override fun coef(n: Int) =
if (n == 0) Frac.ZERO else term[n - 1] * Frac(1, n)
}))
override fun toString() = toString(DISP_TERM)
private fun toString(dpTerm: Int): String {
val sb = StringBuilder()
var c = term[0]
if (c != Frac.ZERO) sb.append(c.toString())
for (i in 1 until dpTerm) {
c = term[i]
if (c != Frac.ZERO) {
if (c > Frac.ZERO && sb.length > 0) sb.append(" + ")
sb.append (when {
c == Frac.ONE -> X_VAR
c == -Frac.ONE -> " - $X_VAR"
c.num < 0 -> " - ${-c}$X_VAR"
else -> "$c$X_VAR"
})
if (i > 1) sb.append("^$i")
}
}
if (sb.length == 0) sb.append("0")
sb.append(" + ...")
return sb.toString()
}
}
fun main(args: Array<String>) {
var cos = FormalPS()
val sin = cos.intg()
cos.copyFrom(FormalPS(listOf(Frac.ONE)) - sin.intg())
println("SIN(x) = $sin")
println("COS(x) = $cos")
}</syntaxhighlight>
{{out}}
<pre>
SIN(x) = x - 1/6x^3 + 1/120x^5 - 1/5040x^7 + 1/362880x^9 - 1/39916800x^11 + ...
COS(x) = 1 - 1/2x^2 + 1/24x^4 - 1/720x^6 + 1/40320x^8 - 1/3628800x^10 + ...
</pre>
=={{header|Lua}}==
Parts of this depend on the formula for integration of a power series: integral(sum(a_n x^n)) = sum(a_n / n * x(n+1))
<
__add = function(z1, z2) return powerseries(function(n) return z1.coeff(n) + z2.coeff(n) end) end,
__sub = function(z1, z2) return powerseries(function(n) return z1.coeff(n) - z2.coeff(n) end) end,
Line 1,087 ⟶ 1,877:
tangent = sine / cosine
print(tangent(math.pi/3), tangent(math.pi/4), tangent(math.pi/6)) --something like 30000 function calls!</
=={{header|Mathematica}}/{{header|Wolfram Language}}==
Mathematica natively supports symbolic power series. For example, this input demonstrates that the integral of the series of Cos minus the series for sin is zero to the order of cancellation.
<syntaxhighlight lang="mathematica">cos = Series[Cos[x], {x, 0, 10}];
sin = Series[Sin[x], {x, 0, 8}];
sin - Integrate[cos, x]</syntaxhighlight>
{{out}}
<pre>O[x]^9</pre>
=={{header|Maxima}}==
<syntaxhighlight 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 + ... * /</syntaxhighlight>
=={{header|Nim}}==
===Using C algorithm===
The program is based on C algorithm but uses rationals instead of floats. There are many other differences in order to use Nim facilities (for instance, function and operator overloading) and some parts were borrowed from Kotlin version.
<syntaxhighlight lang="nim">import rationals, tables
type
Fraction = Rational[int]
FpsKind = enum fpsConst, fpsAdd, fpsSub, fpsMul, fpsDiv, fpsDeriv, fpsInteg
Fps = ref object
kind: FpsKind
s1, s2: Fps
a0: Fraction
cache: Table[Natural, Fraction]
const
Zero: Fraction = 0 // 1
One: Fraction = 1 // 1
DispTerm = 12
XVar = "x"
Super: array['0'..'9', string] = ["⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹"]
#---------------------------------------------------------------------------------------------------
proc `$`(fract: Fraction): string =
## Return the representation of a fraction without the denominator if it is equal to 1.
if fract.den == 1: $fract.num else: rationals.`$`(fract)
#---------------------------------------------------------------------------------------------------
proc exponent(n: Natural): string =
## Return the representation of an exponent using unicode superscript.
if n == 1: return ""
for d in $n: result.add(Super[d])
####################################################################################################
# FPS.
func newFps*(val = 0): Fps =
## Build a FPS of kind fpsConst using the given integer value.
Fps(kind: fpsConst, a0: val // 1)
#---------------------------------------------------------------------------------------------------
func newFps*(val: Fraction): Fps =
## Build a FPS of kind fpsConst using the given fraction.
Fps(kind: fpsConst, a0: val)
#---------------------------------------------------------------------------------------------------
func newFps*(op: FpsKind; x: Fps; y: Fps = nil): Fps =
## Build a FPS for a unary or binary operation.
Fps(kind: op, s1: x, s2: y, a0: Zero)
#---------------------------------------------------------------------------------------------------
func redefine*(fps: Fps; other: Fps) =
## Redefine a FPS, modifying its kind ans its operands.
fps.kind = other.kind
fps.s1 = other.s1
fps.s2 = other.s2
#---------------------------------------------------------------------------------------------------
## Operations on FPS.
func `+`*(x, y: Fps): Fps = newFps(fpsAdd, x, y)
func `-`*(x, y: Fps): Fps = newFps(fpsSub, x, y)
func `*`*(x, y: Fps): Fps = newFps(fpsMul, x, y)
func `/`*(x, y: Fps): Fps = newFps(fpsDiv, x, y)
func derivative*(x: Fps): Fps = newFps(fpsDeriv, x)
func integral*(x: Fps): Fps = newFps(fpsInteg, x)
#---------------------------------------------------------------------------------------------------
func `[]`*(fps: Fps; n: Natural): Fraction =
## Return the nth term of the FPS.
if n in fps.cache: return fps.cache[n]
case fps.kind
of fpsConst:
result = if n > 0: Zero else: fps.a0
of fpsAdd:
result = fps.s1[n] + fps.s2[n]
of fpsSub:
result = fps.s1[n] - fps.s2[n]
of fpsMul:
result = Zero
for i in 0..n: result += fps.s1[i] * fps.s2[n - i]
of fpsDiv:
let d = fps.s2[0]
if d == Zero: raise newException(DivByZeroDefect, "Division by null fraction")
result = fps.s1[n]
for i in 1..n: result -= fps.s2[i] * fps[n - i] / d
of fpsDeriv:
result = fps.s1[n + 1] * (n + 1)
of fpsInteg:
result = if n > 0: fps.s1[n - 1] / n else: fps.a0
fps.cache[n] = result
#---------------------------------------------------------------------------------------------------
proc `$`*(fps: Fps): string =
## Return the representation of a FPS.
var c = fps[0]
if c != Zero: result &= $c
for i in 1..<DispTerm:
c = fps[i]
if c != Zero:
if c > Zero:
if result.len > 0: result &= " + "
else:
result &= " - "
c = -c
result &= (if c == One: XVar else: $c & XVar) & exponent(i)
if result.len == 0: result &= '0'
result &= " + ..."
#———————————————————————————————————————————————————————————————————————————————————————————————————
# Build cos, sin and tan.
var cos = newFps()
let sin = cos.integral()
let tan = sin / cos
cos.redefine(newFps(1) - sin.integral())
echo "sin(x) = ", sin
echo "cos(x) = ", cos
echo "tan(x) = ", tan
# Check that derivative of sin is cos.
echo "derivative of sin(x) = ", sin.derivative()
# Build exp using recursion.
let exp = newFps()
exp.redefine(newFps(1) + exp.integral())
echo "exp(x) = ", exp</syntaxhighlight>
{{out}}
<pre>sin(x) = x - 1/6x³ + 1/120x⁵ - 1/5040x⁷ + 1/362880x⁹ - 1/39916800x¹¹ + ...
cos(x) = 1 - 1/2x² + 1/24x⁴ - 1/720x⁶ + 1/40320x⁸ - 1/3628800x¹⁰ + ...
tan(x) = x + 1/3x³ + 2/15x⁵ + 17/315x⁷ + 62/2835x⁹ + 1382/155925x¹¹ + ...
derivative of sin(x) = 1 - 1/2x² + 1/24x⁴ - 1/720x⁶ + 1/40320x⁸ - 1/3628800x¹⁰ + ...
exp(x) = 1 + x + 1/2x² + 1/6x³ + 1/24x⁴ + 1/120x⁵ + 1/720x⁶ + 1/5040x⁷ + 1/40320x⁸ + 1/362880x⁹ + 1/3628800x¹⁰ + 1/39916800x¹¹ + ...</pre>
===Using closure functions===
This version is largely inspired from Kotlin version even if the way to implement the algorithm is pretty different.
<syntaxhighlight lang="nim">import rationals, sequtils
type
Fraction = Rational[int]
# Function to compute coefficients.
CoeffFunc = proc(n: int): Fraction
# Formal power series.
Fps = ref object
cache: seq[Fraction] # Cache to store values.
coeffs: CoeffFunc # Function to compute coefficients.
const
Zero: Fraction = 0 // 1
One: Fraction = 1 // 1
XVar = "x"
DispTerm = 12
Super: array['0'..'9', string] = ["⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹"]
#---------------------------------------------------------------------------------------------------
proc `$`(fract: Fraction): string =
## Return the representation of a fraction without the denominator if it is equal to 1.
if fract.den == 1: $fract.num else: rationals.`$`(fract)
#---------------------------------------------------------------------------------------------------
proc exponent(n: Natural): string =
## Return the representation of an exponent using unicode superscript.
if n == 1: return ""
for d in $n: result.add(Super[d])
####################################################################################################
# FPS.
func newFps*(coeffs: CoeffFunc): Fps =
## Create a FPS using the given "coeffs" function.
Fps(coeffs: coeffs)
#---------------------------------------------------------------------------------------------------
func newFps*(coeffs: seq[Fraction]): Fps =
## Create a FPS using a list of fractions to initialize coefficients.
Fps(coeffs: proc(n: int): Fraction = (if n in 0..coeffs.high: coeffs[n] else: Zero))
#---------------------------------------------------------------------------------------------------
func newFps*(coeffs: seq[int]): Fps =
## Create a FPS using a list of integer values to initialize coefficients.
Fps(coeffs: proc(n: int): Fraction = (if n in 0..coeffs.high: coeffs[n] // 1 else: Zero))
#---------------------------------------------------------------------------------------------------
func copyFrom(dest, src: Fps) {.inline.} =
## Copy a FPS into another.
dest[] = src[]
#---------------------------------------------------------------------------------------------------
proc `[]`*(fps: Fps; n: int): Fraction =
## Return the element of degree "n" from a FPS.
if n < 0: return Zero
for i in fps.cache.len..n:
fps.cache.add(fps.coeffs(i))
result = fps.cache[n]
#---------------------------------------------------------------------------------------------------
proc inverseCoeff*(fps: FPS; n: int): Fraction =
## Return the inverse coefficient of coefficient of degree "n".
var res = repeat(Zero, n + 1)
res[0] = fps[0].reciprocal
for i in 1..n:
for j in 0..<i: res[i] += fps[i - j] * res[j]
res[i] *= -res[0]
result = res[n]
#---------------------------------------------------------------------------------------------------
proc `+`*(a, b: Fps): Fps =
## Build the FPS sum of two FPS.
Fps(coeffs: proc(n: int): Fraction = a[n] + b[n])
#---------------------------------------------------------------------------------------------------
proc `-`*(a, b: Fps): Fps =
## Build the FPS difference of two FPS.
Fps(coeffs: proc(n: int): Fraction = a[n] - b[n])
#---------------------------------------------------------------------------------------------------
proc `*`*(a, b: Fps): Fps =
## Build the FPS product of two FPS.
Fps(coeffs: proc(n: int): Fraction =
result = Zero
for i in 0..n: result += a[i] * b[n - i])
#---------------------------------------------------------------------------------------------------
proc `/`*(a, b: Fps): Fps =
## Build the FPS quotient of two FPS.
Fps(coeffs: proc(n: int): Fraction =
result = Zero
for i in 0..n: result += a[i] * b.inverseCoeff(n - i))
#---------------------------------------------------------------------------------------------------
proc derivative*(fps: Fps): Fps =
## Build the FPS derivative of a FPS.
Fps(coeffs: proc(n: int): Fraction = fps[n + 1] * (n + 1))
#---------------------------------------------------------------------------------------------------
proc integral*(fps: Fps): Fps =
## Build the FPS integral of a FPS.
Fps(coeffs: proc(n: int): Fraction = (if n == 0: Zero else: fps[n - 1] / n))
#---------------------------------------------------------------------------------------------------
proc `$`*(fps: Fps): string =
## Return the representation of a FPS.
var c = fps[0]
if c != Zero: result &= $c
for i in 1..<DispTerm:
c = fps[i]
if c != Zero:
if c > Zero:
if result.len > 0: result &= " + "
else:
result &= " - "
c = -c
result &= (if c == One: XVar else: $c & XVar) & exponent(i)
if result.len == 0: result &= '0'
result &= " + ..."
#———————————————————————————————————————————————————————————————————————————————————————————————————
# Build cos, sin and tan.
var cos = Fps()
let sin = cos.integral()
cos.copyFrom(newFps(@[1]) - sin.integral())
let tan = sin / cos
echo "sin(x) = ", sin
echo "cos(x) = ", cos
echo "tan(x) = ", tan
# Check that derivative of sin is cos.
echo "derivative of sin(x) = ", sin.derivative()
# Build exp using recursion.
let exp = Fps()
exp.copyFrom(newFps(@[1]) + exp.integral())
echo "exp(x) = ", exp</syntaxhighlight>
{{out}}
Same as with the other version.
=={{header|PARI/GP}}==
Uses the built-in power series handling. Change default(seriesprecision) to get more terms.
<
cos('x)</
=={{header|Perl
Although true Lazy Lists *can* be implemented using perl (using the builtin "tie" function), I felt that doing so would make the example harder to understand.
Instead, I chose to implement delayed evaluation with a generator function and a cache.
Creating a new arithmetic type in perl is relatively easy, using the "overload" module which comes with perl.
This was partly inspired by the Raku example, but is far from being a direct translation.
<syntaxhighlight lang="perl">
package FPS;
use strict;
use warnings;
use Math::BigRat;
sub new {
my
return bless {@_}, $class unless @_ == 1;
my $arg = shift;
return bless { more => $arg }, $class if 'CODE' eq ref $arg;
return bless { coeff => $arg }, $class if 'ARRAY' eq ref $arg;
bless { coeff => [$arg] }, $class;
}
sub coeff {
my $cache = ($self->{coeff} ||= []);
my $more = $self->{more};
for my $j ( @$cache .. $i ) {
last unless $more;
$cache->[$j] = $more->($j, $self);
}
$cache->[$i] or 0;
}
sub invert {
my
ref($orig)->new( sub {
my ($i, $self) = @_;
unless( $i ) {
my $a0 = $orig->coeff(0);
die "Cannot invert power series with zero constant term."
unless $a0;
(Math::BigRat->new(1) / $a0);
} else {
my $sum = 0;
my $terms = $self->{coeff};
for my $j (1 .. $i) {
$sum += $orig->coeff($j) * $terms->[$i - $j];
}
-$terms->[0] * $sum;
}
} );
}
sub fixargs {
my
my
$y = $class->new($y) unless UNIVERSAL::isa($y, $class);
($x, $y) = ($y, $x)
($class, $x, $y);
}
use overload '+' => sub {
my
$class->new(
}, '-' => sub {
my ($class, $x, $y) = &fixargs;
$class->new( sub { $x->coeff($_[0]) - $y->coeff($_[0]) } );
}, '*' => sub {
my ($class, $x, $y) = &fixargs;
$class->new( sub {
my $i = shift;
my $sum = 0;
$sum += $x->coeff($_) * $y->coeff($i-$_) for 0..$i;
$sum;
} );
}, '/' => sub {
my ($class, $x, $y) = &fixargs;
$x * $y->invert;
}, '""' => sub {
my $self = shift;
my $str = $self->coeff(0);
for my $i (1..10) {
my $c = $self->coeff($i);
next unless $c;
$str .= ($c < 0) ? (" - " . (-$c)) : (" + ".$c);
$str .= "x^$i";
}
$str;
};
sub differentiate {
ref($orig)->new( sub {
my $i = shift;
($i+1) * $orig->coeff($i);
} );
}
sub integrate {
ref($orig)->new( coeff => [0], more => sub {
my $i = shift;
$orig->coeff($i-1) / Math::BigRat->new($i);
} );
}
my $sin = __PACKAGE__->new;
my $cos = 1 - $sin->integrate;
%$sin = %{$cos->integrate};
my $tan = $sin / $cos;
my $exp = __PACKAGE__->new();
%$exp = (%{$exp->integrate}, coeff => [1]);
print "sin(x) ~= $sin\n";
print "cos(x) ~= $cos\n";
print "tan(x) ~= $tan\n";
print "exp(x) ~= $exp\n";
print "sin^2 + cos^2 = ", $sin*$sin + $cos*$cos, "\n";
1;
__END__
</syntaxhighlight>
{{out}}
<pre>
sin(x) ~= 0 + 1x^1 - 1/6x^3 + 1/120x^5 - 1/5040x^7 + 1/362880x^9
cos(x) ~= 1 - 1/2x^2 + 1/24x^4 - 1/720x^6 + 1/40320x^8 - 1/3628800x^10
tan(x) ~= 0 + 1x^1 + 1/3x^3 + 2/15x^5 + 17/315x^7 + 62/2835x^9
exp(x) ~= 1 + 1x^1 + 1/2x^2 + 1/6x^3 + 1/24x^4 + 1/120x^5 + 1/720x^6 + 1/5040x^7 + 1/40320x^8 + 1/362880x^9 + 1/3628800x^10
sin^2 + cos^2 = 1
</pre>
For a version which *does* use proper lazy lists, see [[Formal power series/Perl]]
=={{header|Phix}}==
{{trans|C}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">enum</span> <span style="color: #000000;">FPS_UNDEF</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">FPS_CONST</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">FPS_ADD</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">FPS_SUB</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">FPS_MUL</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">FPS_DIV</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">FPS_DERIV</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">FPS_INT</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">FPS_MAX</span><span style="color: #0000FF;">=$</span>
<span style="color: #008080;">type</span> <span style="color: #000000;">fps_type</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">FPS_UNDEF</span> <span style="color: #008080;">and</span> <span style="color: #000000;">f</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">FPS_MAX</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">type</span>
<span style="color: #008080;">enum</span> <span style="color: #000000;">FPS_TYPE</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">FPS_S1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">FPS_S2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">FPS_A0</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">fpss</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
<span style="color: #008080;">type</span> <span style="color: #000000;">fps</span><span style="color: #0000FF;">(</span><span style="color: #004080;">object</span> <span style="color: #000000;">id</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #004080;">integer</span><span style="color: #0000FF;">(</span><span style="color: #000000;">id</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">and</span> <span style="color: #000000;">id</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">1</span> <span style="color: #008080;">and</span> <span style="color: #000000;">id</span><span style="color: #0000FF;"><=</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fpss</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">type</span>
<span style="color: #008080;">type</span> <span style="color: #000000;">fpsn</span><span style="color: #0000FF;">(</span><span style="color: #004080;">object</span> <span style="color: #000000;">id</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">id</span><span style="color: #0000FF;">=</span><span style="color: #004600;">NULL</span> <span style="color: #008080;">or</span> <span style="color: #000000;">fps</span><span style="color: #0000FF;">(</span><span style="color: #000000;">id</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">type</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">fps_new</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fps_type</span> <span style="color: #000000;">ft</span><span style="color: #0000FF;">=</span><span style="color: #000000;">FPS_UNDEF</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fpsn</span> <span style="color: #000000;">s1</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s2</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">a0</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">fpss</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fpss</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">ft</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a0</span><span style="color: #0000FF;">})</span>
<span style="color: #000000;">fps</span> <span style="color: #000000;">fpsid</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fpss</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">fpsid</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #000080;font-style:italic;">-- as per C, for (eg) self or mutually recursive definitions.</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">fps_redefine</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fps</span> <span style="color: #000000;">fpsid</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fps_type</span> <span style="color: #000000;">ft</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fpsn</span> <span style="color: #000000;">s1id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s2id</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">object</span> <span style="color: #000000;">a0</span><span style="color: #0000FF;">=</span><span style="color: #008000;">""</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">fpss</span><span style="color: #0000FF;">[</span><span style="color: #000000;">fpsid</span><span style="color: #0000FF;">][</span><span style="color: #000000;">FPS_TYPE</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ft</span>
<span style="color: #000000;">fpss</span><span style="color: #0000FF;">[</span><span style="color: #000000;">fpsid</span><span style="color: #0000FF;">][</span><span style="color: #000000;">FPS_S1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s1id</span>
<span style="color: #000000;">fpss</span><span style="color: #0000FF;">[</span><span style="color: #000000;">fpsid</span><span style="color: #0000FF;">][</span><span style="color: #000000;">FPS_S2</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s2id</span>
<span style="color: #008080;">if</span> <span style="color: #004080;">atom</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a0</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">fpss</span><span style="color: #0000FF;">[</span><span style="color: #000000;">fpsid</span><span style="color: #0000FF;">][</span><span style="color: #000000;">FPS_A0</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">a0</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">fps_const</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">a0</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">fps</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">fps_new</span><span style="color: #0000FF;">(</span><span style="color: #000000;">FPS_CONST</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a0</span><span style="color: #0000FF;">:=</span><span style="color: #000000;">a0</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- (aside: in the above, the ":=a0" refers to the local namespace
-- as usual, whereas "a0:=" refers to the param namespace
-- /inside/ the () of fps_new(), so there is no conflict.)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">x</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">INF</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1e300</span><span style="color: #0000FF;">*</span><span style="color: #000000;">1e300</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">NAN</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-(</span><span style="color: #000000;">INF</span><span style="color: #0000FF;">/</span><span style="color: #000000;">INF</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">/* Taking the n-th term of series. This is where actual work is done. */</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fps</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">int</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">ret</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">fps_type</span> <span style="color: #000000;">ft</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fpsn</span> <span style="color: #000000;">s1id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fpsn</span> <span style="color: #000000;">s2id</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">a0</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">fpss</span><span style="color: #0000FF;">[</span><span style="color: #000000;">x</span><span style="color: #0000FF;">]</span>
<span style="color: #000080;font-style:italic;">-- FPS_TYPE, FPS_S1, FPS_S2, FPS_A0 <-- nb above must match</span>
<span style="color: #008080;">switch</span> <span style="color: #000000;">ft</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">case</span> <span style="color: #000000;">FPS_CONST</span><span style="color: #0000FF;">:</span> <span style="color: #000000;">ret</span> <span style="color: #0000FF;">:=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">></span><span style="color: #000000;">0</span> <span style="color: #0000FF;">?</span> <span style="color: #000000;">0</span> <span style="color: #0000FF;">:</span> <span style="color: #000000;">a0</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">case</span> <span style="color: #000000;">FPS_ADD</span><span style="color: #0000FF;">:</span> <span style="color: #000000;">ret</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s1id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s2id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">case</span> <span style="color: #000000;">FPS_SUB</span><span style="color: #0000FF;">:</span> <span style="color: #000000;">ret</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s1id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s2id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">case</span> <span style="color: #000000;">FPS_MUL</span><span style="color: #0000FF;">:</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">ret</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s1id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s2id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">case</span> <span style="color: #000000;">FPS_DIV</span><span style="color: #0000FF;">:</span>
<span style="color: #008080;">if</span> <span style="color: #008080;">not</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s2id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">NAN</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">ret</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s1id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">ret</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s2id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s2id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">case</span> <span style="color: #000000;">FPS_DERIV</span><span style="color: #0000FF;">:</span> <span style="color: #000000;">ret</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s1id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">case</span> <span style="color: #000000;">FPS_INT</span><span style="color: #0000FF;">:</span> <span style="color: #000000;">ret</span> <span style="color: #0000FF;">:=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #0000FF;">?</span> <span style="color: #000000;">a0</span> <span style="color: #0000FF;">:</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s1id</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">default</span><span style="color: #0000FF;">:</span> <span style="color: #000000;">ret</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #000080;font-style:italic;">-- (fatal error)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">switch</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">ret</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">term9</span><span style="color: #0000FF;">(</span><span style="color: #004080;">string</span> <span style="color: #000000;">txt</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fps</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s:"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">txt</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">to</span> <span style="color: #000000;">9</span> <span style="color: #008080;">do</span> <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" %g"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">))</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<span style="color: #000000;">fps</span> <span style="color: #000000;">one</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">fps_const</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">fps</span> <span style="color: #000000;">fcos</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">fps_new</span><span style="color: #0000FF;">()</span> <span style="color: #000080;font-style:italic;">/* cosine */</span>
<span style="color: #000000;">fps</span> <span style="color: #000000;">fsin</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">fps_new</span><span style="color: #0000FF;">(</span><span style="color: #000000;">FPS_INT</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fcos</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">/* sine */</span>
<span style="color: #000000;">fps</span> <span style="color: #000000;">ftan</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">fps_new</span><span style="color: #0000FF;">(</span><span style="color: #000000;">FPS_DIV</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fsin</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fcos</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">/* tangent */
/* redefine cos to complete the mutual recursion */</span>
<span style="color: #000000;">fps_redefine</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fcos</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">FPS_SUB</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">one</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fps_new</span><span style="color: #0000FF;">(</span><span style="color: #000000;">FPS_INT</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fsin</span><span style="color: #0000FF;">))</span>
<span style="color: #000000;">fps</span> <span style="color: #000000;">fexp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">fps_const</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">);</span> <span style="color: #000080;font-style:italic;">/* exponential */
/* make exp recurse on self */</span>
<span style="color: #000000;">fps_redefine</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fexp</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">FPS_INT</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fexp</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">);</span>
<span style="color: #000000;">term9</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"Sin"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fsin</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">term9</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"Cos"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fcos</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">term9</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"Tan"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ftan</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">term9</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"Exp"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fexp</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Sin: 0 1 0 -0.166667 0 0.008333 0 -0.000198 0 2.75573e-6
Cos: 1 0 -0.5 0 0.041667 0 -0.001389 0 2.48016e-5 0
Tan: 0 1 0 0.333333 0 0.133333 0 0.053968 0 0.021869
Exp: 1 1 0.5 0.166667 0.041667 0.008333 0.001389 0.000198 2.48016e-5 2.75573e-6
</pre>
=={{header|PicoLisp}}==
With a 'lazy' function, as a frontend to '[http://software-lab.de/doc/refC.html#cache cache]',
<
(def (car Args)
(list (cadr Args)
(cons 'cache (lit (cons))
(cddr Args) ) ) ) )</
we can build a formal power series functionality:
<
(de fpsOne (N)
Line 1,256 ⟶ 2,557:
(if (=0 N)
1.0
((fpsIntegrate fpsExp) N) ) )</
Test:
<
(for N (range 1 11 2)
(prin " " (round (fpsSin N) 9)) )
Line 1,276 ⟶ 2,577:
(for N (range 0 6)
(prin " " (round (fpsExp N) 7)) )
(prinl)</
Output:
<pre>SIN: 1.000000000 -0.166666667 0.008333333 -0.000198413 0.000002756 -0.000000025
Line 1,285 ⟶ 2,586:
=={{header|Python}}==
{{works with|Python|2.6, 3.x}}
<
For a discussion on pipe() and head() see
http://paddy3118.blogspot.com/2009/05/pipe-fitting-with-python-generators.html
Line 1,380 ⟶ 2,681:
assert s == integc, "The integral of cos should be sin"
assert c == integs1, "1 minus the integral of sin should be cos"</
'''Sample output'''
Line 1,392 ⟶ 2,693:
Alternate version that uses a generator function to allow sine and cosine to be defined recursively, following the same method as [[Hamming numbers#Alternate version using "Cyclic Iterators"]]:
{{works with|Python|2.6, 3.x}}
<
from fractions import Fraction
try:
Line 1,449 ⟶ 2,750:
print(list(islice(sinx, 10)))
print("sine")
print(list(islice(cosx, 10)))</
'''Sample output'''
Line 1,458 ⟶ 2,759:
[1, Fraction(0, 1), Fraction(-1, 2), Fraction(0, 1), Fraction(1, 24), Fraction(0, 1), Fraction(-1, 720), Fraction(0, 1), Fraction(1, 40320), Fraction(0, 1)]
</pre>
===Operator overloading===
Define an iterator class as a polynomial to provide overloaded operators and automatic tee-ing. It is kind of overkill.
<syntaxhighlight lang="python">from itertools import count, chain, tee, islice, cycle
from fractions import Fraction
# infinite polynomial class
class Poly:
def __init__(self, gen = None):
self.gen, self.source = (None, gen) if type(gen) is Poly \
else (gen, None)
def __iter__(self):
# We're essentially tee'ing it everytime the iterator
# is, well, iterated. This may be excessive.
return Poly(self)
def getsource(self):
if self.gen == None:
s = self.source
s.getsource()
(a,b) = tee(s.gen, 2)
s.gen = a
self.gen = b
def next(self):
self.getsource()
return next(self.gen)
__next__ = next
# Overload "<<" as stream input operator. Hey, C++ does it.
def __lshift__(self, a): self.gen = a
# The other operators are pretty much what one would expect
def __neg__(self): return Poly(-x for x in self)
def __sub__(a, b): return a + (-b)
def __rsub__(a, n):
a = Poly(a)
def gen():
yield(n - next(a))
for x in a: yield(-x)
return Poly(gen())
def __add__(a, b):
if type(b) is Poly:
return Poly(x + y for (x,y) in zip(a,b))
a = Poly(a)
def gen():
yield(next(a) + b)
for x in a: yield(x)
return Poly(gen())
def __radd__(a,b):
return a + b
def __mul__(a,b):
if not type(b) is Poly:
return Poly(x*b for x in a)
def gen():
s = Poly(cycle([0]))
for y in b:
s += y*a
yield(next(s))
return Poly(gen())
def __rmul__(a,b): return a*b
def __truediv__(a,b):
if not type(b) is Poly:
return Poly(Fraction(x, b) for x in a)
a, b = Poly(a), Poly(b)
def gen():
r, bb = a,next(b)
while True:
aa = next(r)
q = Fraction(aa, bb)
yield(q)
r -= q*b
return Poly(gen())
# these two would probably be better as class methods
def inte(a):
def gen():
yield(0)
for (x,n) in zip(a, count(1)):
yield(Fraction(x,n))
return Poly(gen())
def diff(a):
def gen():
for (x, n) in zip(a, count(0)):
if n: yield(x*n)
return Poly(gen())
# all that for the syntactic sugar
sinx, cosx, tanx, expx = Poly(), Poly(), Poly(), Poly()
sinx << inte(cosx)
cosx << 1 - inte(sinx)
tanx << sinx / cosx # "=" would also work here
expx << 1 + inte(expx)
for n,x in zip(("sin", "cos", "tan", "exp"), (sinx, cosx, tanx, expx)):
print(n, ', '.join(map(str, list(islice(x, 10)))))</syntaxhighlight>
=={{header|Racket}}==
Using '''Lazy Racket''':
<syntaxhighlight lang="racket">
#lang lazy
(require racket/match)
;; element-wise addition and subtraction
(define (<+> s1 s2) (map + s1 s2))
(define (<-> s1 s2) (map - s1 s2))
;; element-wise scaling
(define (scale a s) (map (λ (x) (* a x)) s))
;; series multiplication
(define (<*> fs gs)
(match-let ([(cons f ft) (! fs)]
[(cons g gt) (! gs)])
(cons (* f g) (<+> (scale f gt) (<*> ft gs)))))
;; series division
(define (</> fs gs)
(match-letrec ([(cons f ft) (! fs)]
[(cons g gt) (! gs)]
[qs (cons (/ f g) (scale (/ g) (<-> ft (<*> qs gt))))])
qs))
;; integration and differentiation
(define (int f) (map / f (enum 1)))
(define (diff f) (map * (cdr f) (enum 1)))
;; series of natural numbers greater then n
(define (enum n) (cons n (enum (+ 1 n ))))
</syntaxhighlight>
Examples:
<syntaxhighlight lang="racket">
(define <sin> (cons 0 (int <cos>)))
(define <cos> (cons 1 (scale -1 (int <sin>))))
-> (!! (take 10 <sin>))
'(0 1 0 -1/6 0 1/120 0 -1/5040 0 1/362880)
-> (!! (take 10 <cos>))
'(1 0 -1/2 0 1/24 0 -1/720 0 1/40320 0)
-> (!! (take 10 (diff <sin>)))
'(1 0 -1/2 0 1/24 0 -1/720 0 1/40320 0)
; sin(x)² + cos(x)² = 1
-> (!! (take 10 (<+> (<*> <cos> <cos>) (<*> <sin> <sin>))))
'(1 0 0 0 0 0 0 0 0 0)
; series of (tan x)
-> (!! (take 10 (</> <sin> <cos>)))
'(0 1 0 1/3 0 2/15 0 17/315 0 62/2835)
</syntaxhighlight>
=={{header|Raku}}==
(formerly Perl 6)
<syntaxhighlight lang="raku" line>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);</syntaxhighlight>
{{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|Ruby}}==
{{works with|Ruby|2.3}}
'''Implementation'''
<br />
This uses Ruby's Lazy Enumerators to generate the indefinite length sequences.
<br />
Implements the '<code>+</code>', '<code>-</code>', '<code>*</code>', and
'<code>/</code>' operators, and the <code>deriv()</code> and
<code>integ()</code> methods. Also implements coercion and conversion with
standard Numeric types, and formatting as a string for display. Provides for
evaluating the series at a point.
<br />
Inspirations from the Schame and Raku implementations.
<syntaxhighlight lang="ruby"># Class implementing the Formal Power Series type.
class Fps
# Initialize the FPS instance.
# When nothing specified, all coefficients are 0.
# When const: specifies n, all coefficients are n.
# When delta: specifies n, a[0] = n, then all higher coefficients are zero.
# When iota: specifies n, coefficients are consecutive integers, beginning with a[0] = n.
# When init: specifies an array, coefficients are the array elements, padded with zeroes.
# When enum: specifies a lazy enumerator, that is used for the internal coefficients enum.
def initialize(const: nil, delta: nil, iota: nil, init: nil, enum: nil)
# Create (or save) the specified coefficient enumerator.
case
when const
@coeffenum = make_const(const)
when delta
@coeffenum = make_delta(delta)
when iota
@coeffenum = make_iota(iota)
when init
@coeffenum = make_init(init)
when enum
@coeffenum = enum
else
@coeffenum = make_const(0)
end
# Extend the coefficient enumerator instance with an element accessor.
@coeffenum.instance_eval do
def [](index)
self.drop(index).first
end
end
end
# Return the coefficient at the given index.
def [](index)
@coeffenum.drop(index).first
end
# Return sum: this FPS plus the given FPS.
def +(other)
other = convert(other)
Fps.new(enum:
Enumerator.new do |yielder, inx: 0|
loop do
yielder.yield(@coeffenum[inx] + other[inx])
inx += 1
end
end.lazy)
end
# Return difference: this FPS minus the given FPS.
def -(other)
other = convert(other)
Fps.new(enum:
Enumerator.new do |yielder, inx: 0|
loop do
yielder.yield(@coeffenum[inx] - other[inx])
inx += 1
end
end.lazy)
end
# Return product: this FPS multiplied by the given FPS.
def *(other)
other = convert(other)
Fps.new(enum:
Enumerator.new do |yielder, inx: 0|
loop do
coeff = (0..inx).reduce(0) { |sum, i| sum + (@coeffenum[i] * other[inx - i]) }
yielder.yield(coeff)
inx += 1
end
end.lazy)
end
# Return quotient: this FPS divided by the given FPS.
def /(other)
other = convert(other)
Fps.new(enum:
Enumerator.new do |yielder, inx: 1|
coeffs = [ Rational(@coeffenum[0], other[0]) ]
yielder.yield(coeffs[-1])
loop do
coeffs <<
Rational(
@coeffenum[inx] -
(1..inx).reduce(0) { |sum, i| sum + (other[i] * coeffs[inx - i]) },
other[0])
yielder.yield(coeffs[-1])
inx += 1
end
end.lazy)
end
# Return the derivative of this FPS.
def deriv()
Fps.new(enum:
Enumerator.new do |yielder, inx: 0|
iota = Fps.new(iota: 1)
loop do
yielder.yield(@coeffenum[inx + 1] * iota[inx])
inx += 1
end
end.lazy)
end
# Return the integral of this FPS.
def integ()
Fps.new(enum:
Enumerator.new do |yielder, inx: 0|
iota = Fps.new(iota: 1)
yielder.yield(Rational(0, 1))
loop do
yielder.yield(Rational(@coeffenum[inx], iota[inx]))
inx += 1
end
end.lazy)
end
# Assign a new value to an existing FPS instance.
def assign(other)
other = convert(other)
@coeffenum = other.get_enum
end
# Coerce a Numeric into an FPS instance.
def coerce(other)
if other.kind_of?(Numeric)
[ Fps.new(delta: other), self ]
else
raise TypeError 'non-numeric can\'t be coerced into FPS type'
end
end
# Convert to Integer. (Truncates to 0th coefficient.)
def to_i()
@coeffenum[0].to_i
end
# Convert to Float. (Truncates to 0th coefficient.)
def to_f()
@coeffenum[0].to_f
end
# Convert to Rational. (Truncates to 0th coefficient.)
def to_r()
@coeffenum[0].to_r
end
# Convert to String the first count terms of the FPS.
def to_s(count = 0)
if count <= 0
super()
else
retstr = ''
count.times do |inx|
coeff = (@coeffenum[inx].to_r.denominator == 1) ? @coeffenum[inx].to_i : @coeffenum[inx]
if !(coeff.zero?)
prefix = (retstr != '') ? ' ' : ''
coeffstr =
((coeff.abs == 1) && (inx != 0)) ? '' : "#{coeff.abs.to_s}#{(inx == 0) ? '' : '*'}"
suffix = (inx == 0) ? '' : (inx == 1) ? 'x' : "x^#{inx}"
if coeff < 0
prefix << ((retstr != '') ? '- ' : '-')
else
prefix << ((retstr != '') ? '+ ' : '')
end
retstr << "#{prefix}#{coeffstr}#{suffix}"
end
end
(retstr == '') ? '0' : retstr
end
end
# Evaluate this FPS at the given x value to the given count of terms.
def eval(x, count)
@coeffenum.first(count).each_with_index.reduce(0) { |sum, (coeff, inx) | sum + coeff * x**inx }
end
# Forward method calls to the @coeffenum instance.
def method_missing(name, *args, &block)
@coeffenum.send(name, *args, &block)
end
# Forward respond_to? to the @coeffenum instance.
def respond_to_missing?(name, incl_priv)
@coeffenum.respond_to?(name, incl_priv)
end
protected
# Return reference to the underlying coefficient enumeration.
def get_enum()
@coeffenum
end
private
# Create a "const" lazy enumerator with the given n.
# All elements are n.
def make_const(n)
Enumerator.new do |yielder|
loop { yielder.yield(n) }
end.lazy
end
# Create a "delta" lazy enumerator with the given n.
# First element is n, then all subsequent elements are zero.
def make_delta(n)
Enumerator.new do |yielder|
yielder.yield(n)
loop { yielder.yield(0) }
end.lazy
end
# Create an "iota" lazy enumerator with the given n.
# Elements are consecutive integers, beginning with n.
def make_iota(n)
Enumerator.new do |yielder, i: n|
loop { yielder.yield(i); i += 1 }
end.lazy
end
# Create an "init" lazy enumerator with the given array.
# Elements are the array elements, padded with zeroes.
def make_init(array)
Enumerator.new do |yielder, inx: -1|
loop { yielder.yield((inx < (array.length - 1)) ? array[inx += 1] : 0) }
end.lazy
end
# Convert a Numeric to an FPS instance, if needed.
def convert(other)
if other.kind_of?(Fps)
other
else
if other.kind_of?(Numeric)
Fps.new(delta: other)
else
raise TypeError 'non-numeric can\'t be converted to FPS type'
end
end
end
end</syntaxhighlight>
'''Testing'''
<syntaxhighlight lang="ruby">puts
puts('FPS Creation:')
a = Fps.new(iota: 3)
puts("iota: a = #{a.first(10)}")
b = Fps.new(init: [12, 24, 36, 48])
puts("init: b = #{b.first(10)}")
puts
puts('FPS Arithmetic:')
puts("(a + b) = #{(a + b).first(10)}")
puts("(a - b) = #{(a - b).first(10)}")
puts("(a * b) = #{(a * b).first(10)}")
puts("(a / b) = #{(a / b).first(10)}")
puts("((a + b) - b) = #{((a + b) - b).first(10)}")
puts("((a / b) * b) = #{((a / b) * b).first(10)}")
puts
puts('FPS w/ Other Numerics:')
puts("(a + 3) = #{(a + 3).first(10)}")
puts("(a - 3) = #{(a - 3).first(10)}")
puts("(a * 3) = #{(a * 3).first(10)}")
puts("(a / 3) = #{(a / 3).first(10)}")
puts("(3 + a) = #{(3 + a).first(10)}")
puts("(3 - a) = #{(3 - a).first(10)}")
puts("(3 * a) = #{(3 * a).first(10)}")
puts("(3 / a) = #{(3 / a).first(10)}")
puts("(a + 4.4) = #{(a + 4.5).first(10)}")
puts("(a + Rational(11, 3)) = #{(a + Rational(11, 3)).first(10)}")
puts("(4.4 + a) = #{(4.5 + a).first(10)}")
puts("(Rational(11, 3) + a) = #{(Rational(11, 3) + a).first(10)}")
puts
puts('FPS Differentiation and Integration:')
puts("b.deriv = #{b.deriv.first(10)}")
puts("b.integ = #{b.integ.first(10)}")
puts
puts('Define sin(x) and cos(x) FPSs in terms of each other:')
fpssin = Fps.new
fpscos = 1 - fpssin.integ
fpssin.assign(fpscos.integ)
puts("fpssin = #{fpssin.first(10)}")
puts("fpscos = #{fpscos.first(10)}")
puts
puts('Display sin(x) and cos(x) FPSs as strings:')
puts("sin(x) = #{fpssin.to_s(10)}")
puts("cos(x) = #{fpscos.to_s(10)}")
puts
puts('Define tan(x) FPS as sin(x) / cos(x) from above:')
fpstan = fpssin / fpscos
puts("tan(x) = #{fpstan.to_s(10)}")
puts
puts('Compute sin^2(x)+cos^2(x) FPS from above:')
puts("sin^2(x)+cos^2(x) = #{((fpssin * fpssin) + (fpscos * fpscos)).to_s(5)}")
puts
puts('Define exp(x) in terms of its own integral:')
fpsexp = Fps.new
fpsexp.assign(1 + fpsexp.integ)
puts("exp(x) = #{fpsexp.to_s(10)}")
puts
puts('Evaluate the above at a few points using a few terms:')
puts("sin(0) = #{fpssin.eval(0, 8).to_f}")
puts("cos(0) = #{fpscos.eval(0, 8).to_f}")
puts("sin(pi/2) = #{fpssin.eval(Math::PI / 2, 16).to_f}")
puts("cos(pi/2) = #{fpscos.eval(Math::PI / 2, 16).to_f}")
puts("sin(pi) = #{fpssin.eval(Math::PI, 16).to_f}")
puts("cos(pi) = #{fpscos.eval(Math::PI, 16).to_f}")
puts("tan(0) = #{fpstan.eval(0, 8).to_f}")
puts("exp(1) = #{fpsexp.eval(1, 14).to_f}")</syntaxhighlight>
{{out}}
<pre>FPS Creation:
iota: a = [3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
init: b = [12, 24, 36, 48, 0, 0, 0, 0, 0, 0]
FPS Arithmetic:
(a + b) = [15, 28, 41, 54, 7, 8, 9, 10, 11, 12]
(a - b) = [-9, -20, -31, -42, 7, 8, 9, 10, 11, 12]
(a * b) = [36, 120, 264, 480, 600, 720, 840, 960, 1080, 1200]
(a / b) = [(1/4), (-1/6), (0/1), (0/1), (5/4), (-11/6), (2/3), (0/1), (25/4), (-85/6)]
((a + b) - b) = [3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
((a / b) * b) = [(3/1), (4/1), (5/1), (6/1), (7/1), (8/1), (9/1), (10/1), (11/1), (12/1)]
FPS w/ Other Numerics:
(a + 3) = [6, 4, 5, 6, 7, 8, 9, 10, 11, 12]
(a - 3) = [0, 4, 5, 6, 7, 8, 9, 10, 11, 12]
(a * 3) = [9, 12, 15, 18, 21, 24, 27, 30, 33, 36]
(a / 3) = [(1/1), (4/3), (5/3), (2/1), (7/3), (8/3), (3/1), (10/3), (11/3), (4/1)]
(3 + a) = [6, 4, 5, 6, 7, 8, 9, 10, 11, 12]
(3 - a) = [0, -4, -5, -6, -7, -8, -9, -10, -11, -12]
(3 * a) = [9, 12, 15, 18, 21, 24, 27, 30, 33, 36]
(3 / a) = [(1/1), (-4/3), (1/9), (2/27), (4/81), (8/243), (16/729), (32/2187), (64/6561), (128/19683)]
(a + 4.4) = [7.5, 4, 5, 6, 7, 8, 9, 10, 11, 12]
(a + Rational(11, 3)) = [(20/3), 4, 5, 6, 7, 8, 9, 10, 11, 12]
(4.4 + a) = [7.5, 4, 5, 6, 7, 8, 9, 10, 11, 12]
(Rational(11, 3) + a) = [(20/3), 4, 5, 6, 7, 8, 9, 10, 11, 12]
FPS Differentiation and Integration:
b.deriv = [24, 72, 144, 0, 0, 0, 0, 0, 0, 0]
b.integ = [(0/1), (12/1), (12/1), (12/1), (12/1), (0/1), (0/1), (0/1), (0/1), (0/1)]
Define sin(x) and cos(x) FPSs in terms of each other:
fpssin = [(0/1), (1/1), (0/1), (-1/6), (0/1), (1/120), (0/1), (-1/5040), (0/1), (1/362880)]
fpscos = [(1/1), (0/1), (-1/2), (0/1), (1/24), (0/1), (-1/720), (0/1), (1/40320), (0/1)]
Display sin(x) and cos(x) FPSs as strings:
sin(x) = x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + 1/362880*x^9
cos(x) = 1 - 1/2*x^2 + 1/24*x^4 - 1/720*x^6 + 1/40320*x^8
Define tan(x) FPS as sin(x) / cos(x) from above:
tan(x) = x + 1/3*x^3 + 2/15*x^5 + 17/315*x^7 + 62/2835*x^9
Compute sin^2(x)+cos^2(x) FPS from above:
sin^2(x)+cos^2(x) = 1
Define exp(x) in terms of its own integral:
exp(x) = 1 + x + 1/2*x^2 + 1/6*x^3 + 1/24*x^4 + 1/120*x^5 + 1/720*x^6 + 1/5040*x^7 + 1/40320*x^8 + 1/362880*x^9
Evaluate the above at a few points using a few terms:
sin(0) = 0.0
cos(0) = 1.0
sin(pi/2) = 0.9999999999939768
cos(pi/2) = -6.513362355903657e-11
sin(pi) = -7.727858894168152e-07
cos(pi) = -1.0000041678091434
tan(0) = 0.0
exp(1) = 2.718281828446759</pre>
=={{header|Scheme}}==
Definitions of operations on lazy lists:
<
(syntax-rules ()
((_ lar ldr) (delay (cons lar (delay ldr))))))
Line 1,483 ⟶ 3,441:
(define (repeat n)
(lons n (repeat n)))</
Definitions of operations on formal power series:
<
(apply lap + llists))
Line 1,515 ⟶ 3,473:
(define (fpsdif llist)
(lap * (iota 1) (ldr llist)))</
Now the sine and cosine functions can be defined in terms of eachother using integrals:
<
(fps- (lons 1 (repeat 0)) (fpsint (delay (force fpssin)))))
Line 1,527 ⟶ 3,485:
(display (take 10 fpscos))
(newline)</
Output:
(0 1 0 -1/6 0 1/120 0 -1/5040 0 1/362880)
(1 0 -1/2 0 1/24 0 -1/720 0 1/40320 0)
Now we can do some calculations with these, e.g. show that <math>\sin^2 x + \cos^2 x = 1</math> or define <math>\tan x = \frac{\sin x}{\cos x}</math>:
<
(newline)
Line 1,539 ⟶ 3,497:
(display (take 10 fpstan))
(newline)</
Output:
(1 0 0 0 0 0 0 0 0 0)
Line 1,550 ⟶ 3,508:
This code makes ''extensive'' use of the fact that objects can have methods and variables independent of their class. This greatly reduces the requirement for singleton classes.
<
oo::class create PowerSeries {
Line 1,686 ⟶ 3,644:
rename [cos integrate "sin"] sin
cos destroy ;# remove the dummy to make way for the real one...
rename [[PowerSeries constant 1] subtract [sin integrate]] cos</
Demonstrating:
<
sin(x) == x + -0.1667x³ + 0.008333x⁵ + -0.0001984x⁷ + 2.756e-06x⁹ + -2.505e-08x¹¹ + 1.606e-10x¹³
% cos print 7
Line 1,695 ⟶ 3,653:
1.0000035425842861
% cos evaluate [expr acos(0)]
2.473727636463901e-5</
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-rat}}
<syntaxhighlight lang="wren">import "./rat" for Rat
class Gene {
coef(n) {}
}
class Term {
construct new(gene) {
_gene = gene
_cache = []
}
gene { _gene }
[n] {
if (n < 0) return Rat.zero
if (n >= _cache.count) {
for (i in _cache.count..n) _cache.add(gene.coef(i))
}
return _cache[n]
}
}
class FormalPS {
static DISP_TERM { 12 }
static X_VAR { "x" }
construct new() {
_term = null
}
construct new(term) {
_term = term
}
term { _term }
static fromPolynomial(polynomial) {
class PolyGene is Gene {
construct new (polynomial) { _polynomial = polynomial }
coef(n) { (n < 0 || n >= _polynomial.count) ? Rat.zero : _polynomial[n] }
}
return FormalPS.new(Term.new(PolyGene.new(polynomial)))
}
copyFrom(other) { _term = other.term }
inverseCoef(n) {
var res = List.filled(n+1, null)
res[0] = _term[0].inverse
if (n > 0) {
for (i in 1..n) {
res[i] = Rat.zero
for (j in 0...i) res[i] = res[i] + _term[i - j] * res[j]
res[i] = -res[0] * res[i]
}
}
return res[n]
}
+(other) {
class AddGene is Gene {
construct new(fps, other) {
_fps = fps
_other = other
}
coef(n) { _fps.term[n] + _other.term[n] }
}
return FormalPS.new(Term.new(AddGene.new(this, other)))
}
-(other) {
class SubGene is Gene {
construct new(fps, other) {
_fps = fps
_other = other
}
coef(n) { _fps.term[n] - _other.term[n] }
}
return FormalPS.new(Term.new(SubGene.new(this, other)))
}
*(other) {
class MulGene is Gene {
construct new(fps, other) {
_fps = fps
_other = other
}
coef(n) {
var res = Rat.zero
for (i in 0..n) res = res + _fps.term[i] * _other.term[n-i]
return res
}
}
return FormalPS.new(Term.new(MulGene.new(this, other)))
}
/(other) {
class DivGene is Gene {
construct new(fps, other) {
_fps = fps
_other = other
}
coef(n) {
var res = Rat.zero
for (i in 0..n) res = res + _fps.term[i] * _other.inverseCoef(n-i)
return res
}
}
return FormalPS.new(Term.new(DivGene.new(this, other)))
}
diff() {
class DiffGene is Gene {
construct new(fps) { _fps = fps }
coef(n) { _fps.term[n+1] * Rat.new(n+1) }
}
return FormalPS.new(Term.new(DiffGene.new(this)))
}
intg() {
class IntgGene is Gene {
construct new(fps) { _fps= fps }
coef(n) { (n == 0) ? Rat.zero : _fps.term[n-1] * Rat.new(1, n) }
}
return FormalPS.new(Term.new(IntgGene.new(this)))
}
toString_(dpTerm) {
var sb = ""
var c = _term[0]
Rat.showAsInt = true
var supers = ["⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹", "¹⁰", "¹¹"]
if (c != Rat.zero) sb = sb + c.toString
for (i in 1...dpTerm) {
c = term[i]
if (c != Rat.zero) {
if (c > Rat.zero && sb.count > 0) sb = sb + " + "
var xvar = FormalPS.X_VAR
sb = sb +
((c == Rat.one) ? xvar :
(c == Rat.minusOne) ? " - %(xvar)" :
(c.num < 0) ? " - %(-c)%(xvar)" : "%(c)%(xvar)")
if (i > 1) sb = sb + "%(supers[i])"
}
}
if (sb.count == 0) sb = "0"
sb = sb + " + ..."
return sb
}
toString { toString_(FormalPS.DISP_TERM) }
}
var cos = FormalPS.new()
var sin = cos.intg()
var tan = sin/cos
cos.copyFrom(FormalPS.fromPolynomial([Rat.one]) - sin.intg())
System.print("sin(x) = %(sin)")
System.print("cos(x) = %(cos)")
System.print("tan(x) = %(tan)")
System.print("sin'(x) = %(sin.diff())")
var exp = FormalPS.new()
exp.copyFrom(FormalPS.fromPolynomial([Rat.one]) + exp.intg())
System.print("exp(x) = %(exp)")</syntaxhighlight>
{{out}}
<pre>
sin(x) = x - 1/6x³ + 1/120x⁵ - 1/5040x⁷ + 1/362880x⁹ - 1/39916800x¹¹ + ...
cos(x) = 1 - 1/2x² + 1/24x⁴ - 1/720x⁶ + 1/40320x⁸ - 1/3628800x¹⁰ + ...
tan(x) = x + 1/3x³ + 2/15x⁵ + 17/315x⁷ + 62/2835x⁹ + 1382/155925x¹¹ + ...
sin'(x) = 1 - 1/2x² + 1/24x⁴ - 1/720x⁶ + 1/40320x⁸ - 1/3628800x¹⁰ + ...
exp(x) = 1 + x + 1/2x² + 1/6x³ + 1/24x⁴ + 1/120x⁵ + 1/720x⁶ + 1/5040x⁷ + 1/40320x⁸ + 1/362880x⁹ + 1/3628800x¹⁰ + 1/39916800x¹¹ + ...
</pre>
=={{header|zkl}}==
zkl iterators (aka Walkers) are more versatile than the run-of-the-mill iterator and can be used to represent infinite sequences (eg a finite set can be padded forever or cycled over), which works well here. The Walker tweak method is used to modify iterator behavior (ie how to filter the sequence, what to do if the sequence ends, etc). The Haskell like zipWith Walker method knows how to deal with infinite sequences.
<syntaxhighlight lang="zkl">class IPS{
var [protected] w; // the coefficients of the infinite series
fcn init(w_or_a,b,c,etc){ // IPS(1,2,3) --> (1,2,3,0,0,...)
switch [arglist]{
case(Walker) { w=w_or_a.tweak(Void,0) }
else { w=vm.arglist.walker().tweak(Void,0) }
}
}
fcn __opAdd(ipf){ //IPS(1,2,3)+IPS(4,5)-->IPS(5,6,3,0,...), returns modified self
switch[arglist]{
case(1){ addConst(ipf) } // IPS + int/float
else { w=w.zipWith('+,ipf.w) } // IPS + IPS
}
self
}
fcn __opSub(ipf){ w=w.zipWith('-,ipf.w); self } // IPS - IPSHaskell
fcn __opMul(ipf){ } // stub
fcn __opDiv(x){ w.next().toFloat()/x } // *IPS/x, for integtate()
fcn __opNegate { w=w.tweak(Op("--")); self }
// integtate: b0 = 0 by convention, bn = an-1/n
fcn integrate{ w=w.zipWith('/,[1..]).push(0.0); self }
fcn diff { w=w.zipWith('*,[1..]); self }
fcn facts{ (1).walker(*).tweak(fcn(n){ (1).reduce(n,'*,1) }) } // 1!,2!...
fcn walk(n){ w.walk(n) }
fcn value(x,N=15){ ns:=[1..]; w.reduce(N,'wrap(s,an){ s + an*x.pow(ns.next()) }) }
fcn cons(k){ w.push(k); self } //--> k, a0, a1, a2, ...
// addConst(k) --> k + a0, a1, a2, ..., same as k + IPS
fcn addConst(k){ (w.next() + k) : w.push(_); self }
}</syntaxhighlight>
Add two power series. Add a const to get: 11 - (1 + 2x + 3x^2) ...
<syntaxhighlight lang="zkl">(IPS(1,2,3) + IPS(4,5)).walk(5).println();
(-IPS([1..]) + 11).walk(5).println();</syntaxhighlight>
{{out}}
<pre>
L(5,7,3,0,0)
L(10,-2,-3,-4,-5)
</pre>
Define sine in terms of a Taylor series, cos in terms of sine.
<syntaxhighlight lang="zkl">fcn sine{ // sine Taylor series: 0 + x - x^3/3! + x^5/5! - x^7/7! + x^9/9! - ...
IPS(Utils.Helpers.cycle(1.0, 0.0, -1.0, 0.0).zipWith('/,IPS.facts()))
.cons(0.0)
}
print("Sine Taylor series: "); dostuff(sine,"sin");
fcn cosine{ -sine().integrate() + 1.0 }
print("Cosine power series: "); dostuff(cosine,"cos");
fcn dostuff(ips,name){ // print series, evaluate at various points
f:='wrap(x,xnm){ v:=ips().value(x);
println("%s(%s) \U2192; %f \U394;=%f".fmt(name,xnm,v,x.Method(name)()-v));
};
ips().walk(15).println();
f(0.0,"0"); f((1.0).pi/4,"\Ubc;\U3c0;");
f((1.0).pi/2,"\Ubd;\U3c0;"); f((1.0).pi,"\U3c0;");
}</syntaxhighlight>
{{out}}
<pre>
Sine Taylor series: L(0,1,0,-0.166667,0,0.00833333,0,-0.000198413,0,2.75573e-06,0,-2.50521e-08,0,1.6059e-10,0)
sin(0) → 0.000000 Δ=0.000000
sin(¼π) → 0.707107 Δ=-0.000000
sin(½π) → 1.000000 Δ=-0.000000
sin(π) → 0.000021 Δ=-0.000021
Cosine power series: L(1,0,-0.5,0,0.0416667,0,-0.00138889,0,2.48016e-05,0,-2.75573e-07,0,2.08768e-09,0,-1.14707e-11)
cos(0) → 1.000000 Δ=0.000000
cos(¼π) → 0.707107 Δ=0.000000
cos(½π) → -0.000000 Δ=0.000000
cos(π) → -1.000004 Δ=0.000004
</pre>
|