number.ml

FFTW, a collection of fast C routines to compute the Discrete Fourier Transform in one or more dime

(*
 * Copyright (c) 1997-1999, 2003 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 *)

(* $Id: number.ml,v 1.24 2003/03/16 23:43:46 stevenj Exp $ *)

(* The generator keeps track of numeric constants in symbolic
   expressions using the abstract number type, defined in this file.

   Our implementation of the number type uses arbitrary-precision
   arithmetic from the built-in Num package in order to maintain an
   accurate representation of constants.  This allows us to output
   constants with many decimal places in the generated C code,
   ensuring that we will take advantage of the full precision
   available on current and future machines.

   Note that we have to write our own routine to compute roots of
   unity, since the Num package only supplies simple arithmetic.  The
   arbitrary-precision operations in Num look like the normal
   operations except that they have an appended slash (e.g. +/ -/ */
   // etcetera). *)

open Num

type number = N of num

let makeNum n = N n

(* decimal digits of precision to maintain internally, and to print out: *)
let precision = 50
let print_precision = 45

let inveps = (Int 10) **/ (Int precision)
let epsilon = (Int 1) // inveps

let pinveps = (Int 10) **/ (Int print_precision)
let pepsilon = (Int 1) // pinveps

let round x = epsilon */ (round_num (x */ inveps))

let of_int n = N (Int n)
let zero = of_int 0
let one = of_int 1
let two = of_int 2
let mone = of_int (-1)

(* comparison predicate for real numbers *)
let equal (N x) (N y) = (* use both relative and absolute error *)
  let absdiff = abs_num (x -/ y) in
  absdiff <=/ pepsilon or
  absdiff <=/ pepsilon */ (abs_num x +/ abs_num y)

let is_zero = equal zero
let is_one = equal one
let is_mone = equal mone
let is_two = equal two


(* Note that, in the following computations, it is important to round
   to precision epsilon after each operation.  Otherwise, since the
   Num package uses exact rational arithmetic, the number of digits
   quickly blows up. *)
let mul (N a) (N b) = makeNum (round (a */ b)) 
let div (N a) (N b) = makeNum (round (a // b))
let add (N a) (N b) = makeNum (round (a +/ b)) 
let sub (N a) (N b) = makeNum (round (a -/ b))

let negative (N a) = (a </ (Int 0))
let negate (N a) = makeNum (minus_num a)

let greater a b = negative (sub b a)

let csub (xr, xi) (yr, yi) = (round (xr -/ yr), round (xi -/ yi))
let cdiv (xr, xi) r = (round (xr // r), round (xi // r))
let cmul (xr, xi) (yr, yi) = (round (xr */ yr -/ xi */ yi),
                              round (xr */ yi +/ xi */ yr))
let csqr (xr, xi) = (round (xr */ xr -/ xi */ xi), round ((Int 2) */ xr */ xi))
let cabssq (xr, xi) = xr */ xr +/ xi */ xi
let cconj (xr, xi) = (xr, minus_num xi)
let cinv x = cdiv (cconj x) (cabssq x)

let epsilonsq = epsilon */ epsilon
let epsilonsq2 =  (Int 100) */ epsilonsq

let almost_equal_cnum (xr, xi) (yr, yi) =
    (cabssq (xr -/ yr,xi -/ yi)) <=/ epsilonsq2

(* Put a complex number to an integer power by repeated squaring: *)
let rec ipow_cnum x n =
    if (n == 0) then
      (Int 1, Int 0)
    else if (n < 0) then
      cinv (ipow_cnum x (- n))
    else if (n mod 2 == 0) then
      ipow_cnum (csqr x) (n / 2)
    else
      cmul x (ipow_cnum x (n - 1))

let twopi = 6.28318530717958647692528676655900576839433879875021164194989

(* Find the nth (complex) primitive root of unity by Newton's method: *)
let primitive_root_of_unity n =
    let rec root_iter guess =
        let newguess = csub guess (cdiv (csub guess
                                         (ipow_cnum guess (1 - n)))
                                   (Int n)) in
            if (almost_equal_cnum guess newguess) then newguess
            else root_iter newguess
    in let float_to_num f = (Int (truncate (f *. 1.0e9))) // (Int 1000000000)
    in root_iter (float_to_num (cos (twopi /. (float n))),
		  float_to_num (sin (twopi /. (float n)))) 

let cexp n i =
    if ((i mod n) == 0) then
      (one,zero)
    else
      let (n2,i2) = Util.lowest_terms n i
      in let (c,s) = ipow_cnum (primitive_root_of_unity n2) i2
      in (makeNum c, makeNum s)

let unparse (N n) =
  let f = float_of_num n in
  let f2 = if (f >= 1.0) then (f -. (float (truncate f))) else f
  in let q = string_of_int (truncate(f2 *. 1.0E9))
  in let r = "0000000000" ^ q
  in let l = String.length r
  in
  if (f >= 1.0) then
    ("K" ^ (string_of_int (truncate f)) ^ "_" ^ 
     (String.sub r (l - 9) 9))
  else
    ("K" ^ (String.sub r (l - 9) 9))

let to_string (N n) = approx_num_fix print_precision n

let to_float (N n) = float_of_num n