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: 
Chris@42: open Util
Chris@42: open Expr
Chris@42: 
Chris@42: let node_insert x =  Assoctable.insert Expr.hash x
Chris@42: let node_lookup x =  Assoctable.lookup Expr.hash (==) x
Chris@42: 
Chris@42: (*************************************************************
Chris@42:  * Algebraic simplifier/elimination of common subexpressions
Chris@42:  *************************************************************)
Chris@42: module AlgSimp : sig 
Chris@42:   val algsimp : expr list -> expr list
Chris@42: end = struct
Chris@42: 
Chris@42:   open Monads.StateMonad
Chris@42:   open Monads.MemoMonad
Chris@42:   open Assoctable
Chris@42: 
Chris@42:   let fetchSimp = 
Chris@42:     fetchState >>= fun (s, _) -> returnM s
Chris@42:   let storeSimp s =
Chris@42:     fetchState >>= (fun (_, c) -> storeState (s, c))
Chris@42:   let lookupSimpM key =
Chris@42:     fetchSimp >>= fun table ->
Chris@42:       returnM (node_lookup key table)
Chris@42:   let insertSimpM key value =
Chris@42:     fetchSimp >>= fun table ->
Chris@42:       storeSimp (node_insert key value table)
Chris@42: 
Chris@42:   let subset a b =
Chris@42:     List.for_all (fun x -> List.exists (fun y -> x == y) b) a
Chris@42: 
Chris@42:   let structurallyEqualCSE a b = 
Chris@42:     match (a, b) with
Chris@42:     | (Num a, Num b) -> Number.equal a b
Chris@42:     | (NaN a, NaN b) -> a == b
Chris@42:     | (Load a, Load b) -> Variable.same a b
Chris@42:     | (Times (a, a'), Times (b, b')) ->
Chris@42:  	((a == b) && (a' == b')) ||
Chris@42:  	((a == b') && (a' == b))
Chris@42:     | (CTimes (a, a'), CTimes (b, b')) ->
Chris@42:  	((a == b) && (a' == b')) ||
Chris@42:  	((a == b') && (a' == b))
Chris@42:     | (CTimesJ (a, a'), CTimesJ (b, b')) -> ((a == b) && (a' == b'))
Chris@42:     | (Plus a, Plus b) -> subset a b && subset b a
Chris@42:     | (Uminus a, Uminus b) -> (a == b)
Chris@42:     | _ -> false
Chris@42: 
Chris@42:   let hashCSE x = 
Chris@42:     if (!Magic.randomized_cse) then
Chris@42:       Oracle.hash x
Chris@42:     else
Chris@42:       Expr.hash x
Chris@42: 
Chris@42:   let equalCSE a b = 
Chris@42:     if (!Magic.randomized_cse) then
Chris@42:       (structurallyEqualCSE a b || Oracle.likely_equal a b)
Chris@42:     else
Chris@42:       structurallyEqualCSE a b
Chris@42: 
Chris@42:   let fetchCSE = 
Chris@42:     fetchState >>= fun (_, c) -> returnM c
Chris@42:   let storeCSE c =
Chris@42:     fetchState >>= (fun (s, _) -> storeState (s, c))
Chris@42:   let lookupCSEM key =
Chris@42:     fetchCSE >>= fun table ->
Chris@42:       returnM (Assoctable.lookup hashCSE equalCSE key table)
Chris@42:   let insertCSEM key value =
Chris@42:     fetchCSE >>= fun table ->
Chris@42:       storeCSE (Assoctable.insert hashCSE key value table)
Chris@42: 
Chris@42:   (* memoize both x and Uminus x (unless x is already negated) *) 
Chris@42:   let identityM x =
Chris@42:     let memo x = memoizing lookupCSEM insertCSEM returnM x in
Chris@42:     match x with
Chris@42: 	Uminus _ -> memo x 
Chris@42:       |	_ -> memo x >>= fun x' -> memo (Uminus x') >> returnM x'
Chris@42: 
Chris@42:   let makeNode = identityM
Chris@42: 
Chris@42:   (* simplifiers for various kinds of nodes *)
Chris@42:   let rec snumM = function
Chris@42:       n when Number.is_zero n -> 
Chris@42: 	makeNode (Num (Number.zero))
Chris@42:     | n when Number.negative n -> 
Chris@42: 	makeNode (Num (Number.negate n)) >>= suminusM
Chris@42:     | n -> makeNode (Num n)
Chris@42: 
Chris@42:   and suminusM = function
Chris@42:       Uminus x -> makeNode x
Chris@42:     | Num a when (Number.is_zero a) -> snumM Number.zero
Chris@42:     | a -> makeNode (Uminus a)
Chris@42: 
Chris@42:   and stimesM = function 
Chris@42:     | (Uminus a, b) -> stimesM (a, b) >>= suminusM
Chris@42:     | (a, Uminus b) -> stimesM (a, b) >>= suminusM
Chris@42:     | (NaN I, CTimes (a, b)) -> stimesM (NaN I, b) >>= 
Chris@42: 	fun ib -> sctimesM (a, ib)
Chris@42:     | (NaN I, CTimesJ (a, b)) -> stimesM (NaN I, b) >>= 
Chris@42: 	fun ib -> sctimesjM (a, ib)
Chris@42:     | (Num a, Num b) -> snumM (Number.mul a b)
Chris@42:     | (Num a, Times (Num b, c)) -> 
Chris@42: 	snumM (Number.mul a b) >>= fun x -> stimesM (x, c)
Chris@42:     | (Num a, b) when Number.is_zero a -> snumM Number.zero
Chris@42:     | (Num a, b) when Number.is_one a -> makeNode b
Chris@42:     | (Num a, b) when Number.is_mone a -> suminusM b
Chris@42:     | (a, b) when is_known_constant b && not (is_known_constant a) -> 
Chris@42: 	stimesM (b, a)
Chris@42:     | (a, b) -> makeNode (Times (a, b))
Chris@42: 
Chris@42:   and sctimesM = function 
Chris@42:     | (Uminus a, b) -> sctimesM (a, b) >>= suminusM
Chris@42:     | (a, Uminus b) -> sctimesM (a, b) >>= suminusM
Chris@42:     | (a, b) -> makeNode (CTimes (a, b))
Chris@42: 
Chris@42:   and sctimesjM = function 
Chris@42:     | (Uminus a, b) -> sctimesjM (a, b) >>= suminusM
Chris@42:     | (a, Uminus b) -> sctimesjM (a, b) >>= suminusM
Chris@42:     | (a, b) -> makeNode (CTimesJ (a, b))
Chris@42: 
Chris@42:   and reduce_sumM x = match x with
Chris@42:     [] -> returnM []
Chris@42:   | [Num a] -> 
Chris@42:       if (Number.is_zero a) then 
Chris@42: 	returnM [] 
Chris@42:       else returnM x
Chris@42:   | [Uminus (Num a)] -> 
Chris@42:       if (Number.is_zero a) then 
Chris@42: 	returnM [] 
Chris@42:       else returnM x
Chris@42:   | (Num a) :: (Num b) :: s -> 
Chris@42:       snumM (Number.add a b) >>= fun x ->
Chris@42: 	reduce_sumM (x :: s)
Chris@42:   | (Num a) :: (Uminus (Num b)) :: s -> 
Chris@42:       snumM (Number.sub a b) >>= fun x ->
Chris@42: 	reduce_sumM (x :: s)
Chris@42:   | (Uminus (Num a)) :: (Num b) :: s -> 
Chris@42:       snumM (Number.sub b a) >>= fun x ->
Chris@42: 	reduce_sumM (x :: s)
Chris@42:   | (Uminus (Num a)) :: (Uminus (Num b)) :: s -> 
Chris@42:       snumM (Number.add a b) >>= 
Chris@42:       suminusM >>= fun x ->
Chris@42: 	reduce_sumM (x :: s)
Chris@42:   | ((Num _) as a) :: b :: s -> reduce_sumM (b :: a :: s)
Chris@42:   | ((Uminus (Num _)) as a) :: b :: s -> reduce_sumM (b :: a :: s)
Chris@42:   | a :: s -> 
Chris@42:       reduce_sumM s >>= fun s' -> returnM (a :: s')
Chris@42: 
Chris@42:   and collectible1 = function
Chris@42:     | NaN _ -> false
Chris@42:     | Uminus x -> collectible1 x
Chris@42:     | _ -> true
Chris@42:   and collectible (a, b) = collectible1 a
Chris@42: 
Chris@42:   (* collect common factors: ax + bx -> (a+b)x *)
Chris@42:   and collectM which x = 
Chris@42:     let rec findCoeffM which = function
Chris@42:       |	Times (a, b) when collectible (which (a, b)) -> returnM (which (a, b))
Chris@42:       | Uminus x -> 
Chris@42: 	  findCoeffM which x >>= fun (coeff, b) ->
Chris@42: 	    suminusM coeff >>= fun mcoeff ->
Chris@42:  	      returnM (mcoeff, b)
Chris@42:       | x -> snumM Number.one >>= fun one -> returnM (one, x)
Chris@42:     and separateM xpr = function
Chris@42:  	[] -> returnM ([], [])
Chris@42:       |	a :: b ->
Chris@42:  	  separateM xpr b >>= fun (w, wo) ->
Chris@42: 	    (* try first factor *)
Chris@42:  	    findCoeffM (fun (a, b) -> (a, b)) a >>= fun (c, x) ->
Chris@42:  	      if (xpr == x) && collectible (c, x) then returnM (c :: w, wo)
Chris@42:  	      else
Chris@42: 	      (* try second factor *)
Chris@42:  		findCoeffM (fun (a, b) -> (b, a)) a >>= fun (c, x) ->
Chris@42:  		  if (xpr == x) && collectible (c, x) then returnM (c :: w, wo)
Chris@42:  		  else returnM (w, a :: wo)
Chris@42:     in match x with
Chris@42:       [] -> returnM x
Chris@42:     | [a] -> returnM x
Chris@42:     | a :: b ->
Chris@42:  	findCoeffM which a >>= fun (_, xpr) ->
Chris@42:  	  separateM xpr x >>= fun (w, wo) ->
Chris@42:  	    collectM which wo >>= fun wo' ->
Chris@42:  	      splusM w >>= fun w' ->
Chris@42:  		stimesM (w', xpr) >>= fun t' ->
Chris@42:  		  returnM (t':: wo')
Chris@42: 
Chris@42:   and mangleSumM x = returnM x
Chris@42:       >>= reduce_sumM 
Chris@42:       >>= collectM (fun (a, b) -> (a, b))
Chris@42:       >>= collectM (fun (a, b) -> (b, a))
Chris@42:       >>= reduce_sumM 
Chris@42:       >>= deepCollectM !Magic.deep_collect_depth
Chris@42:       >>= reduce_sumM
Chris@42: 
Chris@42:   and reorder_uminus = function  (* push all Uminuses to the end *)
Chris@42:       [] -> []
Chris@42:     | ((Uminus _) as a' :: b) -> (reorder_uminus b) @ [a']
Chris@42:     | (a :: b) -> a :: (reorder_uminus b)                      
Chris@42: 
Chris@42:   and canonicalizeM = function 
Chris@42:       [] -> snumM Number.zero
Chris@42:     | [a] -> makeNode a                    (* one term *)
Chris@42:     | a -> generateFusedMultAddM (reorder_uminus a)
Chris@42: 
Chris@42:   and generateFusedMultAddM = 
Chris@42:     let rec is_multiplication = function
Chris@42:       | Times (Num a, b) -> true
Chris@42:       | Uminus (Times (Num a, b)) -> true
Chris@42:       | _ -> false
Chris@42:     and separate = function
Chris@42: 	[] -> ([], [], Number.zero)
Chris@42:       | (Times (Num a, b)) as this :: c -> 
Chris@42: 	  let (x, y, max) = separate c in
Chris@42: 	  let newmax = if (Number.greater a max) then a else max in
Chris@42: 	  (this :: x, y, newmax)
Chris@42:       | (Uminus (Times (Num a, b))) as this :: c -> 
Chris@42: 	  let (x, y, max) = separate c in
Chris@42: 	  let newmax = if (Number.greater a max) then a else max in
Chris@42: 	  (this :: x, y, newmax)
Chris@42:       | this :: c ->
Chris@42: 	  let (x, y, max) = separate c in
Chris@42: 	  (x, this :: y, max)
Chris@42:     in fun l ->
Chris@42:       if !Magic.enable_fma && count is_multiplication l >= 2 then
Chris@42: 	let (w, wo, max) = separate l in
Chris@42: 	snumM (Number.div Number.one max) >>= fun invmax' ->
Chris@42: 	  snumM max >>= fun max' ->
Chris@42: 	    mapM (fun x -> stimesM (invmax', x)) w >>= splusM >>= fun pw' ->
Chris@42: 	      stimesM (max', pw') >>= fun mw' ->
Chris@42: 		splusM (wo @ [mw'])
Chris@42:       else 
Chris@42: 	makeNode (Plus l)
Chris@42: 
Chris@42: 
Chris@42:   and negative = function
Chris@42:       Uminus _ -> true
Chris@42:     | _ -> false
Chris@42: 
Chris@42:   (*
Chris@42:    * simplify patterns of the form
Chris@42:    *
Chris@42:    *  ((c_1 * a + ...) + ...) +  (c_2 * a + ...)
Chris@42:    *
Chris@42:    * The pattern includes arbitrary coefficients and minus signs.
Chris@42:    * A common case of this pattern is the butterfly
Chris@42:    *   (a + b) + (a - b)
Chris@42:    *   (a + b) - (a - b)
Chris@42:    *)
Chris@42:   (* this whole procedure needs much more thought *)
Chris@42:   and deepCollectM maxdepth l =
Chris@42:     let rec findTerms depth x = match x with
Chris@42:       | Uminus x -> findTerms depth x
Chris@42:       |	Times (Num _, b) -> (findTerms (depth - 1) b)
Chris@42:       |	Plus l when depth > 0 ->
Chris@42: 	  x :: List.flatten (List.map (findTerms (depth - 1)) l)
Chris@42:       |	x -> [x]
Chris@42:     and duplicates = function
Chris@42: 	[] -> []
Chris@42:       |	a :: b -> if List.memq a b then a :: duplicates b
Chris@42:       else duplicates b
Chris@42: 
Chris@42:     in let rec splitDuplicates depth d x =
Chris@42:       if (List.memq x d) then 
Chris@42: 	snumM (Number.zero) >>= fun zero ->
Chris@42: 	  returnM (zero, x)
Chris@42:       else match x with
Chris@42:       |	Times (a, b) ->
Chris@42: 	  splitDuplicates (depth - 1) d a >>= fun (a', xa) ->
Chris@42: 	    splitDuplicates (depth - 1) d b >>= fun (b', xb) ->
Chris@42: 	      stimesM (a', b') >>= fun ab ->
Chris@42: 		stimesM (a, xb) >>= fun xb' ->
Chris@42: 		  stimesM (xa, b) >>= fun xa' ->
Chris@42: 		    stimesM (xa, xb) >>= fun xab ->
Chris@42: 		      splusM [xa'; xb'; xab] >>= fun x ->
Chris@42: 			returnM (ab, x)
Chris@42:       | Uminus a -> 
Chris@42: 	  splitDuplicates depth d a >>= fun (x, y) ->
Chris@42: 	    suminusM x >>= fun ux -> 
Chris@42: 	      suminusM y >>= fun uy -> 
Chris@42: 		returnM (ux, uy)
Chris@42:       |	Plus l when depth > 0 -> 
Chris@42: 	  mapM (splitDuplicates (depth - 1) d) l >>= fun ld ->
Chris@42: 	    let (l', d') = List.split ld in
Chris@42: 	    splusM l' >>= fun p ->
Chris@42: 	      splusM d' >>= fun d'' ->
Chris@42: 	      returnM (p, d'')
Chris@42:       |	x -> 
Chris@42: 	  snumM (Number.zero) >>= fun zero' ->
Chris@42: 	    returnM (x, zero')
Chris@42: 
Chris@42:     in let l' = List.flatten (List.map (findTerms maxdepth) l)
Chris@42:     in match duplicates l' with
Chris@42:     | [] -> returnM l
Chris@42:     | d ->
Chris@42: 	mapM (splitDuplicates maxdepth d) l >>= fun ld ->
Chris@42: 	  let (l', d') = List.split ld in
Chris@42: 	  splusM l' >>= fun l'' ->
Chris@42: 	    let rec flattenPlusM = function
Chris@42: 	      | Plus l -> returnM l
Chris@42: 	      | Uminus x ->
Chris@42: 		  flattenPlusM x >>= mapM suminusM
Chris@42: 	      | x -> returnM [x]
Chris@42: 	    in
Chris@42: 	    mapM flattenPlusM d' >>= fun d'' ->
Chris@42: 	      splusM (List.flatten d'') >>= fun d''' ->
Chris@42: 		mangleSumM [l''; d''']
Chris@42: 
Chris@42:   and splusM l =
Chris@42:     let fma_heuristics x = 
Chris@42:       if !Magic.enable_fma then 
Chris@42: 	match x with
Chris@42: 	| [Uminus (Times _); Times _] -> Some false
Chris@42: 	| [Times _; Uminus (Times _)] -> Some false
Chris@42: 	| [Uminus (_); Times _] -> Some true
Chris@42: 	| [Times _; Uminus (Plus _)] -> Some true
Chris@42: 	| [_; Uminus (Times _)] -> Some false
Chris@42: 	| [Uminus (Times _); _] -> Some false
Chris@42: 	| _ -> None
Chris@42:       else
Chris@42: 	None
Chris@42:     in
Chris@42:     mangleSumM l >>=  fun l' ->
Chris@42:       (* no terms are negative.  Don't do anything *)
Chris@42:       if not (List.exists negative l') then
Chris@42: 	canonicalizeM l'
Chris@42:       (* all terms are negative.  Negate them all and collect the minus sign *)
Chris@42:       else if List.for_all negative l' then
Chris@42: 	mapM suminusM l' >>= splusM >>= suminusM
Chris@42:       else match fma_heuristics l' with
Chris@42:       |	Some true -> mapM suminusM l' >>= splusM >>= suminusM
Chris@42:       |	Some false -> canonicalizeM l'
Chris@42:       |	None ->
Chris@42:          (* Ask the Oracle for the canonical form *)
Chris@42: 	  if (not !Magic.randomized_cse) &&
Chris@42: 	    Oracle.should_flip_sign (Plus l') then
Chris@42: 	    mapM suminusM l' >>= splusM >>= suminusM
Chris@42: 	  else
Chris@42: 	    canonicalizeM l'
Chris@42: 
Chris@42:   (* monadic style algebraic simplifier for the dag *)
Chris@42:   let rec algsimpM x =
Chris@42:     memoizing lookupSimpM insertSimpM 
Chris@42:       (function 
Chris@42:  	| Num a -> snumM a
Chris@42:  	| NaN _ as x -> makeNode x
Chris@42:  	| Plus a -> 
Chris@42:  	    mapM algsimpM a >>= splusM
Chris@42:  	| Times (a, b) -> 
Chris@42:  	    (algsimpM a >>= fun a' ->
Chris@42:  	      algsimpM b >>= fun b' ->
Chris@42:  		stimesM (a', b'))
Chris@42:  	| CTimes (a, b) -> 
Chris@42:  	    (algsimpM a >>= fun a' ->
Chris@42:  	      algsimpM b >>= fun b' ->
Chris@42: 		sctimesM (a', b'))
Chris@42:  	| CTimesJ (a, b) -> 
Chris@42:  	    (algsimpM a >>= fun a' ->
Chris@42:  	      algsimpM b >>= fun b' ->
Chris@42: 		sctimesjM (a', b'))
Chris@42:  	| Uminus a -> 
Chris@42:  	    algsimpM a >>= suminusM 
Chris@42:  	| Store (v, a) ->
Chris@42:  	    algsimpM a >>= fun a' ->
Chris@42:  	      makeNode (Store (v, a'))
Chris@42:  	| Load _ as x -> makeNode x)
Chris@42:       x
Chris@42: 
Chris@42:    let initialTable = (empty, empty)
Chris@42:    let simp_roots = mapM algsimpM
Chris@42:    let algsimp = runM initialTable simp_roots
Chris@42: end
Chris@42: 
Chris@42: (*************************************************************
Chris@42:  * Network transposition algorithm
Chris@42:  *************************************************************)
Chris@42: module Transpose = struct
Chris@42:   open Monads.StateMonad
Chris@42:   open Monads.MemoMonad
Chris@42:   open Littlesimp
Chris@42: 
Chris@42:   let fetchDuals = fetchState
Chris@42:   let storeDuals = storeState
Chris@42: 
Chris@42:   let lookupDualsM key =
Chris@42:     fetchDuals >>= fun table ->
Chris@42:       returnM (node_lookup key table)
Chris@42: 
Chris@42:   let insertDualsM key value =
Chris@42:     fetchDuals >>= fun table ->
Chris@42:       storeDuals (node_insert key value table)
Chris@42: 
Chris@42:   let rec visit visited vtable parent_table = function
Chris@42:       [] -> (visited, parent_table)
Chris@42:     | node :: rest ->
Chris@42: 	match node_lookup node vtable with
Chris@42: 	| Some _ -> visit visited vtable parent_table rest
Chris@42: 	| None ->
Chris@42: 	    let children = match node with
Chris@42: 	    | Store (v, n) -> [n]
Chris@42: 	    | Plus l -> l
Chris@42: 	    | Times (a, b) -> [a; b]
Chris@42: 	    | CTimes (a, b) -> [a; b]
Chris@42: 	    | CTimesJ (a, b) -> [a; b]
Chris@42: 	    | Uminus x -> [x]
Chris@42: 	    | _ -> []
Chris@42: 	    in let rec loop t = function
Chris@42: 		[] -> t
Chris@42: 	      |	a :: rest ->
Chris@42: 		  (match node_lookup a t with
Chris@42: 		    None -> loop (node_insert a [node] t) rest
Chris@42: 		  | Some c -> loop (node_insert a (node :: c) t) rest)
Chris@42: 	    in 
Chris@42: 	    (visit 
Chris@42: 	       (node :: visited)
Chris@42: 	       (node_insert node () vtable)
Chris@42: 	       (loop parent_table children)
Chris@42: 	       (children @ rest))
Chris@42: 
Chris@42:   let make_transposer parent_table =
Chris@42:     let rec termM node candidate_parent = 
Chris@42:       match candidate_parent with
Chris@42:       |	Store (_, n) when n == node -> 
Chris@42: 	  dualM candidate_parent >>= fun x' -> returnM [x']
Chris@42:       | Plus (l) when List.memq node l -> 
Chris@42: 	  dualM candidate_parent >>= fun x' -> returnM [x']
Chris@42:       | Times (a, b) when b == node -> 
Chris@42: 	  dualM candidate_parent >>= fun x' -> 
Chris@42: 	    returnM [makeTimes (a, x')]
Chris@42:       | CTimes (a, b) when b == node -> 
Chris@42: 	  dualM candidate_parent >>= fun x' -> 
Chris@42: 	    returnM [CTimes (a, x')]
Chris@42:       | CTimesJ (a, b) when b == node -> 
Chris@42: 	  dualM candidate_parent >>= fun x' -> 
Chris@42: 	    returnM [CTimesJ (a, x')]
Chris@42:       | Uminus n when n == node -> 
Chris@42: 	  dualM candidate_parent >>= fun x' -> 
Chris@42: 	    returnM [makeUminus x']
Chris@42:       | _ -> returnM []
Chris@42:     
Chris@42:     and dualExpressionM this_node = 
Chris@42:       mapM (termM this_node) 
Chris@42: 	(match node_lookup this_node parent_table with
Chris@42: 	| Some a -> a
Chris@42: 	| None -> failwith "bug in dualExpressionM"
Chris@42: 	) >>= fun l ->
Chris@42: 	returnM (makePlus (List.flatten l))
Chris@42: 
Chris@42:     and dualM this_node =
Chris@42:       memoizing lookupDualsM insertDualsM
Chris@42: 	(function
Chris@42: 	  | Load v as x -> 
Chris@42: 	      if (Variable.is_constant v) then
Chris@42: 		returnM (Load v)
Chris@42: 	      else
Chris@42: 		(dualExpressionM x >>= fun d ->
Chris@42: 		  returnM (Store (v, d)))
Chris@42: 	  | Store (v, x) -> returnM (Load v)
Chris@42: 	  | x -> dualExpressionM x)
Chris@42: 	this_node
Chris@42: 
Chris@42:     in dualM
Chris@42: 
Chris@42:   let is_store = function 
Chris@42:     | Store _ -> true
Chris@42:     | _ -> false
Chris@42: 
Chris@42:   let transpose dag = 
Chris@42:     let _ = Util.info "begin transpose" in
Chris@42:     let (all_nodes, parent_table) = 
Chris@42:       visit [] Assoctable.empty Assoctable.empty dag in
Chris@42:     let transposerM = make_transposer parent_table in
Chris@42:     let mapTransposerM = mapM transposerM in
Chris@42:     let duals = runM Assoctable.empty mapTransposerM all_nodes in
Chris@42:     let roots = List.filter is_store duals in
Chris@42:     let _ = Util.info "end transpose" in
Chris@42:     roots
Chris@42: end
Chris@42: 
Chris@42: 
Chris@42: (*************************************************************
Chris@42:  * Various dag statistics
Chris@42:  *************************************************************)
Chris@42: module Stats : sig
Chris@42:   type complexity
Chris@42:   val complexity : Expr.expr list -> complexity
Chris@42:   val same_complexity : complexity -> complexity -> bool
Chris@42:   val leq_complexity : complexity -> complexity -> bool
Chris@42:   val to_string : complexity -> string
Chris@42: end = struct
Chris@42:   type complexity = int * int * int * int * int * int
Chris@42:   let rec visit visited vtable = function
Chris@42:       [] -> visited
Chris@42:     | node :: rest ->
Chris@42: 	match node_lookup node vtable with
Chris@42: 	  Some _ -> visit visited vtable rest
Chris@42: 	| None ->
Chris@42: 	    let children = match node with
Chris@42: 	      Store (v, n) -> [n]
Chris@42: 	    | Plus l -> l
Chris@42: 	    | Times (a, b) -> [a; b]
Chris@42: 	    | Uminus x -> [x]
Chris@42: 	    | _ -> []
Chris@42: 	    in visit (node :: visited)
Chris@42: 	      (node_insert node () vtable)
Chris@42: 	      (children @ rest)
Chris@42: 
Chris@42:   let complexity dag = 
Chris@42:     let rec loop (load, store, plus, times, uminus, num) = function 
Chris@42:       	[] -> (load, store, plus, times, uminus, num)
Chris@42:       | node :: rest ->
Chris@42: 	  loop
Chris@42: 	    (match node with
Chris@42: 	    | Load _ -> (load + 1, store, plus, times, uminus, num)
Chris@42: 	    | Store _ -> (load, store + 1, plus, times, uminus, num)
Chris@42: 	    | Plus x -> (load, store, plus + (List.length x - 1), times, uminus, num)
Chris@42: 	    | Times _ -> (load, store, plus, times + 1, uminus, num)
Chris@42: 	    | Uminus _ -> (load, store, plus, times, uminus + 1, num)
Chris@42: 	    | Num _ -> (load, store, plus, times, uminus, num + 1)
Chris@42: 	    | CTimes _ -> (load, store, plus, times, uminus, num)
Chris@42: 	    | CTimesJ _ -> (load, store, plus, times, uminus, num)
Chris@42: 	    | NaN _ -> (load, store, plus, times, uminus, num))
Chris@42: 	    rest
Chris@42:     in let (l, s, p, t, u, n) = 
Chris@42:       loop (0, 0, 0, 0, 0, 0) (visit [] Assoctable.empty dag)
Chris@42:     in (l, s, p, t, u, n)
Chris@42: 
Chris@42:   let weight (l, s, p, t, u, n) =
Chris@42:     l + s + 10 * p + 20 * t + u + n
Chris@42: 
Chris@42:   let same_complexity a b = weight a = weight b
Chris@42:   let leq_complexity a b = weight a <= weight b
Chris@42: 
Chris@42:   let to_string (l, s, p, t, u, n) =
Chris@42:     Printf.sprintf "ld=%d st=%d add=%d mul=%d uminus=%d num=%d\n"
Chris@42: 		   l s p t u n
Chris@42: 		   
Chris@42: end    
Chris@42: 
Chris@42: (* simplify the dag *)
Chris@42: let algsimp v = 
Chris@42:   let rec simplification_loop v =
Chris@42:     let () = Util.info "simplification step" in
Chris@42:     let complexity = Stats.complexity v in
Chris@42:     let () = Util.info ("complexity = " ^ (Stats.to_string complexity)) in
Chris@42:     let v = (AlgSimp.algsimp @@ Transpose.transpose @@ 
Chris@42: 	     AlgSimp.algsimp @@ Transpose.transpose) v in
Chris@42:     let complexity' = Stats.complexity v in
Chris@42:     let () = Util.info ("complexity = " ^ (Stats.to_string complexity')) in
Chris@42:     if (Stats.leq_complexity complexity' complexity) then
Chris@42:       let () = Util.info "end algsimp" in
Chris@42:       v
Chris@42:     else
Chris@42:       simplification_loop v
Chris@42: 
Chris@42:   in
Chris@42:   let () = Util.info "begin algsimp" in
Chris@42:   let v = AlgSimp.algsimp v in
Chris@42:   if !Magic.network_transposition then simplification_loop v else v
Chris@42: