7 ms·

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

by nsajko 3y ago

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 https://gitlab.com/nsajko/FindMinimaxPolynomial.jl.

lifthrasiir 3y ago
> 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 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 https://arxiv.org/pdf/2104.04043.pdf
- nsajko 3y ago
  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 https://news.ycombinator.com/item?id=36201197
  - lifthrasiir 3y ago
    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?
    - nsajko 3y ago
      > 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.
sosodev 3y ago
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.
- nsajko 3y ago
  I'm not claiming anything extraordinary, though, most or all of what I said is well known in the field.
  - sosodev 3y ago
    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.
    - nsajko 3y ago
      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.
      - Sesse__ 3y ago
        > 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…
tyg13 3y ago
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.
- nsajko 3y ago
  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.
  - Sesse__ 3y ago
    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)?
    - nsajko 3y ago
      True, the author of the video should probably have done even more range reduction.
thebooktocome 3y ago
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.
- Sesse__ 3y ago
  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.
  - boywitharupee 3y ago
    How are you computing the max error here? What would be the max error for taylor series?
    - nsajko 3y ago
      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 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.

by nsajko 3y ago

lifthrasiir 3y ago

> 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 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 https://arxiv.org/pdf/2104.04043.pdf

nsajko 3y ago
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 https://news.ycombinator.com/item?id=36201197
- lifthrasiir 3y ago
  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?
  - nsajko 3y ago
    > 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.

sosodev 3y ago

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.

nsajko 3y ago
I'm not claiming anything extraordinary, though, most or all of what I said is well known in the field.
- sosodev 3y ago
  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.
  - nsajko 3y ago
    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.
    - Sesse__ 3y ago
      > 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…

tyg13 3y ago

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.

nsajko 3y ago
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.
- Sesse__ 3y ago
  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)?
  - nsajko 3y ago
    True, the author of the video should probably have done even more range reduction.

thebooktocome 3y ago

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.

Sesse__ 3y ago
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.
- boywitharupee 3y ago
  How are you computing the max error here? What would be the max error for taylor series?
  - nsajko 3y ago
    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 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.