furnace/extern/fftw/genfft/number.ml

165 lines
5.6 KiB
OCaml
Raw Normal View History

(*
* Copyright (c) 1997-1999 Massachusetts Institute of Technology
* Copyright (c) 2003, 2007-14 Matteo Frigo
* Copyright (c) 2003, 2007-14 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*
*)
(* 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 ||
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 epsilonsq = epsilon */ epsilon
let epsilonsq2 = (Int 100) */ epsilonsq
let sqr a = a */ a
let almost_equal (N a) (N b) = (sqr (a -/ b)) <=/ epsilonsq2
(* find square root by Newton's method *)
let sqrt a =
let rec sqrt_iter guess =
let newguess = div (add guess (div a guess)) two in
if (almost_equal newguess guess) then newguess
else sqrt_iter newguess
in sqrt_iter (div a two)
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 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 to_konst (N n) =
let f = float_of_num n in
let f' = if f < 0.0 then f *. (-1.0) else f 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 let prefix = if (f < 0.0) then "KN" else "KP" in
if (f' >= 1.0) then
(prefix ^ (string_of_int (truncate f')) ^ "_" ^
(String.sub r (l - 9) 9))
else
(prefix ^ (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