Dual numbers, part 2

Last time on FAIC I introduced a fun variation on complex numbers called "dual numbers" consisting of a value "plus a tiny bit" notated as I, where the really tiny part has the interesting property that I² is zero.

Well, enough chit-chat; let's implement them.

You know what, just for grins I'm going to make this look a bit more like math, and a bit more concise while we're at it:

using a...... = Dual;
using a = System.Double;

(Unfortunately " is not a legal identifier in C#.)

A dual is logically a value that consists of two doubles, so it makes sense to represent them as an immutable value type:

public struct Dual : IComparable<Dual>
{
public a Real { get; }
public a Epsilon { get; }
public a r => this.Real;
public a I => this.Epsilon;

I got so sick of typing out "Real" and "Epsilon" that I had to make these helper properties, even though they do not meet C# naming guidelines. Guidelines are guidelines, not rules!

Notice also that, fun fact, C# lets you use any non-surrogate-pair Unicode character classified as a letter in an identifier.

The constructor is as you'd expect, and I'm going to make a helper property for 0+1I:

public Dual(a real, a epsilon)
{
this.Real = real;
this.Epsilon = epsilon;
}
public static a...... EpsilonOne = new a......(0, 1);

The math is as we've discussed:

public static a...... operator +(a...... x) => x;
public static a...... operator -(a...... x) => new a......(-x.r, -x.I);
public static a...... operator +(a...... x, a...... y) =>
new a......(x.r + y.r, x.I + y.I);
public static a...... operator -(a...... x, a...... y) =>
new a......(x.r - y.r, x.I - y.I);
public static a...... operator *(a...... x, a...... y) =>
new a......(x.r * y.r, x.I * y.r + x.r * y.I);
public static a...... operator /(a...... x, a...... y) =>
new a......(x.r / y.r, (x.I * y.r - x.r * y.I) / (y.r * y.r));

So easy!

UPDATE: I was WRONG WRONG WRONG in the section which follows.

An earlier version of this episode suggested that we implement comparisons on dual numbers as "compare the real parts; if equal, compare the epsilon parts". But upon reflection, I think that's not right. One of the characteristics of dual numbers that we like is that "lifting" a computation on reals to dual numbers produce the same "real part" in the result. Suppose we have:

static bool Foo(a x) => x > 2.0;

Then we reasonably expect that

static bool Foo(a...... x) => x > 2.0;

agrees with it, no matter what value we have for epsilon.

So I'm going to tweak the code so that we compare only real parts.

Resuming now the original episode"

My preference is to implement the logic in one place in a helper, and then call that helper everywhere:

private static int CompareTo(a...... x, a...... y) =>
x.r.CompareTo(y.r);
public int CompareTo(a...... x) =>
CompareTo(this, x);
public static bool operator <(a...... x, a...... y) =>
CompareTo(x, y) < 0;
public static bool operator >(a...... x, a...... y) =>
CompareTo(x, y) > 0;
public static bool operator <=(a...... x, a...... y) =>
CompareTo(x, y) <= 0;
public static bool operator >=(a...... x, a...... y) =>
CompareTo(x, y) >= 0;
public static bool operator ==(a...... x, a...... y) =>
CompareTo(x, y) == 0;
public static bool operator !=(a...... x, a...... y) =>
CompareTo(x, y) != 0;
public bool Equals(a...... x) =>
CompareTo(this, x) == 0;
public override bool Equals(object obj) =>
obj is a...... x && CompareTo(this, x) == 0;

And finally, a few loose ends. It would be nice to be able to convert to Dual from double automatically, and we should also override ToString and GetHashCode just to be good citizens:

public static implicit operator a......(a x) => new a......(x, 0);
public override string ToString() => $"{r}{(I<0.0?"":"+")}{I}I";
public override int GetHashCode() => r.GetHashCode();

Super, that was really easy. And now we can take any old method that does math on doubles, and make it do math on Duals. Suppose we've got this little guy, that computes x⁴+2x³-12x²-2x+6:

static a Compute(a x) =>
x * x * x * x + 2 * x * x * x - 12 * x * x - 2 * x + 6;

If we just turn all the doubles into Duals:

static a...... Compute(a...... x) =>
x * x * x * x + 2 * x * x * x - 12 * x * x - 2 * x + 6;

Then we have the same function, but now implemented in the dual number system. This:

Console.WriteLine(Compute(1.0 + Dual.EpsilonOne));

Produces the output -5-16I, which agrees with the original method in the real part, and is -16 in the epsilon part. Apparently computing that polynomial with one plus a tiny amount gives us -5, plus -16 "more tiny amounts".

Hmm.

You know what, just for grins I'm going to compute this function from -4.5+Ito 3.5+I and then graph the real and epsilon parts of the output of our "dualized" function.

The blue line is the real part, the orange line is the epsilon part. As we expect, the blue line corresponds to the original polynomial. But the orange line looks suspiciously familiar"

HOLY GOODNESS THE ORANGE LINE IS ZERO WHERE THE BLUE LINE IS FLAT.

The epsilon part is the derivative of the polynomial!

Next time on FAIC: How'd that happen?