Ada.Numerics, Accuracy of trigonometric functions

["Markus Schöpflin" <no.spam@spam.spam>]

The following tests have been performed using GNAT 7.4 on a 64bit Linux system.

I'm puzzled by the seemingly bad accuracy of calling COS on a large value of 
type Short_Float.

Given the following procedure which calculates y1 = cos(2**32) and y2 = 
cos(2**32 mod 2*PI):

---%<---
with ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS;

with ADA.TEXT_IO; use ADA.TEXT_IO;
with ADA.SHORT_FLOAT_TEXT_IO; use ADA.SHORT_FLOAT_TEXT_IO;

procedure TEST_TRIG is

    MRE_COS : constant := 2.0 * SHORT_FLOAT'MODEL_EPSILON;
    THRESHOLD : constant := SHORT_FLOAT'MACHINE_MANTISSA / 2;

    package EF is new ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS (SHORT_FLOAT);

    X : constant SHORT_FLOAT := 2.0**32;
    Y : constant SHORT_FLOAT := EF.COS (X);

begin

    PUT_LINE ("Size of SHORT_FLOAT = " & POSITIVE'IMAGE (SHORT_FLOAT'SIZE));
    PUT_LINE ("Maximum relative error of COS = " & SHORT_FLOAT'IMAGE (MRE_COS));
    PUT_LINE ("Angle threshold = " & POSITIVE'IMAGE (THRESHOLD));

    PUT ("X = ");        PUT (X, EXP => 0);
    PUT ("; COS(X) = "); PUT (Y, EXP => 0);
    NEW_LINE;

end TEST_TRIG;
--->%---

The (for me surprising result is):

---%<---
Size of SHORT_FLOAT =  32
Maximum relative error of COS =  2.38419E-07
Angle threshold =  12
X = 4294967296.00000; COS(X) =  1.00000      <-- WHAT?
--->%---

The corresponding C program:

---%<---
#include <math.h>
#include <stdio.h>

int main() {

   float x = exp2f(32.0);

   printf("sizeof(float) = %lu\n", sizeof(float));
   printf("x = %f; cos(x) = %f\n", x, cosf(x));

   return 0;
}
--->%---

gives:

---%<---
sizeof(float) = 4
x = 4294967296.000000; cos(x) = -0.886887    <--- OK
--->%---

Now, why is the error on cos(2**32) so large in this case? I am aware that I'm 
way beyond the required angle threshold mentioned in G.2.4(12), so it's up to 
the compiler to decide on the accuracy.

Looking at the implementation defined characteristics for GNAT, I find that 
strict mode is on by default, so G.2.4 should apply. For G.2.4(10) the 
documentation states (for the value of the angle threshold and the accuracy 
beyond the threshold): "Information on this subject is not yet available.".

So in principle the compiler is free to give me any answer once I'm beyond the 
minimum required angle threshold, which is about 12 for short float.

But if I understand G.2.4(10) correctly, cos(x1, cycle=>2*PI) should give me 
an answer within the maximum relative error, but this is not the case either, 
as EF.COS (X, CYCLE => TWO_PI) gives 0.33575.

Now for the questions: Is my reasoning above correct? Why does using the cycle 
parameter not help? Do other Ada compilers give the same bad result in this 
case? Should I report this as an error to AdaCore?

TIA,
Markus