Trig Tests

Overhauling the trig tests turned out to be much more epic than I had intended. Guess that’s one of the hazards of writing the test generating code in Perl 6 while Rakudo is still under heavy development.

The old test code was very, very thorough in testing each configuration. Basically, we looped through a list of angles, and tested each angle with every possible way it could be called for each type we wanted to test. So say the angle is 45 degrees. Then for the Num type, we’d call 45.Num.sin(Degrees), sin(45.Num, Degrees), and sin(:x(45.Num), Degrees), and then we’d call the same for the equivalent of 45 in Radians, Gradians, and Circles. Then we’d start it all over again using Rat instead of Num. And so on and so forth.

In practice, any implementation is likely to redispatch the trig functions to two basic sets of methods, one for reals and one for complex. So that’s exactly what the tests do now (leaving out the complex case because it is long:

for TrigTest::cosines() -> $angle
{
    my $desired-result = $angle.result;

    # Num.cos tests -- very thorough
    is_approx($angle.num(Radians).cos, $desired-result, 
              "Num.cos - {$angle.num(Radians)} default");
    for TrigTest::official_bases() -> $base {
        is_approx($angle.num($base).cos($base), $desired-result, 
                  "Num.cos - {$angle.num($base)} $base");
    }
}

Then, instead of trying to exhaustively test the rest of the cases, we try to generate one test for each case, and leave it at that. For instance, here’s the Rat tests for cos:

# Rat tests
is_approx((-0.785398163404734).Rat(1e-9).cos, 0.707106781186548, "Rat.cos - -0.785398163404734");
is_approx((0).Rat(1e-9).cos(Circles), 1, "Rat.cos(Circles) - 0");
is_approx((0.785398163404734).Rat(1e-9).cos(:base(Radians)), 0.707106781186548, "Rat.cos(:base(Radians)) - 0.785398163404734");
is_approx(cos((1.57079632680947).Rat(1e-9)), 0, "cos(Rat) - 1.57079632680947");
is_approx(cos((135).Rat(1e-9), Degrees), -0.707106781186548, "cos(Rat, Degrees) - 135");
is_approx(cos(:x((3.14159265361894).Rat(1e-9))), -1, "cos(:x(Rat)) - 3.14159265361894");
is_approx(cos(:x((250).Rat(1e-9)), :base(Gradians)), -0.707106781186548, "cos(:x(Rat), :base(Gradians)) - 250");

The use of .Rat(1e-9) is to ensure we have a Rat argument — longer decimal numbers may become Nums otherwise. This tests an additional case (Rat.cos(:base(Radians))) that was never actually tested before. Despite that, the number of Rat tests overall goes from something like 160 to 7.

That’s good, because I added three more types to be tested. Two represent generic Cool types: Str and NotComplex. Str is used to test real numbers in a non-Numeric Cool type. NotComplex is a non-Numeric Cool type which returns a Complex which you call .Numeric on it. The final type, DifferentReal, is a Real type which is not built-in.

Much to my surprise, I actually turned up a couple of bugs in Rakudo’s trig implementation when I tried this testing approach with atan2, the oddball trig function which normally takes two real arguments. First, we were not being generous enough on what the second argument would accept, which was easily corrected. Second, our fake enum standing in for TrigBase looks exactly like an Int to Rakudo. So 2.atan2(Degrees) is indistinguishable from 2.atan2(1), which is definitely not what the user wants. Luckily that’s a pretty silly way to call atan2 anyway.

The process of getting this to work turned up a few Rakudo bugs in the non-trig code, which have been reported. I won’t go into the details here.

By my reckoning, this means there is only one thing left on my list for the numeric grant, srand and rand. I will try to finish them off tomorrow. Thanks to The Perl Foundation and Ian Hague for supporting this work.

Detour from Sin

I spent the first two days of YAPC furiously working on overhauling the trig tests every chance I got. On Wednesday, however, I suddenly changed my focus, and I’ve been working on something else ever since. What happened?

Simple: my improved trig tests ran up hard against a major Rakudo limitation. So I’ve been working on fixing it so I can get back to working on those tests. The issue was this: on a 32-bit system, you could only enter about nine digits in your Num constants. The tenth digit would make the action for handling a Rat or Num blow up, and either give an error or incorrect results.

This is something of a problem when you are automatically generating code, because writing a Num usually gets you 15 digits.

So, I set out to fix Rakudo’s Rat / Num constant handling. First I hand wrote a parser, which I ended up not using. Then I coded up a a beautiful functional approach to handling converting the string components to a number. Something like this (for Rat):

$result = [+] $int-part.comb(/\d/).reverse.kv.map(-> $i, $d { $d.Int * 10 ** $i });
$result += [+] $frac-part.comb(/\d/).kv.map(-> $i, $d { $d.Int / 10 ** ($i + 1) }); 

(Note: reconstructed from memory, may not actually work as given.) This worked great, and had the happy property of properly returning a Rat if possible, and promoting to a Num if needed for additional precision. (Basically, because all the basic components of the math are simple ints, and the math operations correctly promote as needed.)

There was only one problem with this approach: it took eight times as long to do its job as the original. I tried rewriting it in a more iterative fashion. Then it only took seven times as long. Still not acceptable.

I toyed around with other ideas for reworking the solution, but finally, at about the fourth hour of my drive home from YAPC, I hit on the problem with this approach. Every one of those simple math operations was necessarily a multi call. Three of those were required for each digit. No wonder it was slow!

That in mind, I decided on a new approach. The basic problem with the original Rakudo approach was that it treated the various components of the numeric constant as Ints, which can easily (and sometimes ungracefully) overflow. But there is a standard Rakudo trick for handling overflow in Ints gracefully: do the math internally as a Num, and then cast the result to Int only if possible. This can be done entirely in PIR for speed — no worrying about multis for the math, it’s all done explicitly in num.

Then the resulting bits are combined in Perl 6, so that multi dispatch on the math operations gets you the correct type results. It took a bit of doing, but I just pushed the resulting code to github.

Now, this solution isn’t perfect. I’m pretty sure it’s possible to construct extreme cases where the step-by-step algorithm would have returned a Rat, but the new one returns a Num. And if you actually have 64-bit Ints available, you might lose a few digits of precision, because Num have less — but this is a problem throughout Rakudo, not just here.

But in general, this solution should be acceptably fast and come much closer to doing the right thing than the previous version did.

Trig Tests

When I originally wrote the expanded trig tests, the tests were very thorough. Given an angle, it would (say), test that angle in the default trig base, and then in radians, degrees, gradians, and “circles”. For three different ways to call it for Num, then for Rat, Complex, and Int. This led to a lot of trig tests: sin/asin had 1782 combined tests, for instance. But I still needed to add more tests! So I’ve known for a while I was going to have to reduce the number of tests. I just wasn’t sure how.

I played around with several different methods for reducing the count. I settled early on on the notion that the full tests would still be performed for the method forms of Num and Complex. That way we have confidence the underlying algorithms do work. My current effort is to get the code to automatically generate one test for each form for each type. So for instance

# Rat tests
is_approx((-3.9269908).Rat.sin, 0.7071067, "Rat.sin - -3.9269908 default");
is_approx((-0.625).Rat.sin(Circles), 0.7071067, "Rat.sin(Circles) - -0.625 default");
is_approx((-3.9269908).Rat.sin(:base(Radians)), 0.7071067, "Rat.sin(:base(Radians)) - -3.9269908 default");
is_approx(sin((-3.9269908).Rat), 0.7071067, "sin(Rat) - -3.9269908 default");
is_approx(sin((-225).Rat, Degrees), 0.7071067, "sin(Rat) - -225 default");
is_approx(sin(:x((-3.9269908).Rat)), 0.7071067, "sin(:x(Rat)) - -3.9269908 default");
is_approx(sin(:x((-250).Rat), :base(Gradians)), 0.7071067, "sin(:x(Rat), :base) - -250 default");

Except this is wrong: I want to test a different angle for each form. I need to do a bit of refactoring to make that possible, but I think I should be able to get it going today.

One trick I am using here might qualify as a new Perl 6 idiom. I have a list of angles to test, and a list of trig bases. I want to come up with lots of combinations. My solution has been to do something like this:

 my $base_list = (<Radians Degrees Gradians Circles> xx *).flat;

That gets you an infinite list that just repeats those four strings over and over again — you just .shift them off as needed. You could do the same thing with array indexing and modular arithmetic, but IMO this is more elegant. (Though now that I think about it, with the array I believe you could say something like @array[$count++ % *] and it would work.)

Thanks to The Perl Foundation and Ian Hague for supporting this work.

Accepting Equality

Recently I realized that Numeric and Real need to implement ACCEPTS. As far as I know, the Spec has nothing to say about what this ACCEPTS should do, so I looked at the current implementation. There are two versions of Num.ACCEPTS, one which takes anything and uses == to compare, and one which takes a Complex and checks that the imaginary part is zero, and the real part is equal to the Num. Both ACCEPTS have special cases for NaN, needed because NaN != NaN.

Talking about this with Jonathan, I was initially stumped. If every class had to have an ACCEPTS for its own type, for Real, and for Numeric, things would be a mess. Then it occurred to me: why not implement it by simply calling to ==? As long as == is doing the right thing, ACCEPTS will do the right things as well. (Errr, modulo the NaN issue, which might get tricky.)

Of course, that depended on a notion of “right thing” for ==. Quick check: does a Real number equal the Complex version of the same number?

> 1 == 1 + 0i
0

Huh. Well, let’s get a second opinion:

colomon: alpha: say 1 == 1 + 0i
p6eval: alpha 30e0ed: OUTPUT«1␤»

Whoops. The False answer in current Rakudo is because of my reworking of Numeric comparison. At the time, I was thinking of sorting, and figured it would do no harm if 1 < 1+0i. But of course that looks silly in the context of ==! Apparently the test suite doesn’t have any tests for this case.

As frequently happens with these posts, I started this one without any clear sense of where I should go, but now I’ve got a notion:

1) Fix == so that Reals and equivalent Complexes are equal again.

2) Make a generic ACCEPTS for Numeric that works based on ==.

3) Somehow fold NaN handling into the mix.

Hmmm. Okay, my actual first step is going to be to completely muck up Numeric.ACCEPTS (intentionally), and see what, if any, tests fail, so I can find out the state of testing ACCEPTS.

Thanks to the Perl Foundation and Ian Hague for supporting this work.

div and mod

After consideration and looking at what Rakudo actually does, I decided the right solution to the div and mod issue was to change them in the Spec. In particular, I see them as having one interesting use: with integer types, they can always return the same integer type as their result, because either’s result will never have greater magnitude than its operands. I’ve rewritten their definitions to correspondingly tighten them as only working on integer types.

Real Operators

I decided to work on the operators for Real tonight. First up was adding tests I skipped writing before, for the Real / Bridge versions of infix +, *, and /. Then I implemented a completely generic version of infix % for Real, and added tests for it.

And then I ran into the tricky bit: div and mod. div is division which generates a result of the same type as the two operands. So for instance, 5 / 4 == 1.25 — two Ints divided to get a Rat. But 5 div 4 == 1. This is perfectly straightforward and sensible, but how do you implement a generic version that does this? It seems like each type descended from Real will have to define its own div. And it won’t necessarily be easy, either — I’m not sure how to sensibly do it for Rat, for instance. You can’t fall back to outputting a Num if things get out of hand.

It gets even worse with mod. The spec says that for built-in types, $x mod $y == $x - ($x div $y) * $y. That makes mod the equivalent of % for integer types. But for non-integer types, div should be closer to normal division, in which case mod will always return zero or possibly a small error value. IMO the spec is poorly defined here, but I’m not yet seeing how to improve it.

Thanks to the Perl Foundation and Ian Hague for supporting this work.

Comparing Two Numerics

I’ve mentioned before that TimToady and I had ideas for comparing two arbitrary Numeric values. The scheme may not be mathematically ideal, but it does define a handy ordering which works predictably for all Numeric types and is easy to understand.

The idea is that any Numeric value should be able to return a list of real values to use for the purposes of comparison. For a Real, the list’s only element is the value itself. For a Complex, the list has two parts, the real and imaginary values of the Complex, in that order. For your user-defined Numeric type, its whatever list of Reals seems to fit well with the above scheme.

Then, when asked to compare to Numeric values, we generate that list of Reals for each and do a lexicographical compare on the two lists. (We don’t go into an infinite recursion here because comparing two Reals calls different code.)

I’ve been kicking around this notion for a long time now without being able to think of an appropriate name for the method to return a list of Reals. I realized yesterday that the problem was I was thinking of it as being akin to the Real Bridge method — should I call it Bridges, or something like that. But in fact, there’s no need to get Bridge involved at all, we really are just returning a list of Reals. So the perfect name (IMO) is reals.

With that out of the way, I started implementing the code as soon as I had some tuits. Here’s my first working stab at it:

multi sub infix:«<=>»(Numeric $a, Numeric $b) {
    my @a = $a.reals;
    my @b = $b.reals;
    for ^(+@a max +@b) -> $i {
        return -1 if $i >= +@a;
        return +1 if $i >= +@b;
        my $c = @a[$i] <=> @b[$i];
        return $c if $c != 0;
    }
    0;
}

Kind of messy, isn’t it? Then I realized that there was a simple meta-op-based version:

multi sub infix:«<=>»(Numeric $a, Numeric $b) {
    my @a = $a.reals;
    my @b = $b.reals;
    [||] (@a Z<=> @b), (+@a <=> +@b);
}

That right there is one of the major reasons I love Perl 6. Perl 5 had a lovely idiom for chained comparisons, but it required you fix the length of the chain in advance. Perl 6 lets you elegantly extend that idea to handle chains of varying length. (If you’re wondering, the Z comparison does pairwise comparison on the two arrays until one of them runs out, then we tag on one additional -1/0/+1 to indicate which array was shorter, and then run a chained || on the entire lot.

Warning: Neither code example here is fully tested yet! But I was so excited about the meta-op version I had to share.

Thanks to the Perl Foundation and Ian Hague for supporting this work.

Pred to sin some more

I’ve accomplished a good bit of work on Numeric in the last eleven days. Soon after my last post I started work on comparison operators. First I provided Real versions of the numeric comparison operators, all of which delegate to Bridge versions (implemented as Num for now). Then I removed the Num, Num version of the cmp operator, replacing it with a Numeric, Numeric version. Note that despite the presence of Numeric, Numeric cmp, numeric comparisons only work on Real types at the moment.

TimToady++ and I have a theory for supporting comparisons on generic Numeric types, but I have not started to implement it yet. My basic idea is similar to the Bridge idea for Reals. There would be a Numeric method (don’t have a name yet) which converts your generic Numeric type to a list of Real numbers. (Obviously Real numbers would convert to a single Real, Complex would convert to two, etc.) You would then compare the Reals lexicographically. This probably isn’t a mathematically sensible ordering, but it does provide a reasonable ordering that would put Reals in the order you expect and still allow you to sort Complex numbers.

After that I got sqrt working properly in the new system, adding Numeric.sqrt and Real.sqrt, and cleaning up the Num and Complex versions.

Then the big effort (at least in terms of lines of code) was getting the trig functions working properly in the new system. It was just a matter of cloning the approach I’d taken for sin for all the other functions, though there was one bizarre catch. It turned out that several of the trig functions had two different implementations in Complex. When I took away the “multi” from the first implementation, the second one kicked in — and the second one was wrong. Since I didn’t actually realize the second one was there, this lead to about an hour of extreme confusion on my part before I sorted things out.

I finished up with tweaks for ceiling, floor, and log. Then I implemented Numeric versions of postfix i, succ, and pred. These are expressed in terms of simple addition and multiplication, so Numeric’s implementations should work for any reasonably robust Numeric class. I’m guessing it will frequently be worthwhile to still have class-specific versions of succ and pred, both for efficiency and to help ensure that $x and $x.pred usually have the same type.

What’s next? As I mentioned about, generic Numeric comparison code. The rand function — it’s currently listed as Numeric, but it doesn’t take any arguments and returns a Num, so I’m not sure why. srand at least takes a Real, but again, it doesn’t really seem like it fits, IMO. Numeric.roots. Work on more tests. Fixing up the spec a bit. And I may take a look in other math libraries to see if there are obvious functions we are missing.

Thanks to the Perl Foundation and Ian Hague for supporting this work.

Numeric Progress Report

Since my last report on the Numeric grant, most of my time has actually gone to moving my family. We are now mostly settled in two hours’ north of where we lived a month ago. In that time I made less progress than I hoped to on Numeric, but now things should be well positioned to make some real progress.

So far I’ve focused on the unpolar methods, the rounding methods, and the trig support methods. The unpolar methods, cis and unpolar, I have moved to Real, both in the spec and in the code. Implemented using the new Bridge method (my current method name for the lingua franca type I was calling RealX in previous posts, not yet specced), cis and unpolar should work for any type that does the Real role and implements a working version of Bridge. Similarly, the rounding methods, ceiling, floor, truncate, and round, are also all implemented generically in Real now. Both ceiling and floor forward to the bridge type versions of those functions, while truncate and round are built on ceiling and floor.

At Moritz’s prodding, I also moved to-radians and from-radians from Any to Numeric, made them public, and specced them. In my experience, if you’re working with trig functions, sooner or later you will need those functions. If we already have to implement them for ourselves, no reason we shouldn’t save every other Perl 6 programmer from having to do the same.

At the moment, there are several Rakudo bugs concerning roles which need to be sorted out before Numeric and Real can be fully implemented. I’m hoping to have time to work on this with Jonathan this week.

Thanks to the Perl Foundation and Ian Hague for supporting this work.

Numeric Method Rethinking

Just looking over the spec, here are four methods I think might need to relocate. I’m posting this here just to see if anyone objects to me changing the spec.

cis and unpolar are both listed as being in Numeric, but their invocants are specified as being Real. I think they both belong in Real, because conceptually they are for mapping Reals to Complex.

sqrt is listed in Real, but of course there is a version for Complex as well. roots on the other hand, is Numeric. It seems to me conceptually these are basically the same thing. I’m not sure what the right answer is here. I’ve been thinking of Quaternions as my canonical example of a non-Real, non-Complex Numeric type. I believe they have a solid concept of square roots. I don’t know if they have a concept of N-root that corresponds to the roots function.

Anyone else have thoughts on these?