Formal power series: Difference between revisions

=={{header|Ada}}==

The Taylor series package is generic to be instantiated with any rational type implementation provided by the task [[Rational Arithmetic]]:

 

generic

Zero : constant Taylor_Series := (0 => Rational_Numbers.Zero);

One : constant Taylor_Series := (0 => Rational_Numbers.One);

The package implementation:

function Normalize (A : Taylor_Series) return Taylor_Series is

begin

end Value;

 

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;

Put ("Cos ="); Put (Cos (5)); Put_Line (" ...");

Put ("Sin ="); Put (Integral (Cos (5))); Put_Line (" ...");

Sample output:

<pre>

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|Clojure}}==

 

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))

(take-while notz? (map z+ (+zs ps0) (+zs 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]]

(let [mrest1 (if (or (nil? rest1) (zero? p0)) nil, (map #(* p0 %) rest1))

accum (cond (nil? resta) mrest1, (nil? mrest1) resta, :else (ps+ resta mrest1))]

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]

 

(defn integrate [ps]

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]))

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))

 

 

(println (take 20 (ps+ (ps* sin sin) (ps* cos cos))))

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))

 

(println (take 10 sinx))

 

=={{header|Common Lisp}}==

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")

#:+ #:- #:* #:/))

 

 

===Lazy Primitives===

 

thunk value)

 

(t (let ((val (funcall (promise-thunk object))))

(setf (promise-thunk object) nil

 

===Lazy Lists===

 

lar

ldr)

 

(defmacro lons (lar ldr)

 

A few utilities to make working with lazy lists easier.

 

(funcall function (lar lons) (ldr lons)))

 

 

(defun up-from (n)

 

===Formal Power Series===

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)

 

 

(defun diff (f)

 

===Example===

 

(locally (declare (special *cosx*))

(delay (int (force *cosx*)))))

 

(defparameter *cosx*

 

<pre>FPS > (force-list (take 10 (coeffs (force *sinx*))))

(* 3 (force *sinx*))))))

(3 -3 -1/2 1/2 1/24 -1/40 -1/720 1/1680 1/40320 -1/120960)</pre>

 

=={{header|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.

(require 'math)

;; converts a finite polynomial (a_0 a_1 .. a_n) to an infinite serie (a_0 ..a_n 0 0 0 ...)

 

 

{{out}}

(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)

→ 2

 

 

=={{header|Go}}==

 

import (

fmt.Println("sin:", partialSeries(sin))

fmt.Println("cos:", partialSeries(cos))

Output:

<pre>

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

 

cosx = 1 - int sinx

 

 

Output:

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.@#

mult=: +//.@(*/)

plus=: +/@,:

 

See also section 2 at http://www.jsoftware.com/papers/tot1.htm

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

represents the power series 1 + x/2 + x/3 + ... because 1/(1+.) maps i to 1/(i+1).

 

 

The exponential power series, Σ (x^i)/i!, can be represented in jq by the filter:

 

assuming "factorial" is defined in the usual way:

reduce range(1; . + 1) as $i

 

For ease of reference, we shall also define ps_exp as 1/factorial:

 

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:

should evaluate to the number e approximately; using the version of ps_evaluate defined below, we find:

evaluates to 2.7182818284590455.

 

The following function definitions are useful for other power series:

. as $x | n as $n

For example, the power series 1 + Σ ( (x/i)^n ) where the summation is over i>0 can be written:

 

The power series representation of ln(1 + x) is as follows:

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;

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:

 

For example, poly( [1,2,3] ) represents the finite power series: 1 + 2x + 3x^2.

jq's "+" operator can be used to add two power series with the intended semantics;

for example:

 

is equal to poly([]), i.e. 0.

 

====Multiplication, Differentiation and Integration====

def M(s;t):

. as $i | reduce range(0; 1+$i) as $k

. as $i

| if $i == 0 then 0 else (($i-1)|s) /$i end;

====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:

def abs: if . < 0 then -. else . end;

reduce range(0;k) as $i

if . then ((($i|s) - ($i|t))|abs) <= eps

else .

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.

def ps_eval(p; k):

. as $x

. as $x

| [0, 1, 0, 1, 1]

====Examples====

 

def abs: if . < 0 then -. else . end;

 

def pi: 4 * (1|atan);

''''cos == I(sin)''''

# Verify that the first 100 terms of I(cos) and of sin are the same:

 

 

((pi | ps_evaluate(I(ps_cos))) - (pi | ps_evaluate(ps_sin))) | abs < 1e-15

''''cos == 1 - I(sin)''''

 

ps_equal( ps_cos; ps_at(0) - I(ps_sin); 100; 1e-5)

 

((pi | ps_evaluate(ps_cos)) - (pi | ps_evaluate(ps_at(0) - I(ps_sin)))) | abs < 1e-15

 

 

 

 

 

 

 

=={{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,

 

tangent = sine / cosine

 

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.

sin = Series[Sin[x], {x, 0, 8}];

<pre>O[x]^9</pre>

 

=={{header|PARI/GP}}==

Uses the built-in power series handling. Change default(seriesprecision) to get more terms.

 

=={{header|Perl}}==

Creating a new arithmetic type in perl is relatively easy, using the "overload" module which comes with perl.

 

 

package FPS;

use strict;

print "exp(x) ~= $exp\n";

 

 

1;

__END__

{{out}}

<pre>

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

</pre>

 

For a version which *does* use proper lazy lists, see [[Formal power series/Perl]]

 

{{out}}

 

=={{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))

(caadr Args)

we can build a formal power series functionality:

 

(de fpsOne (N)

(if (=0 N)

1.0

Test:

(for N (range 1 11 2)

(prin " " (round (fpsSin N) 9)) )

(for N (range 0 6)

(prin " " (round (fpsExp N) 7)) )

Output:

<pre>SIN: 1.000000000 -0.166666667 0.008333333 -0.000198413 0.000002756 -0.000000025

=={{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

 

assert s == integc, "The integral of cos should be sin"

 

'''Sample output'''

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:

print(list(islice(sinx, 10)))

print("sine")

 

'''Sample output'''

===Operator overloading===

Define an iterator class as a polynomial to provide overloaded operators and automatic tee-ing. It is kind of overkill.

from fractions import Fraction

 

 

for n,x in zip(("sin", "cos", "tan", "exp"), (sinx, cosx, tanx, expx)):

 

=={{header|Racket}}==

Using '''Lazy Racket''':

 

#lang lazy

 

;; series of natural numbers greater then n

(define (enum n) (cons n (enum (+ 1 n ))))

 

Examples:

 

(define <sin> (cons 0 (int <cos>)))

(define <cos> (cons 1 (scale -1 (int <sin>))))

-> (!! (take 10 (</> <sin> <cos>)))

'(0 1 0 1/3 0 2/15 0 17/315 0 62/2835)

 

=={{header|Scheme}}==

Definitions of operations on lazy lists:

(syntax-rules ()

((_ lar ldr) (delay (cons lar (delay ldr))))))

 

(define (repeat n)

Definitions of operations on formal power series:

(apply lap + llists))

 

 

(define (fpsdif 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)))))

 

 

(display (take 10 fpscos))

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)

 

 

(display (take 10 fpstan))

Output:

(1 0 0 0 0 0 0 0 0 0)

 

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 {

rename [cos integrate "sin"] sin

cos destroy ;# remove the dummy to make way for the real one...

Demonstrating:

sin(x) == x + -0.1667x³ + 0.008333x⁵ + -0.0001984x⁷ + 2.756e-06x⁹ + -2.505e-08x¹¹ + 1.606e-10x¹³

% cos print 7

1.0000035425842861

% cos evaluate [expr acos(0)]

 

=={{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.

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,...)

// addConst(k) --> k + a0, a1, a2, ..., same as k + IPS

fcn addConst(k){ (w.next() + k) : w.push(_); self }

Add two power series. Add a const to get: 11 - (1 + 2x + 3x^2) ...

{{out}}

<pre>

</pre>

Define sine in terms of a Taylor series, cos in terms of sine.

IPS(Utils.Helpers.cycle(1.0, 0.0, -1.0, 0.0).zipWith('/,IPS.facts()))

.cons(0.0)

f(0.0,"0"); f((1.0).pi/4,"\Ubc;\U3c0;");

f((1.0).pi/2,"\Ubd;\U3c0;"); f((1.0).pi,"\U3c0;");

{{out}}

<pre>