cannam@95: (*
cannam@95:  * Copyright (c) 1997-1999 Massachusetts Institute of Technology
cannam@95:  * Copyright (c) 2003, 2007-11 Matteo Frigo
cannam@95:  * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology
cannam@95:  *
cannam@95:  * This program is free software; you can redistribute it and/or modify
cannam@95:  * it under the terms of the GNU General Public License as published by
cannam@95:  * the Free Software Foundation; either version 2 of the License, or
cannam@95:  * (at your option) any later version.
cannam@95:  *
cannam@95:  * This program is distributed in the hope that it will be useful,
cannam@95:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
cannam@95:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
cannam@95:  * GNU General Public License for more details.
cannam@95:  *
cannam@95:  * You should have received a copy of the GNU General Public License
cannam@95:  * along with this program; if not, write to the Free Software
cannam@95:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
cannam@95:  *
cannam@95:  *)
cannam@95: 
cannam@95: (* The generator keeps track of numeric constants in symbolic
cannam@95:    expressions using the abstract number type, defined in this file.
cannam@95: 
cannam@95:    Our implementation of the number type uses arbitrary-precision
cannam@95:    arithmetic from the built-in Num package in order to maintain an
cannam@95:    accurate representation of constants.  This allows us to output
cannam@95:    constants with many decimal places in the generated C code,
cannam@95:    ensuring that we will take advantage of the full precision
cannam@95:    available on current and future machines.
cannam@95: 
cannam@95:    Note that we have to write our own routine to compute roots of
cannam@95:    unity, since the Num package only supplies simple arithmetic.  The
cannam@95:    arbitrary-precision operations in Num look like the normal
cannam@95:    operations except that they have an appended slash (e.g. +/ -/ */
cannam@95:    // etcetera). *)
cannam@95: 
cannam@95: open Num
cannam@95: 
cannam@95: type number = N of num
cannam@95: 
cannam@95: let makeNum n = N n
cannam@95: 
cannam@95: (* decimal digits of precision to maintain internally, and to print out: *)
cannam@95: let precision = 50
cannam@95: let print_precision = 45
cannam@95: 
cannam@95: let inveps = (Int 10) **/ (Int precision)
cannam@95: let epsilon = (Int 1) // inveps
cannam@95: 
cannam@95: let pinveps = (Int 10) **/ (Int print_precision)
cannam@95: let pepsilon = (Int 1) // pinveps
cannam@95: 
cannam@95: let round x = epsilon */ (round_num (x */ inveps))
cannam@95: 
cannam@95: let of_int n = N (Int n)
cannam@95: let zero = of_int 0
cannam@95: let one = of_int 1
cannam@95: let two = of_int 2
cannam@95: let mone = of_int (-1)
cannam@95: 
cannam@95: (* comparison predicate for real numbers *)
cannam@95: let equal (N x) (N y) = (* use both relative and absolute error *)
cannam@95:   let absdiff = abs_num (x -/ y) in
cannam@95:   absdiff <=/ pepsilon or
cannam@95:   absdiff <=/ pepsilon */ (abs_num x +/ abs_num y)
cannam@95: 
cannam@95: let is_zero = equal zero
cannam@95: let is_one = equal one
cannam@95: let is_mone = equal mone
cannam@95: let is_two = equal two
cannam@95: 
cannam@95: 
cannam@95: (* Note that, in the following computations, it is important to round
cannam@95:    to precision epsilon after each operation.  Otherwise, since the
cannam@95:    Num package uses exact rational arithmetic, the number of digits
cannam@95:    quickly blows up. *)
cannam@95: let mul (N a) (N b) = makeNum (round (a */ b)) 
cannam@95: let div (N a) (N b) = makeNum (round (a // b))
cannam@95: let add (N a) (N b) = makeNum (round (a +/ b)) 
cannam@95: let sub (N a) (N b) = makeNum (round (a -/ b))
cannam@95: 
cannam@95: let negative (N a) = (a </ (Int 0))
cannam@95: let negate (N a) = makeNum (minus_num a)
cannam@95: 
cannam@95: let greater a b = negative (sub b a)
cannam@95: 
cannam@95: let epsilonsq = epsilon */ epsilon
cannam@95: let epsilonsq2 =  (Int 100) */ epsilonsq
cannam@95: 
cannam@95: let sqr a = a */ a
cannam@95: let almost_equal (N a) (N b) = (sqr (a -/ b)) <=/ epsilonsq2
cannam@95: 
cannam@95: (* find square root by Newton's method *)
cannam@95: let sqrt a =
cannam@95:   let rec sqrt_iter guess =
cannam@95:     let newguess = div (add guess (div a guess)) two in
cannam@95:     if (almost_equal newguess guess) then newguess
cannam@95:     else sqrt_iter newguess
cannam@95:   in sqrt_iter (div a two)
cannam@95: 
cannam@95: let csub (xr, xi) (yr, yi) = (round (xr -/ yr), round (xi -/ yi))
cannam@95: let cdiv (xr, xi) r = (round (xr // r), round (xi // r))
cannam@95: let cmul (xr, xi) (yr, yi) = (round (xr */ yr -/ xi */ yi),
cannam@95:                               round (xr */ yi +/ xi */ yr))
cannam@95: let csqr (xr, xi) = (round (xr */ xr -/ xi */ xi), round ((Int 2) */ xr */ xi))
cannam@95: let cabssq (xr, xi) = xr */ xr +/ xi */ xi
cannam@95: let cconj (xr, xi) = (xr, minus_num xi)
cannam@95: let cinv x = cdiv (cconj x) (cabssq x)
cannam@95: 
cannam@95: let almost_equal_cnum (xr, xi) (yr, yi) =
cannam@95:     (cabssq (xr -/ yr,xi -/ yi)) <=/ epsilonsq2
cannam@95: 
cannam@95: (* Put a complex number to an integer power by repeated squaring: *)
cannam@95: let rec ipow_cnum x n =
cannam@95:     if (n == 0) then
cannam@95:       (Int 1, Int 0)
cannam@95:     else if (n < 0) then
cannam@95:       cinv (ipow_cnum x (- n))
cannam@95:     else if (n mod 2 == 0) then
cannam@95:       ipow_cnum (csqr x) (n / 2)
cannam@95:     else
cannam@95:       cmul x (ipow_cnum x (n - 1))
cannam@95: 
cannam@95: let twopi = 6.28318530717958647692528676655900576839433879875021164194989
cannam@95: 
cannam@95: (* Find the nth (complex) primitive root of unity by Newton's method: *)
cannam@95: let primitive_root_of_unity n =
cannam@95:     let rec root_iter guess =
cannam@95:         let newguess = csub guess (cdiv (csub guess
cannam@95:                                          (ipow_cnum guess (1 - n)))
cannam@95:                                    (Int n)) in
cannam@95:             if (almost_equal_cnum guess newguess) then newguess
cannam@95:             else root_iter newguess
cannam@95:     in let float_to_num f = (Int (truncate (f *. 1.0e9))) // (Int 1000000000)
cannam@95:     in root_iter (float_to_num (cos (twopi /. (float n))),
cannam@95: 		  float_to_num (sin (twopi /. (float n)))) 
cannam@95: 
cannam@95: let cexp n i =
cannam@95:     if ((i mod n) == 0) then
cannam@95:       (one,zero)
cannam@95:     else
cannam@95:       let (n2,i2) = Util.lowest_terms n i
cannam@95:       in let (c,s) = ipow_cnum (primitive_root_of_unity n2) i2
cannam@95:       in (makeNum c, makeNum s)
cannam@95: 
cannam@95: let to_konst (N n) =
cannam@95:   let f = float_of_num n in
cannam@95:   let f' = if f < 0.0 then f *. (-1.0) else f in
cannam@95:   let f2 = if (f' >= 1.0) then (f' -. (float (truncate f'))) else f'
cannam@95:   in let q = string_of_int (truncate(f2 *. 1.0E9))
cannam@95:   in let r = "0000000000" ^ q
cannam@95:   in let l = String.length r 
cannam@95:   in let prefix = if (f < 0.0) then "KN" else "KP" in
cannam@95:   if (f' >= 1.0) then
cannam@95:     (prefix ^ (string_of_int (truncate f')) ^ "_" ^ 
cannam@95:      (String.sub r (l - 9) 9))
cannam@95:   else
cannam@95:     (prefix ^ (String.sub r (l - 9) 9))
cannam@95: 
cannam@95: let to_string (N n) = approx_num_fix print_precision n
cannam@95: 
cannam@95: let to_float (N n) = float_of_num n
cannam@95: