Questions on a Haskell -> C# conversion

Question

Background:

I was "dragged" into seeing this question: Fibonacci's Closed-form expression in Haskell
when the author initially tagged with many other languages but later focused to a Haskell question. Unfortunately I have no experience whatsoever with Haskell so I couldn't really participate in the question. However one of the answers caught my eye where the answerer turned it into a pure integer math problem. That sounded awesome to me so I had to figure out how it worked and compare this to a recursive Fibonacci implementation to see how accurate it was. I have a feeling that if I just remembered the relevant math involving irrational numbers, I might be able to work everything out myself (but I don't). So the first step for me was to port it to a language I am familiar with. In this case, I am doing C#.

I am not completely in the dark fortunately. I have plenty experience in another functional language (OCaml) so a lot of it looked somewhat familiar to me. Starting out with the conversion, everything seemed straightforward since it basically defined a new numeric type to help with the calculations. However I've hit a couple of roadblocks in the translation and am having trouble finishing it. I'm getting completely wrong results.

Analysis:

Here's the code that I'm translating:

data Ext = Ext !Integer !Integer
    deriving (Eq, Show)

instance Num Ext where
    fromInteger a = Ext a 0
    negate (Ext a b) = Ext (-a) (-b)
    (Ext a b) + (Ext c d) = Ext (a+c) (b+d)
    (Ext a b) * (Ext c d) = Ext (a*c + 5*b*d) (a*d + b*c) -- easy to work out on paper
    -- remaining instance methods are not needed

fib n = divide $ twoPhi^n - (2-twoPhi)^n
  where twoPhi = Ext 1 1
        divide (Ext 0 b) = b `div` 2^n -- effectively divides by 2^n * sqrt 5

So based on my research and what I can deduce (correct me if I'm wrong anywhere), the first part declares type Ext with a constructor that will have two Integer parameters (and I guess will inherit the Eq and Show types/modules).

Next is the implementation of Ext which "derives" from Num. fromInteger performs a conversion from an Integer. negate is the unary negation and then there's the binary addition and multiplication operators.

The last part is the actual Fibonacci implementation.

Questions:

In the answer, hammar (the answerer) mentions that exponentiation is handled by the default implementation in Num. But what does that mean and how is that actually applied to this type? Is there an implicit number "field" that I'm missing? Does it just apply the exponentiation to each corresponding number it contains? I assume it does the latter and end up with this C# code:

public static Ext operator ^(Ext x, int p) // "exponent"
{
    // just apply across both parts of Ext?
    return new Ext(BigInt.Pow(x.a, p), BigInt.Pow(x.b, p));
    //     Ext     (a^p)               (b^p)
}

However this conflicts with how I perceive why negate is needed, it wouldn't need it if this actually happens.

---

Now the meat of the code. I read the first part divide $ twoPhi^n - (2-twoPhi)^n as:

divide the result of the following expression: twoPhi^n - (2-twoPhi)^n.

Pretty simple. Raise twoPhi to the nth power. Subtract from that the result of the rest. Here we're doing binary subtraction but we only implemented unary negation. Or did we not? Or can binary subtraction be implied because it could be made up combining addition and negation (which we have)? I assume the latter. And this eases my uncertainty about the negation.

---

The last part is the actual division: divide (Ext 0 b) = b `div` 2^n. Two concerns here. From what I've found, there is no division operator, only a `div` function. So I would just have to divide the numbers here. Is this correct? Or is there a division operator but a separate `div` function that does something else special?

I'm not sure how to interpret the beginning of the line. Is it just a simple pattern match? In other words, would this only apply if the first parameter was a 0? What would the result be if it didn't match (the first was non-zero)? Or should I be interpreting it as we don't care about the first parameter and apply the function unconditionally? This seems to be the biggest hurdle and using either interpretation still yields the incorrect results.

Did I make any wrong assumptions anywhere? Or is it all right and I just implemented the C# incorrectly?

Code:

Here's the (non-working) translation and the full source (including tests) so far just in case anyone is interested.

// code removed to keep post size down
// full source still available through link above

Progress:

Ok so looking at the answers and comments so far, I think I know where to go from here and why.

The exponentiation just needed to do what it normally does, multiply p times given that we've implemented the multiply operation. It never crossed my mind that we should do what math class has always told us to do. The implied subtraction from having addition and negation is a pretty handy feature too.

Also spotted a typo in my implementation. I added when I should have multiplied.

// (Ext a b) * (Ext c d) = Ext (a*c + 5*b*d) (a*d + b*c)
public static Ext operator *(Ext x, Ext y)
{
    return new Ext(x.a * y.a + 5*x.b*y.b, x.a*y.b + x.b*y.a);
    //                 ^ oops!
}

Conclusion:

So now it's completed. I only implemented to essential operators and renamed it a bit. Named in a similar manner as complex numbers. So far, consistent with the recursive implementation, even at really large inputs. Here's the final code.

static readonly Complicated TWO_PHI = new Complicated(1, 1);
static BigInt Fib_x(int n)
{
    var x = Complicated.Pow(TWO_PHI, n) - Complicated.Pow(2 - TWO_PHI, n);
    System.Diagnostics.Debug.Assert(x.Real == 0);
    return x.Bogus / BigInt.Pow(2, n);
}

struct Complicated
{
    private BigInt real;
    private BigInt bogus;

    public Complicated(BigInt real, BigInt bogus)
    {
        this.real = real;
        this.bogus = bogus;
    }
    public BigInt Real { get { return real; } }
    public BigInt Bogus { get { return bogus; } }

    public static Complicated Pow(Complicated value, int exponent)
    {
        if (exponent < 0)
            throw new ArgumentException(
                "only non-negative exponents supported",
                "exponent");

        Complicated result = 1;
        Complicated factor = value;
        for (int mask = exponent; mask != 0; mask >>= 1)
        {
            if ((mask & 0x1) != 0)
                result *= factor;
            factor *= factor;
        }
        return result;
    }

    public static implicit operator Complicated(int real)
    {
        return new Complicated(real, 0);
    }

    public static Complicated operator -(Complicated l, Complicated r)
    {
        var real = l.real - r.real;
        var bogus = l.bogus - r.bogus;
        return new Complicated(real, bogus);
    }

    public static Complicated operator *(Complicated l, Complicated r)
    {
        var real = l.real * r.real + 5 * l.bogus * r.bogus;
        var bogus = l.real * r.bogus + l.bogus * r.real;
        return new Complicated(real, bogus);
    }
}

And here's the fully updated source.

Answer 1

[...], the first part declares type Ext with a constructor that will have two Integer parameters (and I guess will inherit the Eq and Show types/modules).

Eq and Show are type classes. You can think of them as similar to interfaces in C#, only more powerful. deriving is a construct that can be used to automatically generate implementations for a handful of standard type classes, including Eq, Show, Ord and others. This reduces the amount of boilerplate you have to write.

The instance Num Ext part provides an explicit implementation of the Num type class. You got most of this part right.

[the answerer] mentions that exponentiation is handled by the default implementation in Num. But what does that mean and how is that actually applied to this type? Is there an implicit number "field" that I'm missing? Does it just apply the exponentiation to each corresponding number it contains?

This was a bit unclear on my part. ^ is not in the type class Num, but it is an auxilliary function defined entirely in terms of the Num methods, sort of like an extension method. It implements exponentiation to positive integral powers through binary exponentiation. This is the main "trick" of the code.

[...] we're doing binary subtraction but we only implemented unary negation. Or did we not? Or can binary subtraction be implied because it could be made up combinding addition and negation (which we have)?

Correct. The default implementation of binary minus is x - y = x + (negate y).

The last part is the actual division: divide (Ext 0 b) = b `div` 2^n. Two concerns here. From what I've found, there is no division operator, only a div function. So I would just have to divide the numbers here. Is this correct? Or is there a division operator but a separate div function that does something else special?

There is only a syntactic difference between operators and functions in Haskell. One can treat an operator as a function by writing it parenthesis (+), or treat a function as a binary operator by writing it in `backticks`.

div is integer division and belongs to the type class Integral, so it is defined for all integer-like types, including Int (machine-sized integers) and Integer (arbitrary-size integers).

I'm not sure how to interpret the beginning of the line. Is it just a simple pattern match? In other words, would this only apply if the first parameter was a 0? What would the result be if it didn't match (the first was non-zero)? Or should I be interpreting it as we don't care about the first parameter and apply the function unconditionally?

It is indeed just a simple pattern match to extract the coefficient of √5. The integral part is matched against a zero to express to readers that we indeed expect it to always be zero, and to make the program crash if some bug in the code was causing it not to be.

---

A small improvement

Replacing Integer with Rational in the original code, you can write fib n even closer to Binet's formula:

fib n = divSq5 $ phi^n - (1-phi)^n
  where divSq5 (Ext 0 b) = numerator b
        phi = Ext (1/2) (1/2)

This performs the divisions throughout the computation, instead of saving it all up for the end. This results in smaller intermediate numbers and about 20% speedup when calculating fib (10^6).

Answer 2

First, Num, Show, Eq are type classes, not types nor modules. They are a bit similar to interfaces in C#, but are resolved statically rather than dynamically.

Second, exponentiation is performed via multiplication with the implementation of ^, which is not a member of the Num typeclass, but a separate function.

The implementation is the following:

(^) :: (Num a, Integral b) => a -> b -> a
x0 ^ y0 | y0 < 0    = error "Negative exponent"
        | y0 == 0   = 1
        | otherwise = f x0 y0
    where -- f : x0 ^ y0 = x ^ y
          f x y | even y    = f (x * x) (y `quot` 2)
                | y == 1    = x
                | otherwise = g (x * x) ((y - 1) `quot` 2) x
          -- g : x0 ^ y0 = (x ^ y) * z
          g x y z | even y = g (x * x) (y `quot` 2) z
                  | y == 1 = x * z
                  | otherwise = g (x * x) ((y - 1) `quot` 2) (x * z)

This seems to be the missing part of solution.

You are right about subtraction. It is implemented via addition and negation.

Now, the divide function divides only if a equals to 0. Otherwise we get a pattern match failure, indicating a bug in the program.

The div function is a simple integer division, equivalent to / applied to integral types in C#. There is also an operator / in Haskell, but it indicates real number division.

Answer 3

A quick implementation in C#. I implemented exponentiation using the square-and-multiply algorithm.

It is enlightening to compare this type which has the form a+b*Sqrt(5) with the complex numbers which take the form a+b*Sqrt(-1). Addition and subtraction work just the same. Multiplication is slightly different, because i^2 isn't -1 but +5 here. Division is slightly more complicated, but shouldn't be too hard either.

Exponentiation is defined as multiplying a number with itself n times. But of course that's slow. So we use the fact that ((a*a)*a)*a is identical to (a*a)*(a*a) and rewrite using the square-and-multiply algorithm. So we just need log(n) multiplications instead of n multiplications.

Just calculating the exponential of the individual components doesn't work. That's because the matrix underlying your type isn't diagonal. Compare this to the property of complex numbers. You can't simply calculate the exponential of the real and imaginary part separately.

struct MyNumber
{
    public readonly BigInteger Real;
    public readonly BigInteger Sqrt5;

    public MyNumber(BigInteger real,BigInteger sqrt5)
    {
        Real=real;
        Sqrt5=sqrt5;
    }

    public static MyNumber operator -(MyNumber left,MyNumber right)
    {
        return new MyNumber(left.Real-right.Real, left.Sqrt5-right.Sqrt5);
    }

    public static MyNumber operator*(MyNumber left,MyNumber right)
    {
        BigInteger real=left.Real*right.Real + left.Sqrt5*right.Sqrt5*5;
        BigInteger sqrt5=left.Real*right.Sqrt5 + right.Real*left.Sqrt5;
        return new MyNumber(real,sqrt5);
    }

    public static MyNumber Power(MyNumber b,int exponent)
    {
        if(!(exponent>=0))
            throw new ArgumentException();
        MyNumber result=new MyNumber(1,0);
        MyNumber multiplier=b;
        while(exponent!=0)
        {
            if((exponent&1)==1)//exponent is odd
                result*=multiplier;
            multiplier=multiplier*multiplier;
            exponent/=2;
        }
        return result;
    }

    public override string ToString()
    {
        return Real.ToString()+"+"+Sqrt5.ToString()+"*Sqrt(5)";
    }
}

BigInteger Fibo(int n)
{
    MyNumber num = MyNumber.Power(new MyNumber(1,1),n)-MyNumber.Power(new MyNumber(1,-1),n);
    num.Dump();
    if(num.Real!=0)
      throw new Exception("Asser failed");
    return num.Sqrt5/BigInteger.Pow(2,n);
}

void Main()
{
  MyNumber num=new MyNumber(1,2);
  MyNumber.Power(num,2).Dump();
  Fibo(5).Dump();
}