Power serious

Power serious

[Jon Harrop]

I stumbled upon the following article that describes a remarkably simple 
implementation of arithmetic over power series in Haskell:

  http://www.cs.dartmouth.edu/~doug/powser.html

This is the only compelling example of Haskell I have ever seen and I'd like 
to see this rewritten in other languages for comparison. Has anyone tried to 
translate this into OCaml?

-- 
Dr Jon Harrop, Flying Frog Consultancy Ltd.
http://www.ffconsultancy.com/?e

Re: [Caml-list] Power serious

[David McClain]

I believe that you'll find this stuff in SICP too... I know I have plenty of Lisp and Scheme examples of this sort of thing... e.g., computing tan from the ratio of the sin and cos series, done in lazy streams along with series convergence acceleration, a la John Houston's paper...

- DM

On Feb 7, 2010, at 20:20 PM, Jon Harrop wrote:

> 
> I stumbled upon the following article that describes a remarkably simple 
> implementation of arithmetic over power series in Haskell:
> 
>  http://www.cs.dartmouth.edu/~doug/powser.html
> 
> This is the only compelling example of Haskell I have ever seen and I'd like 
> to see this rewritten in other languages for comparison. Has anyone tried to 
> translate this into OCaml?
> 
> -- 
> Dr Jon Harrop, Flying Frog Consultancy Ltd.
> http://www.ffconsultancy.com/?e
> 
> _______________________________________________
> Caml-list mailing list. Subscription management:
> http://yquem.inria.fr/cgi-bin/mailman/listinfo/caml-list
> Archives: http://caml.inria.fr
> Beginner's list: http://groups.yahoo.com/group/ocaml_beginners
> Bug reports: http://caml.inria.fr/bin/caml-bugs
> 

Dr. David McClain
dbm@refined-audiometrics.com

RE: [Caml-list] Power serious

[Harrison, John R]

[-- Attachment #1: Type: text/plain, Size: 1427 bytes --]

Hi Jon,

| I stumbled upon the following article that describes a remarkably simple
| implementation of arithmetic over power series in Haskell:
|
|   http://www.cs.dartmouth.edu/~doug/powser.html
|
| This is the only compelling example of Haskell I have ever seen and I'd like
| to see this rewritten in other languages for comparison. Has anyone tried to
| translate this into OCaml?

I once hacked up an OCaml implementation of the algorithms in a similar
paper by Doug McIlroy, "Squinting at Power Series". See the attached
code, which redefines the usual arithmetic operators like "+" over power
series (note that ** is composition, not exponentiation). For example:

  # let atan = integrate(const 1 / (const 1 + x * x));;
  val atan : int -> Num.num = <fun>
  # let tan = revert atan;;
  val tan : int -> Num.num = <fun>
  # List.map (fun i -> string_of_num(tan i)) [0;1;2;3;4;5;6;7;8];;
  - : string list = ["0"; "1"; "0"; "1/3"; "0"; "2/15"; "0"; "17/315"; "0"]

I wanted this for some real power series computations on Bessel function
asymptotic expansions (http://www.cl.cam.ac.uk/~jrh13/papers/bessel.html).
To do anything serious, I found, not surprisingly, that I needed to stick
memoization in various places. Even the trivial example above takes an
appreciable time without it. Still, it doesn't uglify the code much to
add a simple cacheing function all over the place.

John.


[-- Attachment #2: powser.ml --]
[-- Type: application/octet-stream, Size: 2088 bytes --]

(* ========================================================================= *)
(* Power series expansions, following McIlroy's "Squinting at Power Series". *)
(* ========================================================================= *)

open Num;;

let num_0 = num_of_int 0
and num_1 = num_of_int 1;;

(* ------------------------------------------------------------------------- *)
(* Auxiliary functions.                                                      *)
(* ------------------------------------------------------------------------- *)

let sum (m,n) f =
  let rec sumf i acc =
    if n < i then acc else sumf (i + 1) (f(i) +/ acc) in
  sumf m num_0;;

let rest a n = a(n + 1);;

(* ------------------------------------------------------------------------- *)
(* Basic power series.                                                       *)
(* ------------------------------------------------------------------------- *)

let x n = if n = 1 then num_1 else num_0;;

let constant k = fun n -> if n = 0 then k else num_0;;

let const n = constant(num_of_int n);;

let xx a n = if n = 0 then num_0 else a (n - 1);;

(* ------------------------------------------------------------------------- *)
(* Constructors for power series.                                            *)
(* ------------------------------------------------------------------------- *)

let differentiate f n = num_of_int(n + 1) */ f(n + 1);;

let integrate f n = if n = 0 then num_0 else f(n - 1) // num_of_int n;;

let rec ( * ) a b n = a(0) */ b(n) +/ xx(rest a * b) n;;

let (+) a b n = a(n) +/ b(n);;

let (-) a b n = a(n) -/ b(n);;

let rec invert a n =
  (constant(num_1 // a 0) * (const 1 - xx(rest a * invert a))) n;;

let rec (/) a b =
  if a(0) =/ num_0 & b(0) =/ num_0 then rest a / rest b else a * invert(b);;

let ( ** ) a b n =
  let rec compfn i =
    if i = n then constant (a i) else (constant(a i)) + (b * rest compfn i) in
  compfn 0 n;;

let revert f =
  let f1' = constant(num_1 // f(1)) and f'' = rest(rest f) in
  let rec p n = (f1' * (const 1 - xx(p * p * (f'' ** xx p)))) n in
  xx p;;