anchor.diff:      /// Original value and gradient-vector product (directional derivative) of a vector-to-scalar function `f`, at point `x`, along vector `v`
       let inline gradv' f x v =
           let i = GlobalTagger.Next
  -        Array.map2 (dualIPT i) x v |> f |> tuple
  +        Array.map2 (makeDFIPT i) x v |> f |> tuple
   
       /// Gradient-vector product (directional derivative) of a vector-to-scalar function `f`, at point `x`, along vector `v`
       let inline gradv f x v =
  
positive.code:         if (float a) <= 0. then invalidArgLog()
          match a with
          | D(a) -> D(log a)
          | DD(ap, at, ai) -> DD(log ap, at / ap, ai)
      static member Log10 (a:D) =
          if (float a) <= 0. then invalidArgLog10()
          match a with
          | D(a) -> D(log10 a)
          | DD(ap, at, ai) -> DD(log10 ap, at / (ap * log10val), ai)
      static member Exp (a:D) =
          match a with
          | D(a) -> D(exp a)
          | DD(ap, at, ai) -> let expa = exp ap in DD(expa, at * expa, ai)
      static member Sin (a:D) =
          match a with
          | D(a) -> D(sin a)
          | DD(ap, at, ai) -> DD(sin ap, at * cos ap, ai)
      static member Cos (a:D) =
          match a with
          | D(a) -> D(cos a)
          | DD(ap, at, ai) -> DD(cos ap, -at * sin ap, ai)
      static member Tan (a:D) =
          match a with
          | D(a) -> 
              if cos a = 0. then invalidArgTan()
              D(tan a)
          | DD(ap, at, ai) ->
              let cosa = cos ap
              if (float cosa) = 0. then invalidArgTan()
              DD(tan ap, at / (cosa * cosa), ai)
      static member (~-) (a:D) =
          match a with
          | D(a) -> D(-a)
          | DD(ap, at, ai) -> DD(-ap, -at, ai)
      static member Sqrt (a:D) =
          if (float a) <= 0. then invalidArgSqrt()
          match a with
          | D(a) -> D(sqrt a)
          | DD(ap, at, ai) -> let sqrta = sqrt ap in DD(sqrta, at / (2. * sqrta), ai)
      static member Sinh (a:D) =
          match a with
          | D(a) -> D(sinh a)
          | DD(ap, at, ai) -> DD(sinh ap, at * cosh ap, ai)
      static member Cosh (a:D) =
          match a with
          | D(a) -> D(cosh a)
          | DD(ap, at, ai) -> DD(cosh ap, at * sinh ap, ai)
      static member Tanh (a:D) =
          match a with
          | D(a) -> D(tanh a)
          | DD(ap, at, ai) -> let cosha = cosh ap in DD(tanh ap, at / (cosha * cosha), ai)
      static member Asin (a:D) =
          if abs (float a) >= 1. then invalidArgAsin()
          match a with
          | D(a) -> D(asin a)
          | DD(ap, at, ai) -> DD(asin ap, at / sqrt (1. - ap * ap), ai)
      static member Acos (a:D) =
          if abs (float a) >= 1. then invalidArgAcos()
          match a with
          | D(a) -> D(acos a)
          | DD(ap, at, ai) -> DD(acos ap, -at / sqrt (1. - ap * ap), ai)
      static member Atan (a:D) =
          match a with
          | D(a) -> D(atan a)
          | DD(ap, at, ai) -> DD(atan ap, at / (1. + ap * ap), ai)
      static member Abs (a:D) =
          if (float a) = 0. then invalidArgAbs()
          match a with
          | D(a) -> D(abs a)
          | DD(ap, at, ai) -> DD(abs ap, at * float (sign (float ap)), ai)
      static member Floor (a:D) =
          if isInteger (float a) then invalidArgFloor()
          match a with
          | D(a) -> D(floor a)
          | DD(ap, _, ai) -> DD(floor ap, D 0., ai)
      static member Ceiling (a:D) =
          if isInteger (float a) then invalidArgCeil()
          match a with
          | D(a) -> D(ceil a)
          | DD(ap, _, ai) -> DD(ceil ap, D 0., ai)
      static member Round (a:D) =
          if isHalfway (float a) then invalidArgRound()
          match a with
          | D(a) -> D(round a)
          | DD(ap, _, ai) -> DD(round ap, D 0., ai)
  
  
  /// Tagger for generating incremental integers
  
negative.code:     override d.GetHashCode() =
          match d with
          | D(a) -> hash [| a |]
          | DF(ap, at, ai) -> hash [|ap; at; ai|]
      // D - D binary operations
      static member (+) (a:D, b:D) =
          match a, b with
          | D(a),           D(b)           -> D(a + b)
          | D(a),           DF(bp, bt, bi) -> DF(a + bp, bt, bi)
          | DF(ap, at, ai), D(b)           -> DF(ap + b, at, ai)
          | DF(_ , _ , ai), DF(bp, bt, bi) when ai < bi -> DF(a + bp, bt, bi)
          | DF(ap, at, ai), DF(bp, bt, bi) when ai = bi -> DF(ap + bp, at + bt, ai)
          | DF(ap, at, ai), DF(_ , _ , bi) when ai > bi -> DF(ap + b, at, ai)
      static member (-) (a:D, b:D) =
          match a, b with
          | D(a),           D(b)           -> D(a - b)
          | D(a),           DF(bp, bt, bi) -> DF(a - bp, -bt, bi)
          | DF(ap, at, ai), D(b)           -> DF(ap - b, at, ai)
          | DF(_ , _ , ai), DF(bp, bt, bi) when ai < bi -> DF(a - bp, -bt, bi)
          | DF(ap, at, ai), DF(bp, bt, bi) when ai = bi -> DF(ap - bp, at - bt, ai)
          | DF(ap, at, ai), DF(_ , _ , bi) when ai > bi -> DF(ap - b, at, ai)
      static member (*) (a:D, b:D) =
          match a, b with
          | D(a),           D(b)           -> D(a * b)
          | D(a),           DF(bp, bt, bi) -> DF(a * bp, a * bt, bi)
          | DF(ap, at, ai), D(b)           -> DF(ap * b, at * b, ai)
          | DF(_ , _ , ai), DF(bp, bt, bi) when ai < bi -> DF(a * bp, bt * a, bi)
          | DF(ap, at, ai), DF(bp, bt, bi) when ai = bi -> DF(ap * bp, at * bp + bt * ap, ai)
          | DF(ap, at, ai), DF(_ , _ , bi) when ai > bi -> DF(ap * b, at * b, ai)
      static member (/) (a:D, b:D) =
          match a, b with
          | D(a),           D(b)           -> D(a / b)
          | D(a),           DF(bp, bt, bi) -> DF(a / bp, -(bt * a) / (bp * bp), bi)
          | DF(ap, at, ai), D(b)           -> DF(ap / b, at / b, ai)
          | DF(_ , _ , ai), DF(bp, bt, bi) when ai < bi -> DF(a / bp, -(bt * a) / (bp * bp), bi)
          | DF(ap, at, ai), DF(bp, bt, bi) when ai = bi -> DF(ap / bp, (at * bp - bt * ap) / (bp * bp), ai)
          | DF(ap, at, ai), DF(_ , _ , bi) when ai > bi -> DF(ap / b, at / b, ai)
      static member Pow (a:D, b:D) =
          match a, b with
          | D(a),           D(b)           -> D(a ** b)
          | D(a),           DF(bp, bt, bi) -> let apb = D.Pow(a, bp) in DF(apb, apb * (log a) * bt, bi)
          | DF(ap, at, ai), D(b)           -> DF(ap ** b, b * (ap ** (b - 1.)) * at, ai)
          | DF(_ , _ , ai), DF(bp, bt, bi) when ai < bi -> let apb = a ** bp in DF(apb, apb * (log a) * bt, bi)
          | DF(ap, at, ai), DF(bp, bt, bi) when ai = bi -> let apb = ap ** bp in DF(apb, apb * ((bp * at / ap) + ((log ap) * bt)), ai)
          | DF(ap, at, ai), DF(_ , _ , bi) when ai > bi -> let apb = ap ** b in DF(apb, apb * (b * at / ap), ai)
      static member Atan2 (a:D, b:D) =
          match a, b with
          | D(a),           D(b)           -> D(atan2 a b)
          | D(a),           DF(bp, bt, bi) -> DF(D.Atan2(a, bp), -(a * bt) / (a * a + bp * bp), bi)
          | DF(ap, at, ai), D(b)           -> DF(D.Atan2(ap, b), (b * at) / (ap * ap + b * b), ai)
          | DF(_ , _ , ai), DF(bp, bt, bi) when ai < bi -> DF(atan2 a bp, -(a * bt) / (a * a + bp * bp), bi)
          | DF(ap, at, ai), DF(bp, bt, bi) when ai = bi -> DF(atan2 ap bp, (at * bp - ap * bt) / (ap * ap + bp * bp), ai)
          | DF(ap, at, ai), DF(_ , _ , bi) when ai > bi -> DF(atan2 ap b, (at * b) / (ap * ap + b * b), ai)
      // D - float binary operations
      static member (+) (a:D, b:float) =
          match a with
          | D(a) -> D(a + b)
          | DF(ap, at, ai) -> DF(ap + b, at, ai)
      static member (-) (a:D, b:float) =
          match a with
          | D(a) -> D(a - b)
          | DF(ap, at, ai) -> DF(ap - b, at, ai)
      static member (*) (a:D, b:float) =
          match a with
          | D(a) -> D(a * b)
          | DF(ap, at, ai) -> DF(ap * b, at * b, ai)
      static member (/) (a:D, b:float) =
          match a with
          | D(a) -> D(a / b)
          | DF(ap, at, ai) -> DF(ap / b, at / b, ai)
      static member Pow (a:D, b:float) =
          match a with
          | D(a) -> D(a ** b)
          | DF(ap, at, ai) -> DF(ap ** b, b * (ap ** (b - 1.)) * at, ai)
      static member Atan2 (a:D, b:float) =
          match a with
          | D(a) -> D(atan2 a b)
          | DF(ap, at, ai) -> DF(D.Atan2(ap, b), (b * at) / (b * b + ap * ap), ai)
      // float - D binary operations
      static member (+) (a:float, b:D) =
          match b with
          | D(b) -> D(a + b)
          | DF(bp, bt, bi) -> DF(a + bp, bt, bi)
      static member (-) (a:float, b:D) =
          match b with
          | D(b) -> D(a - b)
          | DF(bp, bt, bi) -> DF(a - bp, -bt, bi)
      static member (*) (a:float, b:D) =
          match b with
          | D(b) -> D(a * b)
          | DF(bp, bt, bi) -> DF(a * bp, a * bt, bi)
      static member (/) (a:float, b:D) =
          match b with
          | D(b) -> D(a / b)
          | DF(bp, bt, bi) -> DF(a / bp, -a * bt / (bp * bp), bi)
      static member Pow (a:float, b:D) =
          match b with
          | D(b) -> D(a ** b)
          | DF(bp, bt, bi) -> let apb = D.Pow(a, bp) in DF(apb, apb * (log a) * bt , bi)
      static member Atan2 (a:float, b:D) =
          match b with
          | D(b) -> D(atan2 a b)
          | DF(bp, bt, bi) -> DF(D.Atan2(a, bp), -(a * bt) / (a * a + bp * bp), bi)
      // D - int binary operations
      static member (+) (a:D, b:int) = a + float b
      static member (-) (a:D, b:int) = a - float b
  