8

How could I define trigonometric functions that take arguments in degrees instead of the usual radians, and compute correctly rounded results for these arguments?

Multiplying the argument by M_PI/180.0 before passing it to the corresponding function in radians does not work, because M_PI/180.0 is not π/180. Section 5.5 of the Handbook of Floating-Point Arithmetic offers a method to compute the correctly rounded product of the argument by π/180, but some arguments will still be such that this product is close to the midpoint between two consecutive representable floats, and then applying even a correctly rounded function in radians can produce the wrong final result.

Two strategies that may work alone or in combination are using higher precision and using the sinpi,cospi, tanpi trigonometric functions from CRlibm, that compute respectively sin(πx), cos(πx) and tan(πx).

For the latter strategy, there remains the problem of the division by 180, which is not exact for many arguments.

Regarding the higher-precision strategy (multiplying the argument by an extended-precision representation of π/180, then applying the extended-precision function in radians), there may remain a problem with “exact” cases. The theorem that states that the only rational results of sin, cos and tan for a rational argument are obtained in 0 only applies to the radian versions. It obviously does not apply to the degree versions, and if for some floating-point input x, sindeg(x) is exactly the midpoint between two consecutive representable floating-point numbers, then no amount of intermediate precision is enough to guarantee that the final result is correctly rounded.

Pascal Cuoq
  • 79,187
  • 7
  • 161
  • 281

4 Answers4

3

The only rationals q for which cosdeg(360q) is rational have 1, 2, 3, 4, or 6 as the denominator. This paper by Joerg Jahnel contains a short and beautiful proof using field theory in section 6. (Indeed, the author characterises the degree of the algebraic number cosdeg(360q) using Euler's totient function.) So there is no floating-point q such that cosdeg(360q) is halfway between two adjacent floating-point numbers.

So I guess the answer is "about the same way you implement sin and friends for radians," though @gnasher729 makes the excellent point that argument reduction for degrees is much, much nicer.

tmyklebu
  • 13,915
  • 3
  • 28
  • 57
  • I'll get on with reading the paper immediately, but shouldn't the result from which to infer that there are no nontrivial exact cases for floating-point `sindeg` be “there is no floating-point `q` such that sindeg(`q`) is halfway between two adjacent floating-point numbers”? – Pascal Cuoq Aug 03 '14 at 15:37
  • @PascalCuoq: Yes. The result I (mis)quoted says that the only rational `q` that lead to a rational (never mind halfway-between-two-adjacent-floating-points) value of `sindeg(q)` or `cosdeg(q)` are the multiples of 30 degrees, where the answer is always either irrational or a multiple of 1/2. This is a stronger result than what you need since every floating-point and every halfway-between-floating-points is rational. – tmyklebu Aug 03 '14 at 15:51
  • Actually I appreciate the elementary proof in section 3 much more (though it only covers sindeg and cosdeg). – Pascal Cuoq Aug 03 '14 at 23:19
  • 2
    @PascalCuoq: If you're after a similar result for tan, http://oberlin.edu/faculty/jcalcut/tanpap.pdf looks good. – Mark Dickinson Aug 06 '14 at 19:20
2

It's difficult. On the positive side, you can reduce the argument to +/- 45 degrees exactly. So you need correctly rounded results between +/- 45 degrees. For very small x, sin (x) is about x * (pi / 180) which is hard enough to get rounded exactly.

To get mostly correctly rounded results for the sine function for example, take -45 <= x <= 45. Split x into xhi = round (512 x) / 512 and xlo = x - xhi. Let sin (x degrees) ≈ ax - bx^3. Round a and b so that s (x) a*xhi - b * (xhi^3) is calculated exactly. Calculate the remainder sin (x degrees) - s (x) carefully; the rounding error should be quite small because the result is small. Add to s (x), this will most of the time give the correctly rounded result.

gnasher729
  • 51,477
  • 5
  • 75
  • 98
1

Well, this is a difficult question. Let me clarify some points:

  • What precision is required for the the output? Is it IEEE 754 single or double precision or non standard? Moreover, I assume, that input, i.e. the one represented in degrees, should represented in the same precision as outputs, as this is the case for the normal radian inputs.
  • What is your performance metrics? CRlibm is optimized to produce correctly rounded double precision results. On the other side, MPFR is used for arbitrary precision but it is much slower than CRlibm when you need only double precision output.
  • What is your working range? i.e. [min argument, max argmunet]? This matters for CRlibm as it works for double precision ranges. However, it won't really matter for MPFR.

I basically recommend to use MPFR if you are obliged to use inputs only in degrees. Let me remind you that any argument in degrees, when it is multiplied by (Pi/180), it produces a transcendental number. However, What is passed to the trigonometric function is the floating point representation rounded, preferably rounded to nearest integer, to the working precision.

I recommend you to do as follows:

  1. Use MPFR, use the C library whenever possible as it offers much better performance than its bindings.
  2. Set the MPFR precision much higher than your target precision. For example (target precision + 300). By doing this, you avoid any loss of accuracy for the operation ((Argument*Pi)/180). This can be done easily in MPFR C library by mpfr_set_default_prec().
  3. Do the operation: X_n=(Argument*Pi)/180, and then execute Sin(X_n) or whatever function you want. There is a constant Pi in MPFR, which is represented within your working precision
  4. Round your result to the target precision.

"Elementary functions", by Muller shows statistically that most, NOT ALL, of the hard cases are correctly rounded if the working precision is slightly larger than twice the target precision. But in your case, as the input is theoretically transcendental, to be safe, at the expense of the performance, make the working precision much higher than the target. Actually 10x is totally sufficient for almost 100% of cases, if you require up to double precision final result.

If you need a low precision, i.e. single precision or lower, it is feasible to do exhaustive test to decide on the lowest working precision which produces all cases correctly rounded.

Garp
  • 109
  • 5
  • 1
    1- “Set the MPFR precision much higher than your target precision” does NOT work if there are representable floating-point arguments the image of which is exactly a midpoint between two representable consecutive floating-point values. I'm sorry if I didn't make that clear in my question, but this is what “there may remain a problem with “exact” cases” meant. It has already been pointed out that one need not worry about such exact cases for sin, cos and tan, however. – Pascal Cuoq Aug 03 '14 at 22:46
  • 1
    2- “By doing this, you avoid any loss of accuracy for the operation ((Argument*Pi)/180)” MPFR being binary, there will be loss of accuracy. For instance, tandeg(90) gives a large but finite value with your proposal, rounded to some infinity in the target precision. The correct result is NaN. – Pascal Cuoq Aug 03 '14 at 22:48
  • The output of a transcendental function, e.g sin, cos, log, can NOT be exactly a midpoint. The real output must be a transcendental number. – Garp Aug 03 '14 at 22:48
  • Yes, that was part of my question, and tmyklebu has answered it. That sindeg, cosdeg and tandeg didn't have any exact cases was not obvious before his answer, however. Note that I am interest in sindeg, cosdeg and tandeg. Theorems about sin, cos, tan do not apply directly. – Pascal Cuoq Aug 03 '14 at 22:50
  • Yes, it would give a large finite number as the output didnot exceed the maximum value. I remind you that if you use MPFR, it sets very high maximum values by default. Therefore I assume you didnot exceed the maximum. – Garp Aug 03 '14 at 22:51
  • So you are saying that tandeg(90) would produce a large finite number with your method. Subsequently rounded to an infinity (either positive or negative depending on the width of the significand used in MPFR). And that does not strike you as incorrect? – Pascal Cuoq Aug 03 '14 at 22:54
  • If you want the special outputs to be done right. e.g. tan(90) you should deal with it manually, as tan(90) is a special case just like cos(90) or sin(0). the real 90 is a transcendental which gives a special output when it passed to tan() or cos() – Garp Aug 03 '14 at 22:54
  • Oh, so tandeg **does** have exact cases to take care of before launching an extended-precision computation, then. How do you know that it does not have more exact cases that should be treated specially? – Pascal Cuoq Aug 03 '14 at 22:56
  • Yes, tandeg(90) even for round to inf would give the next above representation which is finite also in MPFR. try to set the maximum exponents to a lower value. They are set to incredibly large values by default. – Garp Aug 03 '14 at 22:57
  • This is was proved in the 19th century actually, you can find it in Elementary functions by Muller. I will get you the page. – Garp Aug 03 '14 at 22:58
  • Page 19, it shows that the transcendental function of a ration argument, which is the case for floating point representations, except for very special cases, such as cos(90) or sin(0), it should produce a transcendental number. – Garp Aug 03 '14 at 23:02
  • And you are talking about implementing sindeg from sin, when I already said in my question that it was not obviously possible at all. – Pascal Cuoq Aug 03 '14 at 23:02
  • No am talking on sindeg, cosdeg, ... I just want you to make the operation (deg*Pi)/180 in much higher precision than the target to avoid a possible loss of accuracy. – Garp Aug 03 '14 at 23:04
  • You should first convert it do radians, then do a correctly rounded operation. – Garp Aug 03 '14 at 23:06
  • And my question clearly states that this approach **does not obviously work at ANY higher precision** because of “exact cases”. So if you want to answer this, you should explain why exact cases are not a problem. You haven't. Fortunately, tmyklebu has: you should read his answer. – Pascal Cuoq Aug 03 '14 at 23:07
  • You should handle these cases yourself, this is not a problem. The real problem is to get otherwise correctly rounded results, that is by passing it to MPFR – Garp Aug 03 '14 at 23:14
  • “You should handle these cases yourself” But what are these cases? Your answer doesn't say. Actually, your answer claims that there aren't any, which is obviously false. – Pascal Cuoq Aug 03 '14 at 23:17
  • I didn't claim such thing, floating point representations are finite, thus there is no really 90 for example. But you cannot obtain correctly rounded results without converting to radians. As MPFR uses polynomials applied only for real numbers. – Garp Aug 03 '14 at 23:25
  • I am pretty sure that 90 can be represented as a binary64 floating-point number (and thus can be passed as argument to a floating-point function, such as the function `tandeg` considered in this question). – Pascal Cuoq Aug 03 '14 at 23:26
  • 90 can be represented by round-to-nearest(pi/2). Anyways, I advise you to lower your dynamic range to get the inf cases for example. Any other rational output or zero output can be obtained if you are working initially in much larger precision, and when you finally round the result to the target precision it will be your correct result. – Garp Aug 03 '14 at 23:32
  • “90 can be represented by round-to-nearest(pi/2).” No, it can't. Doing that will result in an infinite result, either +inf or -inf depending on the direction of the approximation. The **correct** result for tandeg(90) is not infinity, it is NaN. Anyway, I just realized that tmyklebu's answer only proved the absence of exact cases for sindeg and cosdeg, so if you have a convincing argument for the absence of exact cases for tandeg, you should put that in your answer. – Pascal Cuoq Aug 03 '14 at 23:35
  • I mean it can represented in FP finitely by round-to-nearest(pi/2). Anyways, I am sorry for wasting your time. I'll delete the answer if you want to. – Garp Aug 03 '14 at 23:39
  • Right, let me say it once again (this is my last comment on this question): your answer suggests to compute `tandeg(f)` by computing, in extended precision, tan(f * π / 180) and then rounding to the target precision. **This does not work if tandeg(`f`) is a midpoint between two consecutive target floating-point numbers, or an otherwise exceptional value such as `NaN`.** Therefore, the question is whether tandeg(`f`) can be such an exact midpoint, and for what values of `f` it is one. This difficulty is spelt out in the question and you have not addressed it. – Pascal Cuoq Aug 04 '14 at 11:43
  • You are right, I am sorry. However, the paper in tmyklebu says that it is only the multiples of (180/{1, 2, 3, 4, or 6}) that should be considered. Thus you can handle these cases by yourself, and what is left which is the vast majority is fed to MPFR in a much higher working precision then a final rounding to your target precision. – Garp Aug 04 '14 at 13:09
  • Yes. Also as far as I can see the article only handles sine and cosine, but there is a hint at http://math.stackexchange.com/questions/380940/tangent-function-of-rational-angle for proving the same property for the tangent. – Pascal Cuoq Aug 04 '14 at 13:13
1

You first need to detect the exact cases, and this has already been answered. Now, for the other cases, there's the well-known problem of the table maker's dilemma. If your arithmetic has a fixed (and small) precision and you want a certified bound on the intermediate precision that might be needed, there are two known solutions:

  • Obtain a bound based on Nesterenko and Waldschmidt's theorem, as described in Section 4.3 of my PhD thesis (BTW, I think this would also give you the form of the exact cases). But you will get very large bounds on the precision (at least several millions bits?).
  • Find the hardest-to-round case. It suffices to do the search in [0,180], since any larger argument will reduce to a value in [0,180] with the same fractional part (since the period is an integer).
vinc17
  • 2,829
  • 17
  • 23