cannam@95: (*
cannam@95:  * Copyright (c) 1997-1999 Massachusetts Institute of Technology
cannam@95:  * Copyright (c) 2003, 2007-11 Matteo Frigo
cannam@95:  * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology
cannam@95:  *
cannam@95:  * This program is free software; you can redistribute it and/or modify
cannam@95:  * it under the terms of the GNU General Public License as published by
cannam@95:  * the Free Software Foundation; either version 2 of the License, or
cannam@95:  * (at your option) any later version.
cannam@95:  *
cannam@95:  * This program is distributed in the hope that it will be useful,
cannam@95:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
cannam@95:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
cannam@95:  * GNU General Public License for more details.
cannam@95:  *
cannam@95:  * You should have received a copy of the GNU General Public License
cannam@95:  * along with this program; if not, write to the Free Software
cannam@95:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
cannam@95:  *
cannam@95:  *)
cannam@95: 
cannam@95: (*
cannam@95:  * the oracle decrees whether the sign of an expression should
cannam@95:  * be changed.
cannam@95:  *
cannam@95:  * Say the expression (A - B) appears somewhere.  Elsewhere in the
cannam@95:  * expression dag the expression (B - A) may appear.
cannam@95:  * The oracle determines which of the two forms is canonical.
cannam@95:  *
cannam@95:  * Algorithm: evaluate the expression at a random input, and
cannam@95:  * keep the expression with the positive sign.
cannam@95:  *)
cannam@95: 
cannam@95: let make_memoizer hash equal =
cannam@95:   let table = ref Assoctable.empty 
cannam@95:   in 
cannam@95:   (fun f k ->
cannam@95:     match Assoctable.lookup hash equal k !table with
cannam@95:       Some value -> value
cannam@95:     | None ->
cannam@95:         let value = f k in
cannam@95:         begin	
cannam@95:           table := Assoctable.insert hash k value !table;
cannam@95:           value
cannam@95:         end)
cannam@95: 
cannam@95: let almost_equal x y = 
cannam@95:   let epsilon = 1.0E-8 in
cannam@95:   (abs_float (x -. y) < epsilon) ||
cannam@95:   (abs_float (x -. y) < epsilon *. (abs_float x +. abs_float y)) 
cannam@95: 
cannam@95: let absid = make_memoizer
cannam@95:     (fun x -> Expr.hash_float (abs_float x))
cannam@95:     (fun a b -> almost_equal a b || almost_equal (-. a) b)
cannam@95:     (fun x -> x)
cannam@95: 
cannam@95: let make_random_oracle () = make_memoizer 
cannam@95:     Variable.hash 
cannam@95:     Variable.same
cannam@95:     (fun _ -> (float (Random.bits())) /. 1073741824.0)
cannam@95: 
cannam@95: let the_random_oracle = make_random_oracle ()
cannam@95: 
cannam@95: let sum_list l = List.fold_right (+.) l 0.0
cannam@95: 
cannam@95: let eval_aux random_oracle =
cannam@95:   let memoizing = make_memoizer Expr.hash (==) in
cannam@95:   let rec eval x = 
cannam@95:     memoizing
cannam@95:       (function
cannam@95: 	| Expr.Num x -> Number.to_float x
cannam@95: 	| Expr.NaN x -> Expr.transcendent_to_float x
cannam@95: 	| Expr.Load v -> random_oracle v
cannam@95: 	| Expr.Store (v, x) -> eval x
cannam@95: 	| Expr.Plus l -> sum_list (List.map eval l)
cannam@95: 	| Expr.Times (a, b) -> (eval a) *. (eval b)
cannam@95: 	| Expr.CTimes (a, b) -> 
cannam@95: 	    1.098612288668109691395245236 +. 
cannam@95: 	       1.609437912434100374600759333 *. (eval a) *. (eval b)
cannam@95: 	| Expr.CTimesJ (a, b) -> 
cannam@95: 	    0.9102392266268373936142401657 +. 
cannam@95: 	      0.6213349345596118107071993881 *. (eval a) *. (eval b)
cannam@95: 	| Expr.Uminus x -> -. (eval x))
cannam@95:       x
cannam@95:   in eval
cannam@95: 
cannam@95: let eval = eval_aux the_random_oracle
cannam@95: 
cannam@95: let should_flip_sign node = 
cannam@95:   let v = eval node in
cannam@95:   let v' = absid v in
cannam@95:   not (almost_equal v v')
cannam@95: 
cannam@95: (*
cannam@95:  * determine with high probability if two expressions are equal.
cannam@95:  *
cannam@95:  * The test is randomized: if the two expressions have the
cannam@95:  * same value for NTESTS random inputs, then they are proclaimed
cannam@95:  * equal.  (Note that two distinct linear functions L1(x0, x1, ..., xn)
cannam@95:  * and L2(x0, x1, ..., xn) have the same value with probability
cannam@95:  * 0 for random x's, and thus this test is way more paranoid than
cannam@95:  * necessary.)
cannam@95:  *)
cannam@95: let likely_equal a b =
cannam@95:   let tolerance = 1.0e-8
cannam@95:   and ntests = 20
cannam@95:   in
cannam@95:   let rec loop n =
cannam@95:     if n = 0 then 
cannam@95:       true
cannam@95:     else
cannam@95:       let r = make_random_oracle () in
cannam@95:       let va = eval_aux r a
cannam@95:       and vb = eval_aux r b
cannam@95:       in
cannam@95:       if (abs_float (va -. vb)) > 
cannam@95: 	   tolerance *. (abs_float va +. abs_float vb +. 0.0001)
cannam@95:       then
cannam@95: 	false
cannam@95:       else
cannam@95: 	loop (n - 1)
cannam@95:   in
cannam@95:   match (a, b) with
cannam@95: 
cannam@95:     (* 
cannam@95:      * Because of the way eval is constructed, we have
cannam@95:      *     eval (Store (v, x)) == eval x
cannam@95:      * However, we never consider the two expressions equal
cannam@95:      *)
cannam@95:   | (Expr.Store _, _) -> false
cannam@95:   | (_, Expr.Store _) -> false
cannam@95: 
cannam@95:     (*
cannam@95:      * Expressions of the form ``Uminus (Store _)''
cannam@95:      * are artifacts of algsimp
cannam@95:      *)
cannam@95:   | ((Expr.Uminus (Expr.Store _)), _) -> false
cannam@95:   | (_, Expr.Uminus (Expr.Store _)) -> false
cannam@95: 
cannam@95:   | _ -> loop ntests
cannam@95: 
cannam@95: let hash x =
cannam@95:   let f = eval x in
cannam@95:   truncate (f *. 65536.0)