Ada.Numerics, Accuracy of trigonometric functions

Ada.Numerics, Accuracy of trigonometric functions

[Markus Schöpflin]

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

Re: Ada.Numerics, Accuracy of trigonometric functions

[Markus Schöpflin]

What I forgot to add is that the results are the same using SHORT_FLOAT, 
LONG_FLOAT, or LONG_LONG_FLOAT when calling cos(x) without the cycle parameter.

But when adding the cycle parameter, LONG_FLOAT and LONG_LONG_FLOAT yield 
correct results, only SHORT_FLOAT is still wrong.

Markus

Re: Ada.Numerics, Accuracy of trigonometric functions

[Brian Drummond]

On Fri, 07 Oct 2016 12:13:48 +0200, Markus Schöpflin wrote:

> 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.
> 
> ---%<---
> 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:
> 
> ---%<---
> 
> 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.

C tends to silently coerce "float" to "double" any time you pass a float 
argument to a function unless you jump through hoops to make such 
coercion impossible. 

So my suspicion is simply that the C program is giving you double 
precision.

-- Brian

Re: Ada.Numerics, Accuracy of trigonometric functions

[Markus Schöpflin]

Am 07.10.2016 um 12:23 schrieb Brian Drummond:

[...]

> C tends to silently coerce "float" to "double" any time you pass a float
> argument to a function unless you jump through hoops to make such
> coercion impossible.
>
> So my suspicion is simply that the C program is giving you double
> precision.

I can see in the generated assembler dump that the libc cosf function is 
called. So I don't think this is the case.

And I get the same result in Ada using Long_Float or Long_Long_Float; see my 
follow-up.

Markus

Re: Ada.Numerics, Accuracy of trigonometric functions

[Dennis Lee Bieber]

On Fri, 7 Oct 2016 12:13:48 +0200, Markus Schöpflin <no.spam@spam.spam>
declaimed the following:


>X = 4294967296.00000; COS(X) =  1.00000      <-- WHAT?

	That kind of surprises me...  Not the cos() result, but the input X

	Short float should only have around 23-24 bits of significance in the
mantissa -- enough for around 7 decimal digits. Yet you are seeing the full
UNSIGNED 32-bit integer value with a stack of fractional zeros tacked on.

	I'm tempted to think the compiler, not the run-time, is computing the
2^32 as a universal (or whatever the preferred term is)  constant using the
smallest data type that holds the full value.
-- 
	Wulfraed                 Dennis Lee Bieber         AF6VN
    wlfraed@ix.netcom.com    HTTP://wlfraed.home.netcom.com/

Re: Ada.Numerics, Accuracy of trigonometric functions

[Markus Schöpflin]

Am 07.10.2016 um 16:52 schrieb Dennis Lee Bieber:
> On Fri, 7 Oct 2016 12:13:48 +0200, Markus Schöpflin <no.spam@spam.spam>
> declaimed the following:
>
>
>> X = 4294967296.00000; COS(X) =  1.00000      <-- WHAT?
>
> 	That kind of surprises me...  Not the cos() result, but the input X
>
> 	Short float should only have around 23-24 bits of significance in the
> mantissa -- enough for around 7 decimal digits. Yet you are seeing the full
> UNSIGNED 32-bit integer value with a stack of fractional zeros tacked on.
>
> 	I'm tempted to think the compiler, not the run-time, is computing the
> 2^32 as a universal (or whatever the preferred term is)  constant using the
> smallest data type that holds the full value.
>

2**32 has an exact representation as short float: 0x4f800000 (in binary this 
is 0 10011111 00000000000000000000000; sign 0, exponent 159, mantissa 0).

Markus

Re: Ada.Numerics, Accuracy of trigonometric functions

[Jeffrey R. Carter]

On 10/07/2016 03:13 AM, Markus Schöpflin wrote:
>
> Size of SHORT_FLOAT =  32
> Maximum relative error of COS =  2.38419E-07
> Angle threshold =  12
> X = 4294967296.00000; COS(X) =  1.00000      <-- WHAT?

For GNAT (I have 4.9.3) Ada.Numerics.Generic_Elementary_Functions is defined in 
terms of Ada.Numerics.Aux, which defines

type Double is new Long_Long_Float;

Cos, after some checks for zero, converts its argument to Double and calls 
Aux.Cos. Aux.Cos, for abs X > Pi/4, reduces abs X and then invokes the FPU 
instruction fcos or fsin on the reduced value, depending on the quadrant.

Reduce is defined as

    procedure Reduce (X : in out Double; Q : out Natural);
    --  Implements reduction of X by Pi/2. Q is the quadrant of the final
    --  result in the range 0 .. 3. The absolute value of X is at most Pi.

and contains the following

       K  : Double := X * Two_Over_Pi;
    begin
       --  For X < 2.0**32, all products below are computed exactly.
       --  Due to cancellation effects all subtractions are exact as well.
       --  As no double extended floating-point number has more than 75
       --  zeros after the binary point, the result will be the correctly
       --  rounded result of X - K * (Pi / 2.0).

Since X = 2**32, this doesn't apply. It looks to me as if Reduce is converting X 
to zero, which is then passed to fcos, giving 1.

This result is independent of the actual type used to instantiate 
Generic_Elementary_Functions.

The C cosf function probably just passes its argument to fcos.

Investigation of what GNAT does for Cos with Cycle is left as an exercise for 
the reader.

-- 
Jeff Carter
"It's symbolic of his struggle against reality."
Monty Python's Life of Brian
78

Re: Ada.Numerics, Accuracy of trigonometric functions

[Dennis Lee Bieber]

On Fri, 7 Oct 2016 17:19:09 +0200, Markus Schöpflin <no.spam@spam.spam>
declaimed the following:

>Am 07.10.2016 um 16:52 schrieb Dennis Lee Bieber:
>> On Fri, 7 Oct 2016 12:13:48 +0200, Markus Schöpflin <no.spam@spam.spam>
>> declaimed the following:
>>
>>
>>> X = 4294967296.00000; COS(X) =  1.00000      <-- WHAT?
>>
>> 	That kind of surprises me...  Not the cos() result, but the input X
>>
>> 	Short float should only have around 23-24 bits of significance in the
>> mantissa -- enough for around 7 decimal digits. Yet you are seeing the full
>> UNSIGNED 32-bit integer value with a stack of fractional zeros tacked on.
>>
>> 	I'm tempted to think the compiler, not the run-time, is computing the
>> 2^32 as a universal (or whatever the preferred term is)  constant using the
>> smallest data type that holds the full value.
>>
>
>2**32 has an exact representation as short float: 0x4f800000 (in binary this 
>is 0 10011111 00000000000000000000000; sign 0, exponent 159, mantissa 0).
>
>Markus

	Ah... One of those magic numbers then...
-- 
	Wulfraed                 Dennis Lee Bieber         AF6VN
    wlfraed@ix.netcom.com    HTTP://wlfraed.home.netcom.com/

Re: Ada.Numerics, Accuracy of trigonometric functions

[Markus Schöpflin]

Am 07.10.2016 um 18:38 schrieb Jeffrey R. Carter:

> It looks to me as if Reduce is converting X to zero, which is then passed
> to fcos, giving 1.
>
> This result is independent of the actual type used to instantiate
> Generic_Elementary_Functions.

This explains a lot, then.

Thank you very much,
Markus

Re: Ada.Numerics, Accuracy of trigonometric functions

[Robert Eachus]

On Friday, October 7, 2016 at 6:13:50 AM UTC-4, Markus Schöpflin wrote:

> 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):
 
> 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?

Yes, you should report this as an error to AdaCore.  The correct fix is to raise Constraint_Error, and mention why in the documentation.   I'd like to raise Argument_Error, but even mentioning that possibility may result in lots more posts than the issue is worth.

I tried your program both with Short and Long_Float, and values of X up to 2**32 with and without scaling.  The unscaled numbers worked up perfectly until 2**32 (see below).  There seems to be some rounding somewhere in the scaled version.  It's fine up to about 2**20.  The unscaled numbers are accurate to 14 digits up until 2**32:
Size of Long_FLOAT =  64
Maximum relative error of COS =  4.44089209850063E-16
Angle threshold =  26
X**I =          1.00; COS(X**I) =  0.54030230586814; Scaled Cos = 0.54030230586814
X**I =          2.00; COS(X**I) = -0.41614683654714; Scaled Cos =-0.41614683654714
X**I =          4.00; COS(X**I) = -0.65364362086361; Scaled Cos =-0.65364362086361
X**I =          8.00; COS(X**I) = -0.14550003380861; Scaled Cos =-0.14550003380861
X**I =         16.00; COS(X**I) = -0.95765948032338; Scaled Cos =-0.95765948032338
X**I =         32.00; COS(X**I) =  0.83422336050651; Scaled Cos = 0.83422336050651
X**I =         64.00; COS(X**I) =  0.39185723042955; Scaled Cos = 0.39185723042955
X**I =        128.00; COS(X**I) = -0.69289582192017; Scaled Cos =-0.69289582192017
X**I =        256.00; COS(X**I) = -0.03979075993116; Scaled Cos =-0.03979075993116
X**I =        512.00; COS(X**I) = -0.99683339084820; Scaled Cos =-0.99683339084820
X**I =       1024.00; COS(X**I) =  0.98735361821985; Scaled Cos = 0.98735361821985
X**I =       2048.00; COS(X**I) =  0.94973433482365; Scaled Cos = 0.94973433482365
X**I =       4096.00; COS(X**I) =  0.80399061348585; Scaled Cos = 0.80399061348585
X**I =       8192.00; COS(X**I) =  0.29280181314670; Scaled Cos = 0.29280181314670
X**I =      16384.00; COS(X**I) = -0.82853419643601; Scaled Cos =-0.82853419643601
X**I =      32768.00; COS(X**I) =  0.37293782932771; Scaled Cos = 0.37293782932771
X**I =      65536.00; COS(X**I) = -0.72183475091266; Scaled Cos =-0.72183475091266
X**I =     131072.00; COS(X**I) =  0.04209081525030; Scaled Cos = 0.04209081525030
X**I =     262144.00; COS(X**I) = -0.99645672654313; Scaled Cos =-0.99645672654313
X**I =     524288.00; COS(X**I) =  0.98585201574610; Scaled Cos = 0.98585201574611
X**I =    1048576.00; COS(X**I) =  0.94380839390131; Scaled Cos = 0.94380839390132
X**I =    2097152.00; COS(X**I) =  0.78154856879715; Scaled Cos = 0.78154856879720
X**I =    4194304.00; COS(X**I) =  0.22163633077774; Scaled Cos = 0.22163633077789
X**I =    8388608.00; COS(X**I) = -0.90175467375876; Scaled Cos =-0.90175467375863
X**I =   16777216.00; COS(X**I) =  0.62632298329153; Scaled Cos = 0.62632298329196
X**I =   33554432.00; COS(X**I) = -0.21543904120159; Scaled Cos =-0.21543904120051
X**I =   67108864.00; COS(X**I) = -0.90717203905228; Scaled Cos =-0.90717203905196
X**I =  134217728.00; COS(X**I) =  0.64592221687655; Scaled Cos = 0.64592221687210
X**I =  268435456.00; COS(X**I) = -0.16556897949058; Scaled Cos =-0.16556897950206
X**I =  536870912.00; COS(X**I) = -0.94517382606090; Scaled Cos =-0.94517382605945
X**I = 1073741824.00; COS(X**I) =  0.78670712294119; Scaled Cos = 0.78670712290611
X**I = 2147483648.00; COS(X**I) =  0.23781619457280; Scaled Cos = 0.23781619446243
X**I = 4294967296.00; COS(X**I) =  1.00000000000000; Scaled Cos =-0.88688691530282
[2016-10-22 18:28:01] process terminated successfully, elapsed time: 00.34s

with ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS; 

with ADA.TEXT_IO; use ADA.TEXT_IO; 
with ADA.Long_FLOAT_TEXT_IO; use ADA.Long_FLOAT_TEXT_IO; 

procedure TEST_TRIG is 

    MRE_COS : constant := 2.0 * Long_FLOAT'MODEL_EPSILON; 
    THRESHOLD : constant := Long_FLOAT'MACHINE_MANTISSA / 2; 
   Pi : constant := Ada.Numerics.Pi;
   Two_Pi : constant := Pi + Pi;
    package EF is new ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS (Long_FLOAT); 

   X : constant Long_FLOAT := 2.0;
   Temp: Integer;
    Y : constant Long_FLOAT := EF.COS (X); 

begin 

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

    for I in 0..32 loop
      PUT ("X**I = ");        PUT (X**I, Fore => 10, Aft => 2, EXP => 0); 
      PUT ("; COS(X**I) = "); PUT (EF.Cos(X**I), EXP => 0);
      Temp := Integer(X**I/(2.0*Pi));
      Put ("; Scaled Cos =");
      PUT (EF.Cos(Long_Float(Long_Long_Float(X**I)
            - Long_Long_Float'(Long_Long_Float(Temp)*Two_Pi))), EXP => 0);
      NEW_LINE;
      
    end loop;


end TEST_TRIG;