(*
* Copyright (c) 1997-1999 Massachusetts Institute of Technology
* Copyright (c) 2003, 2007-11 Matteo Frigo
* Copyright (c) 2003, 2007-11 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
*
*)
(* 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 = 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