(* * 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