Formal power series: Difference between revisions

Line 2,241:

=={{header|Phix}}==

<lang Phix>~~enum type fps_type FPS_UNDEF = 0,~~

with javascript_semantics

FPS_CONST,

enum FPS_UNDEF = 0,

FPS_ADD,

FPS_CONST,

FPS_SUB,

FPS_ADD,

FPS_MUL,

FPS_SUB,

FPS_DIV,

FPS_MUL,

FPS_DERIV,

FPS_DIV,

FPS_INT

FPS_DERIV,

end type

FPS_INT,

FPS_MAX=$

enum FPS_TYPE, FPS_S1, FPS_S2, FPS_A0

sequence fpss = {}

type fps(object id)

return integer(id) and id>=1 and id<=length(fpss)

end type

type fpsn(object id)

return id=NULL or fps(id)

end type

function fps_new(fps_type ft=FPS_UNDEF, fpsn s1=0, s2=0, atom a0=0)

fpss = append(fpss,{ft,s1,s2,a0})

fps fpsid = length(fpss)

return fpsid

end function

-- as per C, for (eg) self or mutually recursive definitions.

procedure fps_redefine(fps fpsid, fps_type ft, fpsn s1id, s2id, object a0="")

fpss[fpsid][FPS_TYPE] = ft

fpss[fpsid][FPS_S1] = s1id

fpss[fpsid][FPS_S2] = s2id

if atom(a0) then

fpss[fpsid][FPS_A0] = a0

end if

end procedure

function fps_const(atom a0)

fps x = fps_new(FPS_CONST,a0:=a0)

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

return x

end function

type fps_type(integer f)

constant INF = 1e300*1e300,

return f>=FPS_UNDEF and f<=FPS_MAX

NAN = -(INF/INF)

end type

/* Taking the n-th term of series. This is where actual work is done. */

enum FPS_TYPE, FPS_S1, FPS_S2, FPS_A0

function term(fps x, int n)

sequence fpss = {}

atom ret = 0

type fps(object id)

{fps_type ft, fpsn s1id, fpsn s2id, atom a0} = fpss[x]

return integer(id) and id>=1 and id<=length(fpss)

-- FPS_TYPE, FPS_S1, FPS_S2, FPS_A0 <-- nb above must match

end type

switch ft do

case FPS_CONST: ret := iff(n>0 ? 0 : a0)

type fpsn(object id)

case FPS_ADD: ret := term(s1id, n) + term(s2id, n)

return id=NULL or fps(id)

case FPS_SUB: ret := term(s1id, n) - term(s2id, n)

end type

case FPS_MUL:

for i=0 to n do

function fps_new(fps_type ft=FPS_UNDEF, fpsn s1=0, s2=0, atom a0=0)

ret += term(s1id, i) * term(s2id, n-i)

fpss = append(fpss,{ft,s1,s2,a0})

end for

fps fpsid = length(fpss)

case FPS_DIV:

return fpsid

if not term(s2id, 0) then return NAN end if

end function

ret = term(s1id, n)

for i=1 to n do

-- as per C, for (eg) self or mutually recursive definitions.

ret -= term(s2id, i) * term(x, n-i) / term(s2id, 0)

procedure fps_redefine(fps fpsid, fps_type ft, fpsn s1id, s2id, object a0="")

end for

fpss[fpsid][FPS_TYPE] = ft

case FPS_DERIV: ret := n * term(s1id, n+1)

fpss[fpsid][FPS_S1] = s1id

case FPS_INT: ret := iff(n=0 ? a0 : term(s1id, n-1)/n)

fpss[fpsid][FPS_S2] = s2id

default: ret := 9/0 -- (fatal error)

if atom(a0) then

end switch

fpss[fpsid][FPS_A0] = a0

return ret

end if

end function

end procedure

procedure term9(string txt, fps x)

function fps_const(atom a0)

printf(1,"%s:",{txt})

fps x = fps_new(FPS_CONST,a0:=a0)

for i=0 to 9 do printf(1," %g", term(x, i)) end for

-- (aside: in the above, the ":=a0" refers to the local namespace

printf(1,"\n")

-- as usual, whereas "a0:=" refers to the param namespace

end procedure

-- /inside/ the () of fps_new(), so there is no conflict.)

return x

procedure main()

end function

fps one = fps_const(1)

fps fcos = fps_new() /* cosine */

constant INF = 1e300*1e300,

fps fsin = fps_new(FPS_INT,fcos) /* sine */

NAN = -(INF/INF)

fps ftan = fps_new(FPS_DIV,fsin,fcos) /* tangent */

/* Taking the n-th term of series. This is where actual work is done. */

/* redefine cos to complete the mutual recursion */

function term(fps x, int n)

fps_redefine(fcos, FPS_SUB, one, fps_new(FPS_INT,fsin))

atom ret = 0

fps fexp = fps_const(1); /* exponential */

{fps_type ft, fpsn s1id, fpsn s2id, atom a0} = fpss[x]

/* make exp recurse on self */

-- FPS_TYPE, FPS_S1, FPS_S2, FPS_A0 <-- nb above must match

fps_redefine(fexp, FPS_INT, fexp, 0);

switch ft do

case FPS_CONST: ret := iff(n>0 ? 0 : a0)

term9("Sin",fsin)

case FPS_ADD: ret := term(s1id, n) + term(s2id, n)

term9("Cos",fcos)

case FPS_SUB: ret := term(s1id, n) - term(s2id, n)

term9("Tan",ftan)

case FPS_MUL:

term9("Exp",fexp)

for i=0 to n do

end procedure

ret += term(s1id, i) * term(s2id, n-i)

main()</lang>

end for

case FPS_DIV:

if not term(s2id, 0) then return NAN end if

ret = term(s1id, n)

for i=1 to n do

ret -= term(s2id, i) * term(x, n-i) / term(s2id, 0)

end for

case FPS_DERIV: ret := n * term(s1id, n+1)

case FPS_INT: ret := iff(n=0 ? a0 : term(s1id, n-1)/n)

default: ret := 9/0 -- (fatal error)

end switch

return ret

end function

procedure term9(string txt, fps x)

printf(1,"%s:",{txt})

for i=0 to 9 do printf(1," %g", term(x, i)) end for

printf(1,"\n")

end procedure

procedure main()

fps one = fps_const(1)

fps fcos = fps_new() /* cosine */

fps fsin = fps_new(FPS_INT,fcos) /* sine */

fps ftan = fps_new(FPS_DIV,fsin,fcos) /* tangent */

/* redefine cos to complete the mutual recursion */

fps_redefine(fcos, FPS_SUB, one, fps_new(FPS_INT,fsin))

fps fexp = fps_const(1); /* exponential */

/* make exp recurse on self */

fps_redefine(fexp, FPS_INT, fexp, 0);

term9("Sin",fsin)

term9("Cos",fcos)

term9("Tan",ftan)

term9("Exp",fexp)

end procedure

main()

<pre>

Sin: 0 1 0 -0.166667 0 0.~~00833333~~ 0 -0.~~000198413~~ 0 2.75573e-6

Sin: 0 1 0 -0.166667 0 0.008333 0 -0.000198 0 2.75573e-6

Cos: 1 0 -0.5 0 0.~~0416667~~ 0 -0.~~00138889~~ 0 2.48016e-5 0

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.~~0539683~~ 0 0.~~0218695~~

Tan: 0 1 0 0.333333 0 0.133333 0 0.053968 0 0.021869

Exp: 1 1 0.5 0.166667 0.~~0416667~~ 0.~~00833333~~ 0.~~00138889~~ 0.~~000198413~~ 2.48016e-5 2.75573e-6

Exp: 1 1 0.5 0.166667 0.041667 0.008333 0.001389 0.000198 2.48016e-5 2.75573e-6

</pre>