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 ε, where the really tiny part has the interesting property that ε2 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  = Dual;
using  = 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  Real { get; }
  public  Epsilon { get; }
  public  r => this.Real;
  public  ε => 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+1ε:

public Dual( real,  epsilon)
{
  this.Real = real;
  this.Epsilon = epsilon;
}

public static  EpsilonOne = new (0, 1);

The math is as we’ve discussed:

public static  operator +( x) => x;
public static  operator ( x) => new (x.r, x.ε);
public static  operator +( x,  y) =>
  new (x.r + y.r, x.ε + y.ε);
public static  operator ( x,  y) =>
  new (x.r  y.r, x.ε  y.ε);
public static  operator *( x,  y) =>
  new (x.r * y.r, x.ε * y.r + x.r * y.ε);
public static  operator /( x,  y) =>
  new (x.r / y.r, (x.ε * y.r  x.r * y.ε) / (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( x) => x > 2.0;

Then we reasonably expect that

static bool  Foo( 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( x,  y) =>
  x.r.CompareTo(y.r);
public
 int CompareTo( x) =>
  CompareTo(this, x);
public static bool operator <( x,  y) =>
  CompareTo(x, y) < 0;
public static bool operator >( x,  y) =>
  CompareTo(x, y) > 0;
public static bool operator <=( x,  y) =>
  CompareTo(x, y) <= 0;
public static bool operator >=( x,  y) =>
  CompareTo(x, y) >= 0;
public static bool operator ==( x,  y) =>
  CompareTo(x, y) == 0;
public static bool operator !=( x,  y) =>
  CompareTo(x, y) != 0;
public bool Equals( x) =>
CompareTo(this, x) == 0;

public override bool Equals(object obj) =>
  obj is  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 ( x) => new (x, 0);
public override string ToString() => $”{r}{(ε<0.0?“”:“+”)}{ε}ε”;
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 x4+2x3-12x2-2x+6:

static  Compute( 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  Compute( 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-16ε, 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+εto 3.5+ε 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…

Screen Shot 2018-12-17 at 5.36.45 PM.png

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?

Dual numbers, part 1

I’ve recently been looking into a fascinating corner of mathematics that at first glance appears a little bit silly, but actually has far-reaching applications, from physics to numerical methods to machine learning. I thought I’d share what I’ve learned over the next few episodes.

I assume you recall what a complex number is, but perhaps not all of the details. A complex number is usually introduced as a pair of real numbers (a, b), where a is called the “real part” and b is called the “imaginary part”.

A brief aside: it has always bugged me that these labels are unnecessarily value-laden. There is no particular “reality” that is associated with the real part; it is every bit as “imaginary” as the imaginary part. They might as well be called the “rezrov part” and the “gnusto part”, but we’re stuck with “real” and “imaginary”. Moving on.

It would be perfectly sensible to notate these as a pair (a, b) — and indeed, when we characterize complex numbers as points on the complex plane, we often do. But it is also convenient to notate complex numbers in an algebraic form a+bi, where the i indicates the imaginary part.

Finally, we have well-known rules for adding and multiplying complex numbers. The fundamental rule is that (0+1i)2 is equal to -1+0i, and the rest of the rules follow from there. For complex numbers:

(a+bi) + (c+di) = (a+c) + (b+d)i
(a+bi)(c+di)    = (ac-bd) + (ad+bc)i

Division is a little bit tricky and I won’t write it all out here; look it up if you do not remember. We can similarly go on to derive rules for complex exponentiation, trigonometry, and so on.

Let’s switch gears slightly. It is very common in mathematics, and especially calculus, to reason about a quantity “plus a little bit”. That is, any old quantity, plus some positive amount that is not zero, but is much smaller than one, and really, really small compared to the original quantity. What “really, really small” means is of course dependent on context; maybe one-millionth of the original quantity is small. Maybe one-millionth is still pretty big, but one-trillionth, or one part in a googolplex, or whatever, is “really small” in our context. I’m being deliberately vague and hand-wavy here; just run with it.

Let’s call our really small positive quantity ε, which is the lower-case Greek letter epsilon; it is traditionally used in calculus to represent a very small quantity. Let’s suppose our value is a, and we have a very small amount added to it, so that would be a+ε. Hmm. What if we doubled that quantity? Plainly that would be 2a+2ε. So scalar multiplication works as you’d expect. It seems reasonable that if we can have then we can have or for that matter 3.14ε or for any real value of b. So once again we have an (a, b) pair of numbers representing a quantity, which we notate as a+bε. What are the algebraic rules of this thing? It certainly makes sense that addition should work as we expect, the same as complex addition:

(a+bε) + (c+dε) = (a+c) + (b+d)ε

This should match our intuitions. The sum of “a plus a bit” and “c plus a bit” should be “a + c plus a bit” if there is justice in the world. What about multiplication? Let’s just multiply it out using ordinary algebra:

(a+bε)(c+dε) = ac + (ad+bc)ε + bdε2

We seem to have introduced an ε2 term there. What are we going to do with that? Remember, ε is not zero but is very, very small indeed, so that thing squared must be super tiny. So tiny that we can just ignore it entirely! (Remember, I am hand-waving egregiously here.) Let’s just say that ε2 is by definition zero, just as we said that i2 is by definition -1, and see what happens. So we have:

(a+bε)(c+dε) = ac + (ad+bc)ε

Again, this should match our intuitions. The product of something a bit bigger than a with something a bit bigger than c should be a bit bigger than the product of a and c.

What about division? Division by zero is illegal, and division by zero plus a tiny amount I think should also be illegal. Seems reasonable. What then is (a+bε)/(c+dε) when c is not zero? Well, if it is not zero then we can multiply top and bottom by c-dε:

a+bε   (c-dε)(a+bε)   ac + (bc-da)ε    a   bc-da
———— = ——————————— = —————————————— = - + ————— ε
c+dε   (c-dε)(c+dε)       c2           c     c2

Once again our intuition is satisfied: the quotient of “something close to a” divided by “something close to c” is close to the quotient of a and c.

This is a pretty neat number system; it is called the dual number system and it was introduced to the world in 1873 by the English mathematician William Clifford. But so far it might look like it is just a curiosity; however, it is extremely useful.

Next time on FAIC: we will implement the dual number system in C#, and discover a surprising property.

Accessibility of nested classes

Happy New Year all and welcome to another year of having fabulous adventures in coding. I thought I’d start this year off by answering a question I get quite frequently:

I have a public outer class with a public nested class. The nested class represents a concept only meaningful in the context of the outer class, and it accesses private members of the outer class. I want some of the members of the nested class to be public to everyone, but I want some of them to be public only to the outer class; how can I do that?

Typically the members they want to be accessible only to the containing class to be things like property setters and constructors:

public class Outer
{
  public class Inner
  {
    // I want the ctor and setter to be accessible in Outer
    // but not to the general public.
    public Inner() { ... }
    public Foo Foo { get; set; } 
    // whatever ... 
  }
  public Inner Frob() { ... }
  // ... and so on ...
}

I can think of four ways to solve this problem offhand, though I’m sure there are more.

Attempt one: my preferred solution is to deny the premise of the question and make the “private” details that need to be used by the container “internal”. It is your assembly! If you don’t want to abuse the internal details outside of a particular class, then don’t do it. If your coworkers do it, educate them in code reviews. Not every authorial intention needs to be expressed in the accessibility system, after all.

Attempt two: again, deny the premise of the question, but this time make the nested type a private implementation detail. The received wisdom is that we want nested classes to be private implementation details of the public outer class, not themselves public classes. In that case you’d simply make an interface to expose the public details:

public interface IFoo
{
  public Foo Foo { get; }
}
public class Outer
{
  private class Inner : IFoo
  {
    public Inner() { ... }
    public Foo Foo { get; private set; } 
    // whatever ... 
  }
  public IFoo Frob() { ... }
  // ... and so on ...
}

This has a couple of down sides though. We’ve denied the premise of the question, which is that the nested class makes sense to be public, and so not answered the question. We’ve also “polluted” the public namespace by adding a new IFoo concept unnecessarily. Also, there is now the possibility that just anyone can create an implementation of IFoo. Typically in these sorts of scenarios, the reason we want to restrict the constructor and setter is to ensure that the invariants of the inner type are being guaranteed by the code in the outer type.

Attempt three: did you know that it is legal for a type to implement an interface that is less public than the type? That is not true for base types, but it is true for interfaces!

public class Outer
{
  private interface IFooSetter
  {
    public Foo Foo { set; }
  }
  public class Inner : IFooSetter
  {
    public Inner() { ... }
    public Foo Foo { get; private set; } 
    Foo IFooSetter.Foo { set { this.Foo = value; } }
    // whatever ... 
  }
  public Inner Frob() { ... }
  // ... and so on ...
}

Now code in Outer can cast an Inner to IFooSetter and call the setter that way, but code outside of Outer cannot cast to IFooSetter because that type is private to Outer.

This looks pretty good, but again we’ve got a few problems. We still haven’t dealt with the constructor problem. And in both the solutions we’ve seen so far, it works poorly in the case where the inner type is a value type because we end up boxing.  Yes, it is a very bad practice to make a mutable value type, but hey, that’s probably why we’re restricting the mutators to be only accessible from Outer! Regardless though, we’d like a way to avoid the boxing penalty.

Attempt four: avoid interfaces altogether and go straight to delegates:

// This code is wrong!
public class Outer
{
  private static Func<Inner> innerFactory = null;
  private static Action<Inner, Foo> fooSetter = null;
  public class Inner
  {
    static Inner() 
    {
      Outer.innerFactory = ()=>new Inner();
      Outer.fooSetter = (i, f)=>i.Foo = f;
    }
    private Inner() { ... }
    public Foo Foo { get; private set; } 
    // ...
  }
  public Inner Frob() => innerFactory();
  // ... and so on ...
}

Now code in Outer can call innerFactory() any time it needs an Inner, and can call fooSetter(inner, foo) any time it needs to call a setter on an Inner.

I’ll leave you with two exercises:

  • First, the code above is wrong because it can be easily made to crash; do you see how? Can you fix it? (Hint: there is no way using just the rules of C# that I know of; you need to use a special method created for compiler writers.)
  • Second, suppose the inner type is a mutable value type; our setter would mutate a copy. Can you fix that too?