Re: Ada.Numerics, Accuracy of trigonometric functions

[Robert Eachus <rieachus@comcast.net>]

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;