4 ms·

This has been the standard algorithm used by every libm for decades. Its not special to Musl.

enriquto 5y ago
But isn't this code rarely called in practice? I guess on intel architectures the compiler just calls the fsin instruction of the cpu.
- adgjlsfhk1 5y ago
  No. The fsin instruction is inaccurate enough to be useless. It gives 0 correct digits when the output is close to 0.
  - enriquto 5y ago
    > 0 correct digits when the output is close to 0 this is an amusing way to describe the precision of sub-normal floating point numbers
    - lifthrasiir 5y ago
      It is much more amusing if you describe it in ulps; for some inputs the error can reach > 2^90 ulps, more than the mantissa size itself.
    - adgjlsfhk1 5y ago
      It's not just sub-normal numbers. As https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/ https://randomascii.wordpress.com/2014/10/09/intel-underesti... shows, fsin only uses 66 bits of pi, which means you have roughly no precision whenever abs(sin(x))<10^-16 which is way bigger than the biggest subnormal (10^-307 or so)
      - ramshorns 5y ago
        In that range, just returning x would be way better. Maybe even perfect actually - if x is less than 10^-16, then the error of x^3/6 is less than the machine precision for x.
        adgjlsfhk1 5y ago
        the error isn't when x is small. it's when sin(x) is small. the problem happens for x near multiples of pi
- stephencanon 5y ago
  No. FSIN has accuracy issues as sibling mentions, but is also much slower than a good software implementation (it varies with uArch, but 1 result every ~100 cycles is common; even mediocre scalar software implementations can produce a result every twenty cycles).
- chrisseaton 5y ago
  > I guess on intel architectures the compiler just calls the fsin instruction of the cpu. Do people do that in practice? It's on the FPU, which is basically legacy emulated these days, and it's inaccurate.
  - enriquto 5y ago
    > the FPU, which is basically legacy you'll pry my long doubles from my cold, dead hands!
    - jhgb 5y ago
      The question is if you wouldn't be better served with double-doubles today. You get ~100 bits of mantissa AND you can still vectorize your computations.
      - enriquto 5y ago
        Sure. There should be a gcc flag to make "long double" become quadruple precision. The thing is, my first programming language was x86 assembler and the fpu was the funniest part. Spent weeks as a teenager writing almost pure 8087 code. I have a lot of emotional investment in that tiny rolling stack of extended precision floats.
        colejohnson66 5y ago
        How is "quad floating point" implemented on x86? Is it software emulation?
    - chrisseaton 5y ago
      You're selectively quoting me - I said it's 'legacy emulated'. It's emulated using very long-form microcode, so basically emulated in software. I didn't say it was simply 'legacy'.
      - eyelidlessness 5y ago
- jhgb 5y ago
  FSIN only works on x87 registers which you will rarely use on AMD64 systems -- you really want to use at least scalar SSE2 today (since that is whence you receive your inputs as per typical AMD64 calling conventions anyway). Moving data from SSE registers to the FP stack just to calculate FSIN and then moving it back to SSE will probably kill your performance even if your FSIN implementation is good. If you're vectorizing your computation over 4 double floats or 8 single floats in an AVX register, it gets even worse for FSIN.
  - stephencanon 5y ago
    Moving between x87 and xmm registers is actually fairly cheap (it's through memory, so it's not free, but it's also not _that_ bad). FSIN itself is catastrophically slow.
    - deleted 5y ago
      [deleted]
    - jhgb 5y ago
      Fair enough, and I imagine there may even be some forwarding going on? There often is when a load follows a store, if I remember correctly. (Of course this will be microarchitecture-dependent.)
    - deleted 5y ago
      [deleted]
- lifthrasiir 5y ago
  x87 trig functions can be very inaccurate due to the Intel's original sloppy implementation and subsequent compatibility requirements [1]. [1] http://nighthacks.com/jag/blog/134/index.html http://nighthacks.com/jag/blog/134/index.html
- jcelerier 5y ago
  Wasn't there some blog article a few years ago which showed how glibc's implementation was faster than fsin ?

But isn't this code rarely called in practice? I guess on intel architectures the compiler just calls the fsin instruction of the cpu.

adgjlsfhk1 5y ago
No. The fsin instruction is inaccurate enough to be useless. It gives 0 correct digits when the output is close to 0.
- enriquto 5y ago
  > 0 correct digits when the output is close to 0 this is an amusing way to describe the precision of sub-normal floating point numbers
  - lifthrasiir 5y ago
    It is much more amusing if you describe it in ulps; for some inputs the error can reach > 2^90 ulps, more than the mantissa size itself.
  - adgjlsfhk1 5y ago
    It's not just sub-normal numbers. As https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/ https://randomascii.wordpress.com/2014/10/09/intel-underesti... shows, fsin only uses 66 bits of pi, which means you have roughly no precision whenever abs(sin(x))<10^-16 which is way bigger than the biggest subnormal (10^-307 or so)
    - ramshorns 5y ago
      In that range, just returning x would be way better. Maybe even perfect actually - if x is less than 10^-16, then the error of x^3/6 is less than the machine precision for x.
      - adgjlsfhk1 5y ago
        the error isn't when x is small. it's when sin(x) is small. the problem happens for x near multiples of pi
stephencanon 5y ago
No. FSIN has accuracy issues as sibling mentions, but is also much slower than a good software implementation (it varies with uArch, but 1 result every ~100 cycles is common; even mediocre scalar software implementations can produce a result every twenty cycles).
chrisseaton 5y ago
> I guess on intel architectures the compiler just calls the fsin instruction of the cpu. Do people do that in practice? It's on the FPU, which is basically legacy emulated these days, and it's inaccurate.
- enriquto 5y ago
  > the FPU, which is basically legacy you'll pry my long doubles from my cold, dead hands!
  - jhgb 5y ago
    The question is if you wouldn't be better served with double-doubles today. You get ~100 bits of mantissa AND you can still vectorize your computations.
    - enriquto 5y ago
      Sure. There should be a gcc flag to make "long double" become quadruple precision. The thing is, my first programming language was x86 assembler and the fpu was the funniest part. Spent weeks as a teenager writing almost pure 8087 code. I have a lot of emotional investment in that tiny rolling stack of extended precision floats.
      - colejohnson66 5y ago
        How is "quad floating point" implemented on x86? Is it software emulation?
  - chrisseaton 5y ago
    You're selectively quoting me - I said it's 'legacy emulated'. It's emulated using very long-form microcode, so basically emulated in software. I didn't say it was simply 'legacy'.
    - eyelidlessness 5y ago
jhgb 5y ago
FSIN only works on x87 registers which you will rarely use on AMD64 systems -- you really want to use at least scalar SSE2 today (since that is whence you receive your inputs as per typical AMD64 calling conventions anyway). Moving data from SSE registers to the FP stack just to calculate FSIN and then moving it back to SSE will probably kill your performance even if your FSIN implementation is good. If you're vectorizing your computation over 4 double floats or 8 single floats in an AVX register, it gets even worse for FSIN.
- stephencanon 5y ago
  Moving between x87 and xmm registers is actually fairly cheap (it's through memory, so it's not free, but it's also not _that_ bad). FSIN itself is catastrophically slow.
  - deleted 5y ago
    [deleted]
  - jhgb 5y ago
    Fair enough, and I imagine there may even be some forwarding going on? There often is when a load follows a store, if I remember correctly. (Of course this will be microarchitecture-dependent.)
  - deleted 5y ago
    [deleted]
lifthrasiir 5y ago
x87 trig functions can be very inaccurate due to the Intel's original sloppy implementation and subsequent compatibility requirements [1]. [1] http://nighthacks.com/jag/blog/134/index.html http://nighthacks.com/jag/blog/134/index.html
jcelerier 5y ago
Wasn't there some blog article a few years ago which showed how glibc's implementation was faster than fsin ?