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