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