diff --git a/src/experimental/graph/core.nim b/src/experimental/graph/core.nim
index 31e57ea..e322488 100644
--- a/src/experimental/graph/core.nim
+++ b/src/experimental/graph/core.nim
@@ -26,6 +26,7 @@ type
     tag: Gtags
     inputs*: seq[Gvalue]
     gfunc*: Gfunc
+    locals*: seq[Gvalue]  ## for sharing values between forward and among backward functions
     epoch: int
 
 type
@@ -75,7 +76,7 @@ proc nodeRepr*(x: Gvalue): string =
 method newOneOf*(x: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("newOneOf(" & $x & ")")  ## Be sure to zero init fields
 method valCopy*(z: Gvalue, x: Gvalue) {.base.} = raiseErrorBaseMethod("valCopy(" & $z & "," & $x & ")")
 
-method isTrue*(x: Gvalue): bool {.base.} = raiseErrorBaseMethod("isTrue(" & $x & ")")
+method isZero*(x: Gvalue): bool {.base.} = raiseErrorBaseMethod("isZero(" & $x & ")")
 method update*(x: Gvalue, y: int) {.base.} = raiseErrorBaseMethod("update(" & $x & "," & $y & ")")
 method update*(x: Gvalue, y: float) {.base.} = raiseErrorBaseMethod("update(" & $x & "," & $y & ")")
 
@@ -182,10 +183,10 @@ proc condb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
     raiseValueError("i must be 0 or 1, got: " & $i)
 
 proc condf(v: Gvalue) =
-  if v.inputs[0].isTrue:
-    v.valCopy v.inputs[1]
-  else:
+  if v.inputs[0].isZero:
     v.valCopy v.inputs[2]
+  else:
+    v.valCopy v.inputs[1]
 
 let gcond = newGfunc(forward = condf, backward = condb, name = "cond")
 
@@ -200,6 +201,14 @@ proc updated*(x: Gvalue) =
   inc epoch
   x.epoch = epoch
 
+proc evaluated*(x: Gvalue) =
+  ## signal up-to-date value, given inputs, useful for update value outside of eval, as in update to locals in forward
+  var maxep = 0
+  for i in x.inputs:
+    if maxep < i.epoch:
+      maxep = i.epoch
+  x.epoch = maxep
+
 proc eval*(v: Gvalue): Gvalue {.discardable.} =
   proc r(x: Gvalue) =
     if gtVisited in x.tag:
@@ -210,14 +219,14 @@ proc eval*(v: Gvalue): Gvalue {.discardable.} =
       x.inputs[0].r
       if maxep < x.inputs[0].epoch:
         maxep = x.inputs[0].epoch
-      if x.inputs[0].isTrue:
-        x.inputs[1].r
-        if maxep < x.inputs[1].epoch:
-          maxep = x.inputs[1].epoch
-      else:
+      if x.inputs[0].isZero:
         x.inputs[2].r
         if maxep < x.inputs[2].epoch:
           maxep = x.inputs[2].epoch
+      else:
+        x.inputs[1].r
+        if maxep < x.inputs[1].epoch:
+          maxep = x.inputs[1].epoch
     else:
       for i in x.inputs:
         i.r
diff --git a/src/experimental/graph/gauge.nim b/src/experimental/graph/gauge.nim
index 62ac32b..84f216a 100644
--- a/src/experimental/graph/gauge.nim
+++ b/src/experimental/graph/gauge.nim
@@ -11,42 +11,45 @@ import layout, ../../gauge, physics/qcdTypes
 type Gauge = seq[DLatticeColorMatrixV]
 
 type Ggauge* {.final.} = ref object of Gvalue
+  isZero: bool = false  ## specialized for zero fields, unrelated to actual gval
   gval: Gauge
 
 proc getgauge*(x: Gvalue): Gauge = Ggauge(x).gval
 
-proc update*(x: Gvalue, g: Gauge) =
-  let u = x.getgauge
-  threads:
-    for mu in 0..<u.len:
-      u[mu] := g[mu]
+proc update*(x: Gvalue, g: Gauge, isZero = false) =
+  let u = Ggauge(x)
+  u.isZero = isZero
+  if not isZero:
+    threads:
+      for mu in 0..<u.gval.len:
+        u.gval[mu] := g[mu]
   x.updated
 
-proc toGvalue*(x: Gauge): Ggauge =
+proc toGvalue*(x: Gauge, isZero = false): Ggauge =
   # proc instead of converter to avoid converting seq
-  result = Ggauge(gval: x)
+  result = Ggauge(isZero: isZero, gval: x)
   result.updated
 
 method newOneOf*(x: Ggauge): Gvalue =
   let g = x.gval.newOneOf
-  threads:
-    for f in g:
-      f := 0.0
-  Ggauge(gval: g)
+  Ggauge(isZero: true, gval: g)
 method valCopy*(z: Ggauge, x: Ggauge) =
   let u = z.gval
   let v = x.gval
-  threads:
-    for mu in 0..<u.len:
-      u[mu] := v[mu]
+  z.isZero = x.isZero
+  if not x.isZero:
+    threads:
+      for mu in 0..<u.len:
+        u[mu] := v[mu]
 
 method `$`*(x: Ggauge): string =
   let v = x.gval[0][0][0,0]
-  "Gauge (" & $v.re[0] & ", " & $v.im[0] & ")"
+  if x.isZero:
+    result = "Gauge=0 (" & $v.re[0] & ", " & $v.im[0] & ")"
+  else:
+    result = "Gauge (" & $v.re[0] & ", " & $v.im[0] & ")"
 
-#
-# basic ops
-#
+method isZero*(x: Ggauge): bool = x.isZero
 
 method retr*(x: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("retr(" & $x & ")")
 method adj*(x: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("adj(" & $x & ")")
@@ -55,8 +58,18 @@ method redot*(x: Gvalue, y: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("red
 method expDeriv*(b: Gvalue, x: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("expDeriv(" & $b & "," & $x & ")")
 method projTAH*(x: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("projTAH(" & $x & ")")
 
+method adjmul*(x: Gvalue, y: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("adjmul(" & $x & "," & $y & ")")
+method muladj*(x: Gvalue, y: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("muladj(" & $x & "," & $y & ")")
+method contractProjTAH*(x: Gvalue, y: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("contractProjTAH(" & $x & "," & $y & ")")
+method axexp*(a: Gvalue, x: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("axexp(" & $a & "," & $x & ")")
+method axexpmuly*(a: Gvalue, x: Gvalue, y: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("axexpmuly(" & $a & "," & $x & "," & $y & ")")
+
 method redot*(x: Gscalar, y: Gscalar): Gvalue = x*y
 
+#
+# basic ops
+#
+
 proc retrgb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
   let g = z.inputs[0].getgauge.newOneOf
   threads:
@@ -75,11 +88,14 @@ proc retrgb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
 proc retrgf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let z = Gscalar(v)
-  threads:
-    var t = 0.0
-    for mu in 0..<x.gval.len:
-      t += x.gval[mu].trace.re
-    threadMaster: z.getfloat = t
+  if x.isZero:
+    z.getfloat = 0.0
+  else:
+    threads:
+      var t = 0.0
+      for mu in 0..<x.gval.len:
+        t += x.gval[mu].trace.re
+      threadMaster: z.getfloat = t
 
 let retrg = newGfunc(forward = retrgf, backward = retrgb, name = "retrg")
 
@@ -97,9 +113,15 @@ proc adjgb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
 proc adjgf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      z.gval[mu] := x.gval[mu].adj
+  if x.isZero:
+    if not z.isZero:
+      z.isZero = true
+  else:
+    threads:
+      for mu in 0..<z.gval.len:
+        z.gval[mu] := x.gval[mu].adj
+    if z.isZero:
+      z.isZero = false
 
 let adjg = newGfunc(forward = adjgf, backward = adjgb, name = "adjg")
 
@@ -118,11 +140,14 @@ proc norm2gb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
 proc norm2gf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let z = Gscalar(v)
-  threads:
-    var t = 0.0
-    for mu in 0..<x.gval.len:
-      t += x.gval[mu].norm2
-    threadMaster: z.getfloat = t
+  if x.isZero:
+    z.getfloat = 0.0
+  else:
+    threads:
+      var t = 0.0
+      for mu in 0..<x.gval.len:
+        t += x.gval[mu].norm2
+      threadMaster: z.getfloat = t
 
 let norm2g = newGfunc(forward = norm2gf, backward = norm2gb, name = "norm2g")
 
@@ -131,9 +156,15 @@ method norm2*(x: Ggauge): Gvalue = Gscalar(inputs: @[Gvalue(x)], gfunc: norm2g)
 proc neggf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      z.gval[mu] := -x.gval[mu]
+  if x.isZero:
+    if not z.isZero:
+      z.isZero = true
+  else:
+    threads:
+      for mu in 0..<z.gval.len:
+        z.gval[mu] := -x.gval[mu]
+    if z.isZero:
+      z.isZero = false
 
 proc neggb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
   if zb == nil:
@@ -163,9 +194,27 @@ proc addsgf(v: Gvalue) =
   let x = Gscalar(v.inputs[0])
   let y = Ggauge(v.inputs[1])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      z.gval[mu] := x.getfloat + y.gval[mu]
+  if x.isZero:
+    if y.isZero:
+      if not z.isZero:
+        z.isZero = true
+    else:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := y.gval[mu]
+      if z.isZero:
+        z.isZero = false
+  else:
+    if y.isZero:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := x.getfloat
+    else:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := x.getfloat + y.gval[mu]
+    if z.isZero:
+      z.isZero = false
 
 let addsg = newGfunc(forward = addsgf, backward = addsgb, name = "s+g")
 
@@ -181,9 +230,27 @@ proc addggf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let y = Ggauge(v.inputs[1])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      z.gval[mu] := x.gval[mu] + y.gval[mu]
+  if x.isZero:
+    if y.isZero:
+      if not z.isZero:
+        z.isZero = true
+    else:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := y.gval[mu]
+      if z.isZero:
+        z.isZero = false
+  else:
+    if y.isZero:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := x.gval[mu]
+    else:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := x.gval[mu] + y.gval[mu]
+    if z.isZero:
+      z.isZero = false
 
 let addgg = newGfunc(forward = addggf, backward = addggb, name = "g+g")
 
@@ -204,9 +271,15 @@ proc mulsgf(v: Gvalue) =
   let x = Gscalar(v.inputs[0])
   let y = Ggauge(v.inputs[1])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      z.gval[mu] := x.getfloat * y.gval[mu]
+  if x.isZero or y.isZero:
+    if not z.isZero:
+      z.isZero = true
+  else:
+    threads:
+      for mu in 0..<z.gval.len:
+        z.gval[mu] := x.getfloat * y.gval[mu]
+    if z.isZero:
+      z.isZero = false
 
 let mulsg = newGfunc(forward = mulsgf, backward = mulsgb, name = "s*g")
 
@@ -217,9 +290,9 @@ proc mulggb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
     raiseValueError("no backward node")
   case i
   of 0:
-    return zb*z.inputs[1].adj
+    return zb.muladj z.inputs[1]
   of 1:
-    return z.inputs[0].adj*zb
+    return z.inputs[0].adjmul zb
   else:
     raiseValueError("i must be 0 or 1, got: " & $i)
 
@@ -227,9 +300,15 @@ proc mulggf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let y = Ggauge(v.inputs[1])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      z.gval[mu] := x.gval[mu] * y.gval[mu]
+  if x.isZero or y.isZero:
+    if not z.isZero:
+      z.isZero = true
+  else:
+    threads:
+      for mu in 0..<z.gval.len:
+        z.gval[mu] := x.gval[mu] * y.gval[mu]
+    if z.isZero:
+      z.isZero = false
 
 let mulgg = newGfunc(forward = mulggf, backward = mulggb, name = "g*g")
 
@@ -254,11 +333,15 @@ proc redotggf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let y = Ggauge(v.inputs[1])
   let z = Gscalar(v)
-  threads:
-    var t = 0.0
-    for mu in 0..<x.gval.len:
-      t += redot(x.gval[mu], y.gval[mu])
-    threadMaster: z.getfloat = t
+  if x.isZero or y.isZero:
+    if not z.isZero:
+      z.getfloat = 0.0
+  else:
+    threads:
+      var t = 0.0
+      for mu in 0..<x.gval.len:
+        t += redot(x.gval[mu], y.gval[mu])
+      threadMaster: z.getfloat = t
 
 let redotgg = newGfunc(forward = redotggf, backward = redotggb, name = "redotgg")
 
@@ -279,9 +362,27 @@ proc subgsf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let y = Gscalar(v.inputs[1])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      z.gval[mu] := x.gval[mu] - y.getfloat
+  if x.isZero:
+    if y.isZero:
+      if not z.isZero:
+        z.isZero = true
+    else:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := -y.getfloat
+      if z.isZero:
+        z.isZero = false
+  else:
+    if y.isZero:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := x.gval[mu]
+    else:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := x.gval[mu] - y.getfloat
+    if z.isZero:
+      z.isZero = false
 
 let subgs = newGfunc(forward = subgsf, backward = subgsb, name = "g-s")
 
@@ -302,9 +403,27 @@ proc subggf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let y = Ggauge(v.inputs[1])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      z.gval[mu] := x.gval[mu] - y.gval[mu]
+  if x.isZero:
+    if y.isZero:
+      if not z.isZero:
+        z.isZero = true
+    else:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := -y.gval[mu]
+      if z.isZero:
+        z.isZero = false
+  else:
+    if y.isZero:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := x.gval[mu]
+    else:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := x.gval[mu] - y.gval[mu]
+    if z.isZero:
+      z.isZero = false
 
 let subgg = newGfunc(forward = subggf, backward = subggb, name = "g-g")
 
@@ -322,10 +441,17 @@ proc expgb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
 proc expgf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      for e in z.gval[mu]:
-        z.gval[mu][e] := exp(x.gval[mu][e])
+  if x.isZero:
+    threads:
+      for mu in 0..<z.gval.len:
+        z.gval[mu] := 1.0
+  else:
+    threads:
+      for mu in 0..<z.gval.len:
+        for e in z.gval[mu]:
+          z.gval[mu][e] := exp(x.gval[mu][e])
+    if z.isZero:
+      z.isZero = false
 
 let expg = newGfunc(forward = expgf, backward = expgb, name = "expg")
 
@@ -344,10 +470,23 @@ proc expDerivgf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let y = Ggauge(v.inputs[1])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      for e in z.gval[mu]:
-        z.gval[mu][e] := expDeriv(y.gval[mu][e], x.gval[mu][e])
+  if x.isZero:
+    if not z.isZero:
+      z.isZero = true
+  elif y.isZero:
+    threads:
+      for mu in 0..<z.gval.len:
+        for e in z.gval[mu]:
+          z.gval[mu][e] := x.gval[mu][e]
+    if z.isZero:
+      z.isZero = false
+  else:
+    threads:
+      for mu in 0..<z.gval.len:
+        for e in z.gval[mu]:
+          z.gval[mu][e] := expDeriv(y.gval[mu][e], x.gval[mu][e])
+    if z.isZero:
+      z.isZero = false
 
 let expDerivg = newGfunc(forward = expDerivgf, backward = expDerivgb, name = "expDerivg")
 
@@ -365,10 +504,16 @@ proc projTAHb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
 proc projTAHf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      for e in z.gval[mu]:
-        z.gval[mu][e].projectTAH(x.gval[mu][e])
+  if x.isZero:
+    if not z.isZero:
+      z.isZero = true
+  else:
+    threads:
+      for mu in 0..<z.gval.len:
+        for e in z.gval[mu]:
+          z.gval[mu][e].projectTAH(x.gval[mu][e])
+    if z.isZero:
+      z.isZero = false
 
 let projTAHg = newGfunc(forward = projTAHf, backward = projTAHb, name = "projTAH")
 
@@ -378,25 +523,74 @@ method projTAH*(x: Ggauge): Gvalue = Ggauge(gval: x.gval.newOneOf, inputs: @[Gva
 # Fused ops
 #
 
-#[ TODO:
-- reuse common terms in the backward
-  - save & reuse local vars attached to a node
-- fuse the backward
-]#
+proc adjmulggb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
+  if zb == nil:
+    raiseValueError("no backward node")
+  case i
+  of 0:
+    return z.inputs[1].muladj zb
+  of 1:
+    return z.inputs[0] * zb
+  else:
+    raiseValueError("i must be 0 or 1, got: " & $i)
 
-method contractProjTAH*(x: Gvalue, y: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("contractProjTAH(" & $x & "," & $y & ")")
-method axexp*(a: Gvalue, x: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("axexp(" & $a & "," & $x & ")")
-method axexpmuly*(a: Gvalue, x: Gvalue, y: Gvalue): Gvalue {.base.} = raiseErrorBaseMethod("axexpmuly(" & $a & "," & $x & "," & $y & ")")
+proc adjmulggf(v: Gvalue) =
+  let x = Ggauge(v.inputs[0])
+  let y = Ggauge(v.inputs[1])
+  let z = Ggauge(v)
+  if x.isZero or y.isZero:
+    if not z.isZero:
+      z.isZero = true
+  else:
+    threads:
+      for mu in 0..<z.gval.len:
+        z.gval[mu] := x.gval[mu].adj * y.gval[mu]
+    if z.isZero:
+      z.isZero = false
+
+let adjmulgg = newGfunc(forward = adjmulggf, backward = adjmulggb, name = "g.adj*g")
+
+method adjmul*(x: Ggauge, y: Ggauge): Gvalue = Ggauge(gval: x.gval.newOneOf, inputs: @[Gvalue(x), y], gfunc: adjmulgg)
+
+proc muladjggb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
+  if zb == nil:
+    raiseValueError("no backward node")
+  case i
+  of 0:
+    return zb * z.inputs[1]
+  of 1:
+    return zb.adjmul z.inputs[0]
+  else:
+    raiseValueError("i must be 0 or 1, got: " & $i)
+
+proc muladjggf(v: Gvalue) =
+  let x = Ggauge(v.inputs[0])
+  let y = Ggauge(v.inputs[1])
+  let z = Ggauge(v)
+  if x.isZero or y.isZero:
+    if not z.isZero:
+      z.isZero = true
+  else:
+    threads:
+      for mu in 0..<z.gval.len:
+        z.gval[mu] := x.gval[mu] * y.gval[mu].adj
+    if z.isZero:
+      z.isZero = false
+
+let muladjgg = newGfunc(forward = muladjggf, backward = muladjggb, name = "g*g.adj")
+
+method muladj*(x: Ggauge, y: Ggauge): Gvalue = Ggauge(gval: x.gval.newOneOf, inputs: @[Gvalue(x), y], gfunc: muladjgg)
 
 proc contractProjTAHb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
   if zb == nil:
     raiseValueError("no backward node")
-  let zbTAH = zb.projTAH
+  if z.locals[0] == nil:
+    z.locals[0] = zb.projTAH
   case i
   of 0:
-    return zbTAH * z.inputs[1]
+    return z.locals[0] * z.inputs[1]
   of 1:
-    return zbTAH.adj * z.inputs[0]
+    return z.locals[0].adjmul z.inputs[0]
   else:
     raiseValueError("i must be 0 or 1, got: " & $i)
 
@@ -404,15 +598,21 @@ proc contractProjTAHf(v: Gvalue) =
   let x = Ggauge(v.inputs[0])
   let y = Ggauge(v.inputs[1])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      for e in z.gval[mu]:
-        let s = x.gval[mu][e]*y.gval[mu][e].adj
-        z.gval[mu][e].projectTAH s
+  if x.isZero or y.isZero:
+    if not z.isZero:
+      z.isZero = true
+  else:
+    threads:
+      for mu in 0..<z.gval.len:
+        for e in z.gval[mu]:
+          let s = x.gval[mu][e]*y.gval[mu][e].adj
+          z.gval[mu][e].projectTAH s
+    if z.isZero:
+      z.isZero = false
 
 let contractProjTAHg = newGfunc(forward = contractProjTAHf, backward = contractProjTAHb, name = "projTAHg")
 
-method contractProjTAH*(x: Ggauge, y: Ggauge): Gvalue = Ggauge(gval: x.gval.newOneOf, inputs: @[Gvalue(x), y], gfunc: contractProjTAHg)
+method contractProjTAH*(x: Ggauge, y: Ggauge): Gvalue = Ggauge(gval: x.gval.newOneOf, inputs: @[Gvalue(x), y], gfunc: contractProjTAHg, locals: newseq[Gvalue](1))
 
 proc axexpb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
   # z = exp(ax)
@@ -421,27 +621,36 @@ proc axexpb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
   # Dx = a exp'(Dz, ax)
   if zb == nil:
     raiseValueError("no backward node")
-  let dexp = expDeriv(zb, z.inputs[0]*z.inputs[1])
+  if z.locals[0] == nil:
+    z.locals[0] = expDeriv(zb, z.inputs[0]*z.inputs[1])
   case i
   of 0:
-    return dexp.redot z.inputs[1]
+    return z.locals[0].redot z.inputs[1]
   of 1:
-    return z.inputs[0]*dexp
+    return z.inputs[0]*z.locals[0]
   else:
     raiseValueError("i must be 0 or 1, got: " & $i)
 
 proc axexpf(v: Gvalue) =
-  let a = Gscalar(v.inputs[0]).getfloat
+  let a = Gscalar(v.inputs[0])
   let x = Ggauge(v.inputs[1])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      for e in z.gval[mu]:
-        z.gval[mu][e] := exp(a*x.gval[mu][e])
+  if a.isZero or x.isZero:
+    threads:
+      for mu in 0..<z.gval.len:
+        z.gval[mu] := 1.0
+  else:
+    let f = a.getfloat
+    threads:
+      for mu in 0..<z.gval.len:
+        for e in z.gval[mu]:
+          z.gval[mu][e] := exp(f*x.gval[mu][e])
+  if z.isZero:
+    z.isZero = false
 
 let axexpg = newGfunc(forward = axexpf, backward = axexpb, name = "axexp")
 
-method axexp*(a: Gscalar, x: Ggauge): Gvalue = Ggauge(gval: x.gval.newOneOf, inputs: @[Gvalue(a), x], gfunc: axexpg)
+method axexp*(a: Gscalar, x: Ggauge): Gvalue = Ggauge(gval: x.gval.newOneOf, inputs: @[Gvalue(a), x], gfunc: axexpg, locals: newseq[Gvalue](1))
 
 proc axexpmulyb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
   # z = exp(ax)y
@@ -452,14 +661,18 @@ proc axexpmulyb(zb: Gvalue, z: Gvalue, i: int, dep: Gvalue): Gvalue =
   # Dy = exp(ax)† Dz
   if zb == nil:
     raiseValueError("no backward node")
-  let dexp = expDeriv(zb*z.inputs[2].adj, z.inputs[0]*z.inputs[1])
+  if i == 0 or i == 1:
+    if z.locals[1] == nil:
+      z.locals[1] = expDeriv(zb.muladj z.inputs[2], z.inputs[0]*z.inputs[1])
   case i
   of 0:
-    return dexp.redot z.inputs[1]
+    return z.locals[1].redot z.inputs[1]
   of 1:
-    return z.inputs[0]*dexp
+    return z.inputs[0]*z.locals[1]
   of 2:
-    return axexp(z.inputs[0], z.inputs[1]).adj*zb
+    if z.locals[0] == nil:
+      z.locals[0] = axexp(z.inputs[0], z.inputs[1])
+    return z.locals[0].adjmul zb
   else:
     raiseValueError("i must be 0, 1, or 2, got: " & $i)
 
@@ -468,16 +681,44 @@ proc axexpmulyf(v: Gvalue) =
   let x = Ggauge(v.inputs[1])
   let y = Ggauge(v.inputs[2])
   let z = Ggauge(v)
-  threads:
-    for mu in 0..<z.gval.len:
-      for e in z.gval[mu]:
-        var t{.noinit.}: evalType(x.gval[mu][e])
-        t := exp(a*x.gval[mu][e])
-        z.gval[mu][e] := t*y.gval[mu][e]
+  if y.isZero:
+    if not z.isZero:
+      z.isZero = true
+  elif a.isZero or x.isZero:
+    threads:
+      for mu in 0..<z.gval.len:
+        z.gval[mu] := y.gval[mu]
+    if z.isZero:
+      z.isZero = false
+  else:
+    let f = a.getfloat
+    if z.locals[0] == nil:
+      threads:
+        for mu in 0..<z.gval.len:
+          for e in z.gval[mu]:
+            var t{.noinit.}: evalType(x.gval[mu][e])
+            t := exp(f*x.gval[mu][e])
+            z.gval[mu][e] := t*y.gval[mu][e]
+      if z.isZero:
+        z.isZero = false
+    else:
+      let laxexp = Ggauge(v.locals[0])
+      threads:
+        for mu in 0..<z.gval.len:
+          for e in z.gval[mu]:
+            var t{.noinit.}: evalType(x.gval[mu][e])
+            t := exp(f*x.gval[mu][e])
+            laxexp.gval[mu][e] := t
+            z.gval[mu][e] := t*y.gval[mu][e]
+      if z.isZero:
+        z.isZero = false
+      if laxexp.isZero:
+        laxexp.isZero = false
+      laxexp.evaluated
 
 let axexpmulyg = newGfunc(forward = axexpmulyf, backward = axexpmulyb, name = "axexpmuly")
 
-method axexpmuly*(a: Gscalar, x: Ggauge, y: Ggauge): Gvalue = Ggauge(gval: x.gval.newOneOf, inputs: @[Gvalue(a), x, y], gfunc: axexpmulyg)
+method axexpmuly*(a: Gscalar, x: Ggauge, y: Ggauge): Gvalue = Ggauge(gval: x.gval.newOneOf, inputs: @[Gvalue(a), x, y], gfunc: axexpmulyg, locals: newseq[Gvalue](2))
 
 #
 # gauge action
@@ -604,12 +845,16 @@ proc gaugeActionf(v: Gvalue) =
   let gc = Gactcoeff(v.inputs[0]).cval
   let g = Ggauge(v.inputs[1])
   let z = Gscalar(v)
-  if gc.adjplaq == 0:
-    z.getfloat = gc.gaugeAction1 g.gval
-  elif gc.rect == 0 and gc.pgm == 0:
-    z.getfloat = gc.actionA g.gval
+  if g.isZero:
+    if not z.isZero:
+      z.getfloat = 0.0
   else:
-    raiseValueError("Gauge coefficient unsupported: " & $gc)
+    if gc.adjplaq == 0:
+      z.getfloat = gc.gaugeAction1 g.gval
+    elif gc.rect == 0 and gc.pgm == 0:
+      z.getfloat = gc.actionA g.gval
+    else:
+      raiseValueError("Gauge coefficient unsupported: " & $gc)
 
 let gaugeActiong = newGfunc(forward = gaugeActionf, backward = gaugeActionb, name = "gaugeAction")
 
@@ -633,14 +878,20 @@ proc gaugeActionDerivf(v: Gvalue) =
   var gc = Gactcoeff(v.inputs[0]).cval
   let g = Ggauge(v.inputs[1])
   let z = Ggauge(v)
-  for f in gc.fields:  # gaugeActionDeriv/gaugeADeriv has a different sign than we want!
-    f = -f
-  if gc.adjplaq == 0:
-    gc.gaugeActionDeriv(g.gval, z.gval)
-  elif gc.rect == 0 and gc.pgm == 0:
-    gc.gaugeADeriv(g.gval, z.gval)
+  if g.isZero:
+    if not z.isZero:
+      z.getfloat = 0.0
   else:
-    raiseValueError("Gauge coefficient unsupported: " & $gc)
+    for f in gc.fields:  # gaugeActionDeriv/gaugeADeriv has a different sign than we want!
+      f = -f
+    if gc.adjplaq == 0:
+      gc.gaugeActionDeriv(g.gval, z.gval)
+    elif gc.rect == 0 and gc.pgm == 0:
+      gc.gaugeADeriv(g.gval, z.gval)
+    else:
+      raiseValueError("Gauge coefficient unsupported: " & $gc)
+    if z.isZero:
+      z.isZero = false
 
 let gaugeActionDerivg = newGfunc(forward = gaugeActionDerivf, backward = gaugeActionDerivb, name = "gaugeActionDeriv")
 
@@ -664,15 +915,21 @@ proc gaugeActionDeriv2f(v: Gvalue) =
   let gc = Gactcoeff(v.inputs[1]).cval
   let g = Ggauge(v.inputs[2])
   let z = Ggauge(v)
-  if gc.adjplaq == 0:
-    threads:
-      for mu in 0..<z.gval.len:
-        z.gval[mu] := 0.0
-    gc.gaugeDerivDeriv2(g.gval, b.gval, z.gval)
-  elif gc.rect == 0 and gc.pgm == 0:
-    raiseValueError("unimplemented")
+  if b.isZero or g.isZero:
+    if not z.isZero:
+      z.isZero = true
   else:
-    raiseValueError("Gauge coefficient unsupported: " & $gc)
+    if gc.adjplaq == 0:
+      threads:
+        for mu in 0..<z.gval.len:
+          z.gval[mu] := 0.0
+      gc.gaugeDerivDeriv2(g.gval, b.gval, z.gval)
+    elif gc.rect == 0 and gc.pgm == 0:
+      raiseValueError("unimplemented")
+    else:
+      raiseValueError("Gauge coefficient unsupported: " & $gc)
+    if z.isZero:
+      z.isZero = false
 
 let gaugeActionDeriv2g = newGfunc(forward = gaugeActionDeriv2f, backward = gaugeActionDeriv2b, name = "gaugeActionDeriv2")
 
diff --git a/src/experimental/graph/scalar.nim b/src/experimental/graph/scalar.nim
index fb126ec..b16eee6 100644
--- a/src/experimental/graph/scalar.nim
+++ b/src/experimental/graph/scalar.nim
@@ -29,7 +29,7 @@ method valCopy*(z: Gscalar, x: Gscalar) = z.sval = x.sval
 
 method `$`*(x: Gscalar): string = $x.sval
 
-method isTrue*(x: Gscalar): bool = x.sval != 0.0
+method isZero*(x: Gscalar): bool = x.sval == 0.0
 
 proc `getfloat=`*(x: Gvalue, y: int) =
   let xs = Gscalar(x)
@@ -58,7 +58,7 @@ method valCopy*(z: Gint, x: Gint) = z.ival = x.ival
 
 method `$`*(x: Gint): string = $x.ival
 
-method isTrue*(x: Gint): bool = x.ival != 0
+method isZero*(x: Gint): bool = x.ival == 0
 
 proc negsf(v: Gvalue) =
   let x = Gscalar(v.inputs[0])
diff --git a/src/experimental/graph/tggauge.nim b/src/experimental/graph/tggauge.nim
index 3b47922..a1c188d 100644
--- a/src/experimental/graph/tggauge.nim
+++ b/src/experimental/graph/tggauge.nim
@@ -233,6 +233,26 @@ suite "gauge fused":
     let x {.used.} = toGvalue a
     let y {.used.} = toGvalue b
 
+  test "adjmul":
+    let rf = gg.adjmul gu
+    let rg = gg.adj * gu
+    norm2(rf - rg) :< 1e-26
+    let srf = rf.norm2
+    let srg = rg.norm2
+    norm2(grad(srf, gg) - grad(srg, gg)) :< 1e-25
+    norm2(grad(srf, gu) - grad(srg, gu)) :< 1e-25
+    ckgradm2(adjmul, gg, gu, gp, gq, gm)
+
+  test "muladj":
+    let rf = gg.muladj gu
+    let rg = gg * gu.adj
+    norm2(rf - rg) :< 1e-26
+    let srf = rf.norm2
+    let srg = rg.norm2
+    norm2(grad(srf, gg) - grad(srg, gg)) :< 1e-25
+    norm2(grad(srf, gu) - grad(srg, gu)) :< 1e-25
+    ckgradm2(muladj, gg, gu, gp, gq, gm)
+
   test "contractProjTAH":
     let rf = contractProjTAH(gg, gu)
     let rg = projTAH(gg * gu.adj)