Simple piSin piCos mathematical function in Squeak

Working with radians is great, because we have the derivatives of sine function being the cosine function, and derivative of cosine being the sine negated. I always ask engineers to work with radians, which avoids many mistakes!

But when it comes to handling the periodicity modulo 2π electronically, radians are not so great because π is irrational and hardly representable as a Float value!

With Attitude and Heading Reference System delivering vehicle heading in degrees or other fractions of π, conversion to radians somehow introduce unecessary numerical noise.

To work around this, many libm provide a sinpi(x) which computes sin(π x), and there is a similar sin_pi in boost library.

I was seeking for simplest implementation that could possibly work, so I decided to give it a try in Squeak Smalltalk: these mathematical functions are called like this:

     1 piCos = -1.
     (1/2) piSin = 1.

See http://www.squeaksource.com/NumberExtensions/Math-Number-Extensions-nice.6.diff.
And the tests http://www.squeaksource.com/NumberExtensions/Math-Number-ExtensionsTests-nice.1.diff.

It's a very simple implementation: reduce receiver in the interval -1/4,1/4 that is reduce the angle argument in the interval -π/4,π/4, then fallback to Float approximation of sin(π x) or cos(π x).

The nice thing is that we can afford to answer an exact value for every case when the sine or cosine value takes a rational value. This could only happen for infinitely many irrational receivers that we cannot represent (so we don't care), and a handful of rational fraction of π once we performed argument (receiver in our case) reduction. That's the Nivens theorem - See http://mathworld.wolfram.com/NivensTheorem.html and I would not have imagined that it could be so cheap before implementing it!

The only thing that bothers me, is that this exactness cost me the weak monotonicity property locally. For example, we may expect that:

    (1/6 + (1<<100) reciprocal) piSin >= (1/6) piSin.

But it's not the case... It happens that:

    (1/6 + (1<<100) reciprocal) asFloat < (1/6) and:
        [ (Float pi) < (1 asArbitraryPrecisionFloatNumBits: 100) pi ].

With such cumulated rounding by default, the floating point approximation falls short:

    ((1/6+(1<<100) reciprocal) * Float pi) sin = 0.5 predecessor.

This is where a bit more precision before the float conversion would help:

    ((1 asArbitraryPrecisionFloatNumBits: 100) pi / 6) asFloat sin = 0.5.

Weak monotonicity is preserved for Float receiver, which is already something, but I see no easy solution for the case of exact arithmetic. If I try with a longer approximation of π, then the weak monotonicity may break for mixed arithmetic, I'll have to come back to it... Fortunately, π/6 is the only case to handle (all reductions should be exact).

In Squeak there is a degreeCos and degreeSin message already that I tweaked some years ago, but those could further benefit from this new piSin and piCos functions (err, messages!). This little experiment might end up in Squeak trunk if there is some interest.