Index: /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs
===================================================================
--- /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs	(revision 16692)
+++ /branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs	(revision 16693)
@@ -4,5 +4,4 @@
 using HeuristicLab.Common;
 using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
-using HEAL.Attic;
 
 namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
@@ -117,5 +116,5 @@
 
 
-  public class VectorEvaluator : Interpreter<DoubleVector> {
+  public sealed class VectorEvaluator : Interpreter<DoubleVector> {
     private const int BATCHSIZE = 128;
     [ThreadStatic]
@@ -192,5 +191,5 @@
   }
 
-  public class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<DoubleVector>> {
+  public sealed class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<DoubleVector>> {
     private const int BATCHSIZE = 128;
     [ThreadStatic]
@@ -318,5 +317,5 @@
 
 
-  public class IntervalEvaluator : Interpreter<AlgebraicInterval> {
+  public sealed class IntervalEvaluator : Interpreter<AlgebraicInterval> {
     [ThreadStatic]
     private Dictionary<string, Interval> intervals;
@@ -352,6 +351,6 @@
       instruction.dblVal = constant.Value;
       instruction.value = new AlgebraicInterval(
-        new MultivariateDual<Double>(constant.Value, constant, 1.0),
-        new MultivariateDual<Double>(constant.Value, constant, 1.0) // use node as key
+        new MultivariateDual<AlgebraicDouble>(constant.Value, constant, 1.0),
+        new MultivariateDual<AlgebraicDouble>(constant.Value, constant, 1.0) // use node as key
         );
     }
@@ -360,6 +359,6 @@
       instruction.dblVal = variable.Weight;
       instruction.value = new AlgebraicInterval(
-        low: new MultivariateDual<Double>(intervals[variable.VariableName].LowerBound, variable, intervals[variable.VariableName].LowerBound),  // bounds change by variable value d/dc (c I(var)) = I(var)
-        high: new MultivariateDual<Double>(intervals[variable.VariableName].UpperBound, variable, intervals[variable.VariableName].UpperBound)
+        low: new MultivariateDual<AlgebraicDouble>(intervals[variable.VariableName].LowerBound, variable, intervals[variable.VariableName].LowerBound),  // bounds change by variable value d/dc (c I(var)) = I(var)
+        high: new MultivariateDual<AlgebraicDouble>(intervals[variable.VariableName].UpperBound, variable, intervals[variable.VariableName].UpperBound)
         );
     }
@@ -371,4 +370,6 @@
 
   public interface IAlgebraicType<T> {
+    T Zero { get; }
+
     T AssignAbs(T a); // set this to assign abs(a)
     T Assign(T a); // assign this to same value as a (copy!)
@@ -414,29 +415,29 @@
 
   // algebraic type wrapper for a double value
-  public class Double : IAlgebraicType<Double> {
-    public static implicit operator Double(double value) { return new Double(value); }
-    public static implicit operator double(Double value) { return value.Value; }
+  public sealed class AlgebraicDouble : IAlgebraicType<AlgebraicDouble> {
+    public static implicit operator AlgebraicDouble(double value) { return new AlgebraicDouble(value); }
+    public static implicit operator double(AlgebraicDouble value) { return value.Value; }
     public double Value;
-    public Double() { }
-    public Double(double value) { this.Value = value; }
-    public Double Add(Double a) { Value += a.Value; return this; }
-    public Double Assign(Double a) { Value = a.Value; return this; }
-
-    public Double AssignAbs(Double a) { Value = Math.Abs(a.Value); return this; }
-    public Double AssignCos(Double a) { Value = Math.Cos(a.Value); return this; }
-    public Double AssignExp(Double a) { Value = Math.Exp(a.Value); return this; }
-    public Double AssignIntPower(Double a, int p) { Value = Math.Pow(a.Value, p); return this; }
-    public Double AssignIntRoot(Double a, int r) { Value = Math.Pow(a.Value, 1.0 / r); return this; }
-    public Double AssignInv(Double a) { Value = 1.0 / a.Value; return this; }
-    public Double AssignLog(Double a) { Value = Math.Log(a.Value); return this; }
-    public Double AssignNeg(Double a) { Value = -a.Value; return this; }
-    public Double AssignSin(Double a) { Value = Math.Sin(a.Value); return this; }
-    public Double AssignSgn(Double a) { Value = Math.Sign(a.Value); return this; }
-    public Double Clone() { return new Double(Value); }
-    public Double Div(Double a) { Value /= a.Value; return this; }
-    public Double Mul(Double a) { Value *= a.Value; return this; }
-    public Double Scale(double s) { Value *= s; return this; }
-    public Double Sub(Double a) { Value -= a.Value; return this; }
-
+
+    public AlgebraicDouble Zero => new AlgebraicDouble(0.0);
+    public AlgebraicDouble() { }
+    public AlgebraicDouble(double value) { this.Value = value; }
+    public AlgebraicDouble Assign(AlgebraicDouble a) { Value = a.Value; return this; }
+    public AlgebraicDouble Add(AlgebraicDouble a) { Value += a.Value; return this; }
+    public AlgebraicDouble Sub(AlgebraicDouble a) { Value -= a.Value; return this; }
+    public AlgebraicDouble Mul(AlgebraicDouble a) { Value *= a.Value; return this; }
+    public AlgebraicDouble Div(AlgebraicDouble a) { Value /= a.Value; return this; }
+    public AlgebraicDouble Scale(double s) { Value *= s; return this; }
+    public AlgebraicDouble AssignInv(AlgebraicDouble a) { Value = 1.0 / a.Value; return this; }
+    public AlgebraicDouble AssignNeg(AlgebraicDouble a) { Value = -a.Value; return this; }
+    public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; }
+    public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; }
+    public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = Math.Log(a.Value); return this; }
+    public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; }
+    public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; }
+    public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = Math.Pow(a.Value, 1.0 / r); return this; }
+    public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; }
+    public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = Math.Sign(a.Value); return this; }
+    public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); }
   }
 
@@ -445,5 +446,8 @@
     private double[] arr;
 
+
     public double this[int idx] { get { return arr[idx]; } set { arr[idx] = value; } }
+
+    public int Length => arr.Length;
 
     public DoubleVector(int length) {
@@ -461,41 +465,22 @@
     }
 
-    public DoubleVector(int length, double constantValue) : this(length) {
-      for (int i = 0; i < length; ++i) arr[i] = constantValue;
-    }
-
-    public void CopyTo(double[] dest, int idx, int length) {
-      Array.Copy(arr, 0, dest, idx, length);
-    }
-
-    public void CopyRowTo(double[,] dest, int row) {
-      for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];
-    }
-    internal void CopyColumnTo(double[,] dest, int column, int row, int len) {
-      for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];
-    }
-
-    public void AssignConstant(double constantValue) {
-      for (int i = 0; i < arr.Length; ++i) {
-        arr[i] = constantValue;
-      }
-    }
-
-    public DoubleVector Assign(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
-    public DoubleVector AssignCos(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
-    public DoubleVector Div(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
-    public DoubleVector AssignExp(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
-    public DoubleVector AssignIntPower(DoubleVector a, int p) { throw new NotImplementedException(); }
-    public DoubleVector AssignIntRoot(DoubleVector a, int r) { throw new NotImplementedException(); }
-    public DoubleVector AssignInv(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
-    public DoubleVector AssignLog(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
+    public DoubleVector Zero => new DoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static)
+    public DoubleVector Assign(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
     public DoubleVector Add(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }
+    public DoubleVector Sub(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
     public DoubleVector Mul(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }
-    public DoubleVector AssignNeg(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
+    public DoubleVector Div(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
+    public DoubleVector AssignNeg(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
+    public DoubleVector AssignInv(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
     public DoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }
-    public DoubleVector AssignSin(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
-    public DoubleVector Sub(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
-    public DoubleVector AssignAbs(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; }
-    public DoubleVector AssignSgn(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; }
+    public DoubleVector AssignLog(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
+    public DoubleVector AssignSin(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
+    public DoubleVector AssignExp(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
+    public DoubleVector AssignCos(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
+    public DoubleVector AssignIntPower(DoubleVector a, int p) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Pow(a.arr[i], p); } return this; }
+    public DoubleVector AssignIntRoot(DoubleVector a, int r) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Pow(a.arr[i], 1.0 / r); } return this; }
+    public DoubleVector AssignAbs(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; }
+    public DoubleVector AssignSgn(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; }
+
     public DoubleVector Clone() {
       var v = new DoubleVector(this.arr.Length);
@@ -504,12 +489,29 @@
     }
 
+    public void AssignConstant(double constantValue) {
+      for (int i = 0; i < arr.Length; ++i) {
+        arr[i] = constantValue;
+      }
+    }
+
+    public void CopyTo(double[] dest, int idx, int length) {
+      Array.Copy(arr, 0, dest, idx, length);
+    }
+
     public void CopyFrom(double[] data, int rowIndex) {
       Array.Copy(data, rowIndex, arr, 0, Math.Min(arr.Length, data.Length - rowIndex));
     }
 
+    public void CopyRowTo(double[,] dest, int row) {
+      for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];
+    }
+
+    internal void CopyColumnTo(double[,] dest, int column, int row, int len) {
+      for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];
+    }
   }
 
   // vectors of algebraic types
-  public class Vector<T> : IAlgebraicType<Vector<T>> where T : IAlgebraicType<T> {
+  public sealed class Vector<T> : IAlgebraicType<Vector<T>> where T : IAlgebraicType<T> {
     private T[] elems;
 
@@ -517,4 +519,5 @@
 
     public int Length => elems.Length;
+
 
     private Vector() { }
@@ -527,5 +530,5 @@
     /// 
     /// </summary>
-    /// <param name="elems">The array is copied</param>
+    /// <param name="elems">The array is copied (element-wise clone)</param>
     public Vector(T[] elems) {
       this.elems = new T[elems.Length];
@@ -562,19 +565,20 @@
     }
 
+    public Vector<T> Zero => new Vector<T>(Length);
+    public Vector<T> Assign(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
     public Vector<T> Add(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }
-    public Vector<T> Assign(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
+    public Vector<T> Sub(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }
+    public Vector<T> Mul(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }
+    public Vector<T> Div(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }
+    public Vector<T> AssignNeg(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }
+    public Vector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }
+    public Vector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
+    public Vector<T> AssignInv(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }
+    public Vector<T> AssignLog(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }
+    public Vector<T> AssignExp(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }
+    public Vector<T> AssignSin(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }
     public Vector<T> AssignCos(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignCos(a.elems[i]); } return this; }
-    public Vector<T> AssignExp(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }
     public Vector<T> AssignIntPower(Vector<T> a, int p) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntPower(a.elems[i], p); } return this; }
     public Vector<T> AssignIntRoot(Vector<T> a, int r) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntRoot(a.elems[i], r); } return this; }
-    public Vector<T> AssignInv(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }
-    public Vector<T> AssignLog(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }
-    public Vector<T> AssignNeg(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }
-    public Vector<T> AssignSin(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }
-    public Vector<T> Div(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }
-    public Vector<T> Mul(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }
-    public Vector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }
-    public Vector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
-    public Vector<T> Sub(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }
     public Vector<T> AssignAbs(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; }
     public Vector<T> AssignSgn(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSgn(a.elems[i]); } return this; }
@@ -587,9 +591,9 @@
   /// <typeparam name="K">Key type</typeparam>
   /// <typeparam name="T">Element type</typeparam>
-  public class SparseVector<K, T> : IAlgebraicType<SparseVector<K, T>> where T : IAlgebraicType<T> {
+  public sealed class SparseVector<K, T> : IAlgebraicType<SparseVector<K, T>> where T : IAlgebraicType<T> {
 
     private Dictionary<K, T> elems;
-
     public IReadOnlyDictionary<K, T> Elements => elems;
+
 
     public SparseVector(SparseVector<K, T> original) {
@@ -597,9 +601,14 @@
     }
 
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="keys"></param>
+    /// <param name="values">values are cloned</param>
     public SparseVector(K[] keys, T[] values) {
       if (keys.Length != values.Length) throw new ArgumentException("lengths of keys and values doesn't match in SparseVector");
       elems = new Dictionary<K, T>(keys.Length);
       for (int i = 0; i < keys.Length; ++i) {
-        elems.Add(keys[i], values[i]);
+        elems.Add(keys[i], values[i].Clone());
       }
     }
@@ -609,119 +618,56 @@
     }
 
-    public SparseVector<K, T> Add(SparseVector<K, T> a) {
+
+
+    private void AssignTransformed(SparseVector<K, T> a, Func<T, T, T> mapAssign) {
       foreach (var kvp in a.elems) {
-        if (elems.TryGetValue(kvp.Key, out var value)) {
-          value.Add(kvp.Value);
-        } else {
-          elems.Add(kvp.Key, kvp.Value.Clone());
+        if (elems.TryGetValue(kvp.Key, out T value))
+          mapAssign(kvp.Value, value);
+        else {
+          var newValue = kvp.Value.Zero;
+          elems.Add(kvp.Key, newValue);
+          mapAssign(kvp.Value, newValue);
         }
       }
-      return this;
-    }
-
-    public SparseVector<K, T> Scale(T s) {
-      foreach (var kvp in elems) {
-        kvp.Value.Mul(s);
-      }
-      return this;
-    }
-
-    public SparseVector<K, T> Scale(double s) {
-      foreach (var kvp in elems) {
-        kvp.Value.Scale(s);
-      }
-      return this;
-    }
-
-
-    public SparseVector<K, T> Assign(SparseVector<K, T> a) {
-      elems.Clear();
-      elems = a.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone());
-      return this;
-    }
-
-    public SparseVector<K, T> AssignCos(SparseVector<K, T> a) {
-      throw new NotImplementedException();
-    }
-
-    public SparseVector<K, T> AssignExp(SparseVector<K, T> a) {
-      throw new NotImplementedException();
-    }
-
-    public SparseVector<K, T> AssignIntPower(SparseVector<K, T> a, int p) {
-      throw new NotImplementedException();
-    }
-
-    public SparseVector<K, T> AssignIntRoot(SparseVector<K, T> a, int r) {
-      throw new NotImplementedException();
-    }
-
-    public SparseVector<K, T> AssignInv(SparseVector<K, T> a) {
-      throw new NotImplementedException();
-    }
-
-    public SparseVector<K, T> AssignLog(SparseVector<K, T> a) {
-      throw new NotImplementedException();
-    }
-
-    public SparseVector<K, T> AssignNeg(SparseVector<K, T> a) {
-      throw new NotImplementedException();
-    }
-
-    public SparseVector<K, T> AssignSin(SparseVector<K, T> a) {
-      throw new NotImplementedException();
-    }
+    }
+
+    public SparseVector<K, T> Zero => new SparseVector<K, T>();
+
+    public SparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; }
+    public SparseVector<K, T> Scale(double s) { foreach (var kvp in elems) { kvp.Value.Scale(s); } return this; }
+
+    public SparseVector<K, T> Assign(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Assign(src)); return this; }
+    public SparseVector<K, T> Add(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Add(src)); return this; }
+    public SparseVector<K, T> Mul(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Mul(src)); return this; }
+    public SparseVector<K, T> Sub(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Sub(src)); return this; }
+    public SparseVector<K, T> Div(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Div(src)); return this; }
+    public SparseVector<K, T> AssignInv(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignInv(src)); return this; }
+    public SparseVector<K, T> AssignNeg(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignNeg(src)); return this; }
+    public SparseVector<K, T> AssignLog(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignLog(src)); return this; }
+    public SparseVector<K, T> AssignExp(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignExp(src)); return this; }
+    public SparseVector<K, T> AssignIntPower(SparseVector<K, T> a, int p) { AssignTransformed(a, (src, dest) => dest.AssignIntPower(src, p)); return this; }
+    public SparseVector<K, T> AssignIntRoot(SparseVector<K, T> a, int r) { AssignTransformed(a, (src, dest) => dest.AssignIntRoot(src, r)); return this; }
+    public SparseVector<K, T> AssignSin(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignSin(src)); return this; }
+    public SparseVector<K, T> AssignCos(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignCos(src)); return this; }
+    public SparseVector<K, T> AssignAbs(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignAbs(src)); return this; }
+    public SparseVector<K, T> AssignSgn(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignSgn(src)); return this; }
 
     public SparseVector<K, T> Clone() {
       return new SparseVector<K, T>(this);
     }
-
-    public SparseVector<K, T> Div(SparseVector<K, T> a) {
-      throw new NotImplementedException();
-    }
-
-    public SparseVector<K, T> Mul(SparseVector<K, T> a) {
-      throw new NotImplementedException();
-    }
-
-    public SparseVector<K, T> Sub(SparseVector<K, T> a) {
-      foreach (var kvp in a.elems) {
-        if (elems.TryGetValue(kvp.Key, out var value)) {
-          value.Sub(kvp.Value);
-        } else {
-          elems.Add(kvp.Key, kvp.Value.Clone().Neg());
-        }
-      }
-      return this;
-    }
-
-    public SparseVector<K, T> AssignAbs(SparseVector<K, T> a) {
-      elems.Clear();
-      foreach (var kvp in a.elems) {
-        elems.Add(kvp.Key, kvp.Value.Clone().Abs());
-      }
-      return this;
-    }
-
-    public SparseVector<K, T> AssignSgn(SparseVector<K, T> a) {
-      elems.Clear();
-      foreach (var kvp in a.elems) {
-        elems.Add(kvp.Key, kvp.Value.Clone().Sgn());
-      }
-      return this;
-    }
   }
 
 
   public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
-    private MultivariateDual<Double> low;
-    private MultivariateDual<Double> high;
-
-    public MultivariateDual<Double> LowerBound => low.Clone();
-    public MultivariateDual<Double> UpperBound => high.Clone();
+    private MultivariateDual<AlgebraicDouble> low;
+    private MultivariateDual<AlgebraicDouble> high;
+
+    public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();
+    public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();
+
 
     public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }
 
-    public AlgebraicInterval(MultivariateDual<Double> low, MultivariateDual<Double> high) {
+    public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {
       this.low = low.Clone();
       this.high = high.Clone();
@@ -729,62 +675,12 @@
 
     public AlgebraicInterval(double low, double high) {
-      this.low = new MultivariateDual<Double>(new Double(low));
-      this.high = new MultivariateDual<Double>(new Double(high));
-    }
-
+      this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));
+      this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));
+    }
+
+    public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);
     public AlgebraicInterval Add(AlgebraicInterval a) {
       low.Add(a.low);
       high.Add(a.high);
-      return this;
-    }
-
-    public AlgebraicInterval Assign(AlgebraicInterval a) {
-      low = a.low;
-      high = a.high;
-      return this;
-    }
-
-    public AlgebraicInterval AssignCos(AlgebraicInterval a) {
-      throw new NotImplementedException();
-    }
-
-    public AlgebraicInterval Div(AlgebraicInterval a) {
-      if (a.Contains(0.0)) {
-        if (a.low.Value.Value.IsAlmost(0.0) && a.high.Value.Value.IsAlmost(0.0)) {
-          low = new MultivariateDual<Double>(double.NegativeInfinity);
-          high = new MultivariateDual<Double>(double.PositiveInfinity);
-        } else if (a.low.Value.Value.IsAlmost(0.0))
-          Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<Double>(double.PositiveInfinity)));
-        else
-          Mul(new AlgebraicInterval(new MultivariateDual<Double>(double.NegativeInfinity), a.low.Clone().Inv()));
-      } else {
-        Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
-      }
-      return this;
-    }
-
-    public AlgebraicInterval AssignExp(AlgebraicInterval a) {
-      low.AssignExp(a.low);
-      high.AssignExp(a.high);
-      return this;
-    }
-
-    public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
-      throw new NotImplementedException();
-    }
-
-    public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
-      throw new NotImplementedException();
-    }
-
-    public AlgebraicInterval AssignInv(AlgebraicInterval a) {
-      low = new MultivariateDual<Double>(1.0);
-      high = new MultivariateDual<Double>(1.0);
-      return Div(a);
-    }
-
-    public AlgebraicInterval AssignLog(AlgebraicInterval a) {
-      low.AssignLog(a.low);
-      high.AssignLog(a.high);
       return this;
     }
@@ -801,6 +697,107 @@
     }
 
+    public AlgebraicInterval Assign(AlgebraicInterval a) {
+      low = a.low;
+      high = a.high;
+      return this;
+    }
+
+    public AlgebraicInterval AssignCos(AlgebraicInterval a) {
+      return AssignSin(a.Clone().Sub(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));
+    }
+
+    public AlgebraicInterval Div(AlgebraicInterval a) {
+      if (a.Contains(0.0)) {
+        if (a.low.Value.Value.IsAlmost(0.0) && a.high.Value.Value.IsAlmost(0.0)) {
+          low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
+          high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
+        } else if (a.low.Value.Value.IsAlmost(0.0))
+          Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));
+        else
+          Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));
+      } else {
+        Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
+      }
+      return this;
+    }
+
+    public AlgebraicInterval AssignExp(AlgebraicInterval a) {
+      low.AssignExp(a.low);
+      high.AssignExp(a.high);
+      return this;
+    }
+
+    public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
+      if (p == 0) {
+        // => 1
+        low = new MultivariateDual<AlgebraicDouble>(1.0);
+        high = new MultivariateDual<AlgebraicDouble>(1.0);
+        return this;
+      }
+      if (p == 1) return this;
+
+      if (p < 0) {  // x^-3 == 1/(x^3)
+        AssignIntPower(a, -p);
+        return AssignInv(this);
+      } else {
+        // p is even => interval must be positive
+        if (p % 2 == 0) {
+          if (a.Contains(0.0)) {
+            low = new MultivariateDual<AlgebraicDouble>(0.0);
+            high = Algebraic.Max(low.Clone().IntPower(p), high.Clone().IntPower(p));
+          } else {
+            var lowPower = low.Clone().IntPower(p);
+            var highPower = high.Clone().IntPower(p);
+            low = Algebraic.Min(lowPower, highPower);
+            high = Algebraic.Max(lowPower, highPower);
+          }
+        } else {
+          // p is uneven
+          var lowPower = low.Clone().IntPower(p);
+          var highPower = high.Clone().IntPower(p);
+          low = Algebraic.Min(lowPower, highPower);
+          high = Algebraic.Max(lowPower, highPower);
+        }
+        return this;
+      }
+    }
+
+    public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
+      if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }
+      if (r == 1) return this;
+      if (r < 0) {
+        // x^ (-1/2) = 1 / (x^(1/2))
+        AssignIntRoot(a, -r);
+        return AssignInv(this);
+      } else {
+        // root only defined for positive arguments
+        if (a.LowerBound.Value.Value < 0) {
+          low = new MultivariateDual<AlgebraicDouble>(double.NaN);
+          high = new MultivariateDual<AlgebraicDouble>(double.NaN);
+          return this;
+        } else {
+          low.AssignIntRoot(a.low, r);
+          high.AssignIntRoot(a.high, r);
+          return this;
+        }
+      }
+    }
+
+    public AlgebraicInterval AssignInv(AlgebraicInterval a) {
+      low = new MultivariateDual<AlgebraicDouble>(1.0);
+      high = new MultivariateDual<AlgebraicDouble>(1.0);
+      return Div(a);
+    }
+
+    public AlgebraicInterval AssignLog(AlgebraicInterval a) {
+      low.AssignLog(a.low);
+      high.AssignLog(a.high);
+      return this;
+    }
+
     public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
-      throw new NotImplementedException();
+      low.AssignNeg(a.high);
+      high.AssignNeg(a.low);
+      return this;
     }
 
@@ -817,5 +814,34 @@
 
     public AlgebraicInterval AssignSin(AlgebraicInterval a) {
-      throw new NotImplementedException();
+      if (Math.Abs(a.UpperBound.Value.Value - a.LowerBound.Value.Value) >= Math.PI * 2) {
+        low = new MultivariateDual<AlgebraicDouble>(-1.0);
+        high = new MultivariateDual<AlgebraicDouble>(1.0);
+      }
+
+      //divide the interval by PI/2 so that the optima lie at x element of N (0,1,2,3,4,...)
+      double Pihalf = Math.PI / 2;
+      var scaled = this.Clone().Scale(1.0 / Pihalf);
+      //move to positive scale
+      if (scaled.LowerBound.Value.Value < 0) {
+        int periodsToMove = Math.Abs((int)scaled.LowerBound.Value.Value / 4) + 1;
+        scaled.Add(new AlgebraicInterval(periodsToMove * 4, periodsToMove * 4));
+      }
+
+      double scaledLowerBound = scaled.LowerBound.Value.Value % 4.0;
+      double scaledUpperBound = scaled.UpperBound.Value.Value % 4.0;
+      if (scaledUpperBound < scaledLowerBound) scaledUpperBound += 4.0;
+      List<double> sinValues = new List<double>();
+      sinValues.Add(Math.Sin(scaledLowerBound * Pihalf));
+      sinValues.Add(Math.Sin(scaledUpperBound * Pihalf));
+
+      int startValue = (int)Math.Ceiling(scaledLowerBound);
+      while (startValue < scaledUpperBound) {
+        sinValues.Add(Math.Sin(startValue * Pihalf));
+        startValue += 1;
+      }
+
+      low = new MultivariateDual<AlgebraicDouble>(sinValues.Min());
+      high = new MultivariateDual<AlgebraicDouble>(sinValues.Max());
+      return this;
     }
 
@@ -840,5 +866,5 @@
         var abshigh = a.high.Clone().Abs();
         a.high.Assign(Algebraic.Max(abslow, abshigh));
-        a.low.Assign(new MultivariateDual<Double>(0.0)); // lost gradient for lower bound
+        a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound
       } else {
         var abslow = a.low.Clone().Abs();
@@ -851,5 +877,7 @@
 
     public AlgebraicInterval AssignSgn(AlgebraicInterval a) {
-      throw new NotImplementedException();
+      low.AssignSgn(a.low);
+      high.AssignSgn(a.high);
+      return this;
     }
   }
@@ -864,4 +892,5 @@
     public V Derivative => dv;
 
+
     public Dual(V v, V dv) {
       this.v = v;
@@ -873,4 +902,6 @@
     }
 
+    public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero);
+
     public Dual<V> Add(Dual<V> a) {
       v.Add(a.v);
@@ -892,5 +923,6 @@
 
     public Dual<V> Div(Dual<V> a) {
-      throw new NotImplementedException();
+      Mul(a.Inv());
+      return this;
     }
 
@@ -902,13 +934,19 @@
 
     public Dual<V> AssignIntPower(Dual<V> a, int p) {
-      throw new NotImplementedException();
+      v.AssignIntPower(a.v, p);
+      dv.Assign(a.dv).Scale(p).Mul(a.v.IntPower(p - 1));
+      return this;
     }
 
     public Dual<V> AssignIntRoot(Dual<V> a, int r) {
-      throw new NotImplementedException();
+      v.AssignIntRoot(a.v, r);
+      dv.Assign(a.dv).Scale(1.0 / r).Mul(a.v.IntRoot(r - 1));
+      return this;
     }
 
     public Dual<V> AssignInv(Dual<V> a) {
-      throw new NotImplementedException();
+      v.AssignInv(a.v);
+      dv.AssignNeg(a.dv).Mul(v).Mul(v); // (1/f(x))' = - f(x)' / f(x)^2
+      return this;
     }
 
@@ -935,5 +973,7 @@
 
     public Dual<V> AssignNeg(Dual<V> a) {
-      throw new NotImplementedException();
+      v.AssignNeg(a.v);
+      dv.AssignNeg(a.dv);
+      return this;
     }
 
@@ -945,9 +985,13 @@
 
     public Dual<V> AssignSin(Dual<V> a) {
-      throw new NotImplementedException();
+      v.AssignSin(a.v);
+      dv.Assign(a.dv).Mul(a.v.Clone().Cos());
+      return this;
     }
 
     public Dual<V> Sub(Dual<V> a) {
-      throw new NotImplementedException();
+      v.Sub(a.v);
+      dv.Sub(a.dv);
+      return this;
     }
 
@@ -1011,7 +1055,21 @@
     }
 
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="v">not cloned</param>
+    /// <param name="gradient">not cloned</param>
+    internal MultivariateDual(V v, SparseVector<object, V> gradient) {
+      this.v = v;
+      this.dv = gradient;
+    }
+
     public MultivariateDual<V> Clone() {
       return new MultivariateDual<V>(this);
     }
+
+
+    public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero);
+
 
     public MultivariateDual<V> Add(MultivariateDual<V> a) {
@@ -1097,5 +1155,5 @@
       v.AssignSgn(a.v);
       // sign(f(x)) = 0;
-      dv = new SparseVector<object, V>();
+      dv = a.dv.Zero;
       return this;
     }
