Hacker News Reader

Comment by nsajko | Hacker News Reader

The video is very well produced and engaging, however (from skimming it), it seems that all the implementation methods it presents are suboptimal.

The N64 CPU had an FPU, so it seems obvious that the best method, from both a speed and accuracy (or whichever is deemed more important) standpoint would be the same method that's used in all the libms in libcs like Glibc, musl or LLVM-libc, or in the standard libraries of languages like Go or Java: approximation using a polynomial. This polynomial is not derived analytically, rather it's found using an iterative algorithm, traditionally using the Remez algorithm, but also see my WIP project here: https://gitlab.com/nsajko/FindMinimaxPolynomial.jl.

View on HN

> This polynomial is not derived analytically, rather it's found using an iterative algorithm, traditionally using the Remez algorithm, but also see my WIP project here: https://gitlab.com/nsajko/FindMinimaxPolynomial.jl.

Is this your independent attempt? Because I think RLIBM did the essentially same thing recently [1].

[1] https://arxiv.org/pdf/2104.04043.pdf

One difference is that their approach doesn't scale to 64-bit floats AFAIK. My code and their approach are actually quite different, with different goals, different places where they're smart and different places where they could be improved.

I hope to improve my code based on their ideas when I finally read their papers, however I still have plenty of my own that I want to implement.

See also this thread from twenty days ago, someone asked me basically the same thing, linking an earlier paper by Lim and Nagarakatte (et al):

https://news.ycombinator.com/item?id=36201197

Cool! I think the RLIBM approach is indeed fresh and has potentials, but I was never sure if it can be generalized to domains past binary32 when we are looking for correctly rounded results. I guess it's a fundamental difference between RLIBM and your algorithm---you don't guarantee correctly rounded results, right?

> I guess it's a fundamental difference between RLIBM and your algorithm---you don't guarantee correctly rounded results, right?

Yes, but note that the linked package doesn't even attempt to implement a polynomial approximation (for use in software), rather it just finds the coefficients of a polynomial. Another WIP package of mine will use FindMinimaxPolynomial to actually implement the approximation (with source code output etc), however even that package was never meant to guarantee correctly rounded results.

In particular their work is more ambitious than mine in that they seem to account for the range reduction and output compensation while looking for the coefficients with LP, which enables them to guarantee correctly rounded results.

I really don't think you can claim that the solution is suboptimal. Kaze has clearly spent a lot of time on this and has a real code to show for it. You're just theorizing that it's bad without providing any real evidence.

I'm not claiming anything extraordinary, though, most or all of what I said is well known in the field.

The burden of proof is on you though. If you think you can do it better then do it. It's ridiculous to post a dismissal without even watching the full video.

A good reference is the Handbook of Floating-Point Arithmetic. If you can be more specific about what is it that's not clear, I'll give more specific "proof". Even though you're being a bit "sealion"-y.

> It's ridiculous to post a dismissal without even watching the full video.

The video is jarring to watch. It's well-produced, engaging, etc.; but also pretentious and overconfident. It seems to be chock-full with irrelevant "developments". So it's not reasonable to expect me to watch it end-to-end.

> The video is jarring to watch. It's well-produced, engaging, etc.; but also pretentious and overconfident. It seems to be chock-full with irrelevant "developments".

Don't even try to read the YouTube comments :-) “Kaze has to be the best assembly programmer of all time”, “this is the greatest genius since the inverse square root trick”, “this SM64 hack is the most optimized Nintendo 64 game of all time”, etc…

Polynomial approximation was used throughout the video, and minimax was one of the final methods. Timestamp is sometime around 18:00. There seemed to be numerical issues with the approximation producing sin values above 1.0.

Yeah, turns out I missed that part of the video. So the author, in fact, knew that minimax approximation is a thing, however he only tried to do it by hand! It's a bit silly to try to optimize a polynomial for minimizing the maximum error and then complain when the result predictably does worse in a completely different metric. What he should have tried is taking some tool for minimax approximation, like Sollya (Remez-based, also offers other functionality), and then try too look for a minimax that would additionally satisfy his other constraints. Sollya allows a weighted approximation, so he could have tried constraining the polynomial to specifically be more accurate around zero, π/2 and π. NB: My WIP package will be even more flexible than Sollya's remez and offer more direct ways to specify such constraints, due to the completely different approach.

Why does he even care about pi/2? Shouldn't the wrapping of the input (due to the use of fixed-point) make sure that sin(pi/2) == sin(0)?

True, the author of the video should probably have done even more range reduction.

By itself, Remez is not usually helpful for these kinds of low-accuracy polynomial approximations, because getting the minmax error beneath the quantization threshold is too strong a constraint. Usually one ends up doing some brute force polishing at the end to try and reduce the polynomial degree and/or the bits necessary to represent coefficients.

FWIW, Sollya's Remez algorithm takes intermediate floating-point accuracy into account. At the very least, one should start by such a polynomial and then try to tweak it by hand, not just drag coefficients around from scratch. Numerics isn't “just messing around with numbers” as he claims; it's a large field with lots of cleverness.

As a quick example (using Maple, though, not Sollya): f(x) = x * (1.020875931 + x * (-0.0683760152 - x * 0.1122036320)) gives f(0) = 0 (exactly), f(Pi/2) = 1 (to at least nine decimal digits) and has a max error around 0.0019, significantly better than his hand-tweaked one (max error 0.007). For comparison, his fourth-order function has 0.002787 and the Bhaskara formula has about 0.0016. I don't know how many bytes it needs in MIPS assembly, but I normally count function speed in cycles :-)

If you allow a divide between two second-degree polynomials, like in the Bhaskara function, you can go down an order of magnitude in error, of course at extra cycle cost.

I looked at the video a bit more closely, and it seems he requires the symmetry around 0 (he doesn't explicitly do range reduction to positive values), so you want only odd-numbered coefficients. If so, you can do a formula of the type x * (c + x*x*(b + x*x*a)) (aka: ax⁵ + bx³ + cx; this is essentially the same as his formula with c=1 and prescaling x by a constant), which Sollya readily gives for single precision:

  > fpminimax(sin(x*2*Pi/65536.0),[|x,x^3,x^5|],[|SG...|],[0;16384]);
  Warning: For at least 2 of the constants displayed in decimal, rounding has happened.
  x * (9.58634263952262699604034423828125e-5 + x^2 * (-1.46252571143166976153082714517950080335140228271484e-13 + x^2 * 6.1585575698811928868593794069233315902067715796875e-23))

This has a maximum error of 0.000108 over that range. This is pretty close to optimal, but you can squeeze it ever so slightly lower (just below 0.0001) if you're willing to spend a lot of CPU time. :-)

How are you computing the max error here? What would be the max error for taylor series?

Take an interval, for sine it might be, for example [0,π/4], I think it was [0,π] in the video (a mistake, most probably). You're interested in the maximum approximation error over the interval, so you could, for example, sample ten thousand points from the interval uniformly, evaluate the error over all of them, and take the maximum.

Technically, what you're really looking for is the supremum, not the maximum. The former is the limit of the latter.

Sollya does something smarter than just sampling uniformly, though, it knows about the derivatives of all the functions it works with (and derivatives of the derivatives, etc), and uses that information to compute an arbitrarily small interval that is guaranteed to contain the supremum: https://www.sollya.org/sollya-weekly/help.php#supnorm

My package also does something more complicated than just sampling uniformly, however it doesn't know anything about the derivatives of the relevant functions, so the maximum it finds is just approximate.

In this case, you know the input is one out of 16384 distinct values (originally 65536, but due to range-reduction and symmetries, you only need to consider one of them). So you can simply test them all.

But in the case of the minimax polynomial, I just got it out of Maple (it can tell you as part of computing it). That's not accounting for floating-point inaccuracies, though.

Reply on news.ycombinator.com