programming the x87 for math.h etc.

[NickJohnson]

Hi again,

I just got around to trying to implement the trig functions, and I was wondering what the "limits specified by the standard" are exactly. Do things just have to be within FLT_EPSILON or DBL_EPSILON, and is it as strict for large numbers? I can't get my Taylor series approximation of sin() to get within DBL_EPSILON of glibc's sin() for inputs greater than about 100.0, regardless of the number of terms (due to inaccuracy from fmod() I think) but it's close.

[Candy]

Can't you just clip for sin(), tan() and cos()?

On a sidenote, ask Solar. He's been working on a functional implementation too.

[NickJohnson]

Yes, that's what I'm doing. I move everything to [-pi, pi], record if it's in [-pi, 0] (because sin(-x) = -sin(x)), move it to [0, pi], then flip everything from [pi/2, pi] to [0, pi/2] (because sin(x) = sin(pi-x)). My Taylor series is built around 0 currently, but I suppose one around pi/4 would be faster for [0, pi/2]. Unfortunately, I think there is inaccuracy introduced by fmod() or by glibc's value for M_PI and M_PI_2, which is why the error increases beyond DBL_EPSILON for large values of x. I'm attempting to write this in C with doubles, so maybe there would be a benefit to using the 80 bit semi-long-doubles that the x87 reportedly uses.

[jal]

I move everything to [-pi, pi], record if it's in [-pi, 0] (because sin(-x) = -sin(x)), move it to [0, pi], then flip everything from [pi/2, pi] to [0, pi/2] (because sin(x) = sin(pi-x)). My Taylor series is built around 0 currently, but I suppose one around pi/4 would be faster for [0, pi/2]. Unfortunately, I think there is inaccuracy introduced by fmod() or by glibc's value for M_PI and M_PI_2, which is why the error increases beyond DBL_EPSILON for large values of x.

Haha, lol, now I know what it feels like for someone not into computers when I talk with my friend about OSes :).


JAL

[Combuster]

What you need to realize, is that for larger numbers precision is lost relative to the domain of sin. While for numbers near 1, you have 22 significant bits, at 128 its 15 significant bits (assuming float). Near values of 32 million, there are no significant bits left: 32 million mod pi - 32 million + 1 mod pi is two; that means that output values may vary over a range of more than one depending on the bits that are lost in rounding.

That means it's theoretically not possible to guarantee precision of 0.000001 (whatever the value of epsilon is) at that magnitude. Instead, it is a more sensible test that the output of any function is correct for any input number from x-epsilon/2 to x+epsilon/2.

Try out what the results are for 100 float, the first number after 100, and the first number before 100 for both glibc and your version, and what the results are when you cast the numbers to long double and compute the sine for the number boundaries, and see if your sine implementation still conforms in that condition.

[Owen]

You're conflating significant bits and fractional bits; significanty bits are constant over the entire normalised number range (Obviously, denormals are excepted)

Otherwise, your points stand.

[NickJohnson]

Does that mean that it is also valid to say that if my function matches glibc's everywhere in [0, 2*pi] to within the type's epsilon, then it is working? I just checked, and it does. The fmod() I'm using is from glibc, so any rounding/precision errors outside of [-2*pi, 2*pi] wouldn't be due to my function.

[gerryg400]

Yes, that's what I'm doing. I move everything to [-pi, pi], record if it's in [-pi, 0] (because sin(-x) = -sin(x)), move it to [0, pi], then flip everything from [pi/2, pi] to [0, pi/2] (because sin(x) = sin(pi-x)). My Taylor series is built around 0 currently, but I suppose one around pi/4 would be faster for [0, pi/2]. Unfortunately, I think there is inaccuracy introduced by fmod() or by glibc's value for M_PI and M_PI_2, which is why the error increases beyond DBL_EPSILON for large values of x. I'm attempting to write this in C with doubles, so maybe there would be a benefit to using the 80 bit semi-long-doubles that the x87 reportedly uses.

You're conflating significant bits and fractional bits; significanty bits are constant over the entire normalised number range (Obviously, denormals are excepted)

Hmm, I thought this was an English only forum...

[Combuster]

You're conflating significant bits and fractional bits; significanty bits are constant over the entire normalised number range (Obviously, denormals are excepted)

I accuse thee of the same. For large numbers fewer bits from the float's fraction cause changes in [-pi, pi], and since only those bits determine the precision sin can reach, they are therefore significant bits. Larger number = less significant bits, the same fractional bits (the fractional part of a float is constant size, except for denormals).

[Solar]

Sorry, I've been missing this thread for some time it seems.

Using the FPU to implement <math.h> is (as I said before) somewhat tricky. I admit I haven't coded for the FPU so far (and haven't started on <math.h> for PDCLib), so I don't really know how it works, but there are special cases to consider, especially when NaN and INF values are involved.

For example, speaking of sin(), when fed an INF value, it's expected to return NaN and raise an "invalid" floating point exception. A naive FPU implementation might miss that. (Then again, it might not. As I said, no hands-on experience so far.)

You're also expected to set errno as appropriate, which involves reading the FPU status, which involves some heavy penalty on FPU performance.

All in all, I believe Plaugher when he writes that it's not worth using the FPU to implement <math.h>, and will go the integer arithmetics path. That way, I won't be bitten by the next FPU bug Intel comes up with, either.

As for what kind of precision the standard demands, I have to admit I couldn't find the passage ad-hoc. (As I haven't ever used <math.h> myself, and haven't started on implementing it myself, my knowledge about that corner of the standard is flakey.)