DegreesToRadians and RadiansToDegrees overflow and inaccuracy

[ed0d2b2ce19451f2]

The following methods have issues with overflow and inaccuracy:

• Single.DegreesToRadians(Single)
• Single.RadiansToDegrees(Single)
• Double.DegreesToRadians(Double)
• Double.RadiansToDegrees(Double)

The current implementations for Single appear to be equivalent to:

float DegreesToRadians(float degrees) {
    return (degrees * Pi) / 180.0f;
}
float RadiansToDegrees(float radians) {
    return (radians * 180.0f) / Pi;
}

Both return infinity when the multiplication yields infinity, even if the end result would otherwise be finite. The Double versions have the same issue.

Overflow aside, these implementations seem unnecessarily inaccurate.

Here is a table comparing the existing and two alternative binary32 implementations of each conversion for the 2^23 = 8388608 values in [1, 2) or any other binade not involving overflow or subnormals. In the implementation column, operators and constants are exact and f(...) indicates rounding to the nearest value representable in the floating point format.

|             |                   |      distance from correctly rounded      |
|             |                   |     2 ULP     |      1 ULP      |  0 ULP  |
|  direction  | implementation    | count |   %   |  count  |   %   |         |
| ----------- | ----------------- | ----- | ----- | ------- | ----- | ------- |
| deg --> rad | f(f(x*f(pi))/180) | 20895 | 0.249 | 3147446 | 37.52 | balance |
| deg --> rad | f(x*f(pi/180))    |     0 | 0     |  792445 |  9.44 | balance |
| deg --> rad | f(x/f(180/pi))    |     0 | 0     | 1194111 | 14.23 | balance |
| rad --> deg | f(f(x*180)/f(pi)) |  3502 | 0.041 | 3183331 | 37.94 | balance |
| rad --> deg | f(x/f(pi/180))    |     0 | 0     |  792441 |  9.44 | balance |
| rad --> deg | f(x*f(180/pi))    |     0 | 0     | 1194116 | 14.23 | balance |

Here is a similar table comparing binary64 implementations for 1e7 samples chosen uniformly at random from [1, 2).

|             |                   | distance from correctly rounded |
|             |                   |   2 ULP   |   1 ULP   |  0 ULP  |
|  direction  | implementation    |     %     |     %     |         |
| ----------- | ----------------- | --------- | --------- | ------- |
| deg --> rad | f(f(x*f(pi))/180) | 0.0001357 | 33.20333  | balance |
| deg --> rad | f(x*f(pi/180))    | 0         | 11.05166  | balance |
| deg --> rad | f(x/f(180/pi))    | 0         | 22.69309  | balance |
| rad --> deg | f(f(x*180)/f(pi)) | 0         | 33.23503  | balance |
| rad --> deg | f(x/f(pi/180))    | 0         | 11.04806  | balance |
| rad --> deg | f(x*f(180/pi))    | 0         | 22.70538  | balance |

For both Single and Double, multiplying (for DegreesToRadians) or dividing (for RadiansToDegrees) by the nearest floating point value to pi/180 is the most accurate of the options compared, and the current implementation is the least accurate. The more accurate implementations also avoid the overflow issue.

Of course it's possible to improve accuracy even further, but going beyond this level will likely require sacrificing performance.

[dotnet-policy-service]

Tagging subscribers to this area: @JulieLeeMSFT, @jakobbotsch
See info in area-owners.md if you want to be subscribed.

[dotnet-policy-service]

Tagging subscribers to this area: @dotnet/area-system-numerics
See info in area-owners.md if you want to be subscribed.

[tannergooding]

To start, this is documented with reasoning as to why here: https://github.com/dotnet/dotnet/blob/main/src/runtime/src/libraries/System.Private.CoreLib/src/System/Numerics/ITrigonometricFunctions.cs#L65-L111

This reasoning is particularly true when you consider that the generic math interfaces can be implemented by any numeric type and so can have any potential number of significand bits. This then factors into being consistent across implementations in light of that and us wanting to get the API shipped in the simplest form at first.

In particular what is best (between using the current algorithm, an in-type computed value from pi and 180, or a direct nearest representable constant to the infinitely precise pi / 180 or 180 / pi value) is essentially random and dependent on P (significant bits). It's just "luck of the draw" and there isn't strictly a way to account for that that I can discern.

Double (binary64, P=53) and Single (binary32, P=24) happen to be such special cases that get "lucky", where the in-type computed pi / 180 still lands on the nearest representable value to the infinitely precise pi / 180 and where that nearest representable is itself unusually close (0.145 and 0.170 of a half-ulp, against an average of 0.5). Quad (binary128, P=113) is an "unlucky" case for DegreesToRadians and where the current algorithm is the better one instead.

That luck is also direction dependent. Single gets it for DegreesToRadians but not for RadiansToDegrees, Quad gets it the other way around, and only Double gets both.

The current fallback algorithm is a reasonable one that balances precision across all potential inputs and where overflow is largely considered an edge case we somewhat need to accept regardless due to perf and complexity here.

---

It's potentially too late for .NET 11 to change anything here, but I will think more on if deviating for float/double is something we're wanting to do and will do a bar check for .NET 11.

There is definitely a few simple improvements, particularly when FMA is available that can improve accuracy and perf for the concrete types (particularly float/double), but its ultimately not somethign we can guarantee in all scenarios.

[jeffhandley]

This would have my support for a fix in RC1