* Ada.Numerics, Accuracy of trigonometric functions
@ 2016-10-07 10:13 Markus Schöpflin
2016-10-07 10:23 ` Markus Schöpflin
` (4 more replies)
0 siblings, 5 replies; 10+ messages in thread
From: Markus Schöpflin @ 2016-10-07 10:13 UTC (permalink / raw)
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
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Ada.Numerics, Accuracy of trigonometric functions
2016-10-07 10:13 Ada.Numerics, Accuracy of trigonometric functions Markus Schöpflin
@ 2016-10-07 10:23 ` Markus Schöpflin
2016-10-07 10:23 ` Brian Drummond
` (3 subsequent siblings)
4 siblings, 0 replies; 10+ messages in thread
From: Markus Schöpflin @ 2016-10-07 10:23 UTC (permalink / raw)
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
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Ada.Numerics, Accuracy of trigonometric functions
2016-10-07 10:13 Ada.Numerics, Accuracy of trigonometric functions Markus Schöpflin
2016-10-07 10:23 ` Markus Schöpflin
@ 2016-10-07 10:23 ` Brian Drummond
2016-10-07 10:40 ` Markus Schöpflin
2016-10-07 14:52 ` Dennis Lee Bieber
` (2 subsequent siblings)
4 siblings, 1 reply; 10+ messages in thread
From: Brian Drummond @ 2016-10-07 10:23 UTC (permalink / raw)
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
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Ada.Numerics, Accuracy of trigonometric functions
2016-10-07 10:23 ` Brian Drummond
@ 2016-10-07 10:40 ` Markus Schöpflin
0 siblings, 0 replies; 10+ messages in thread
From: Markus Schöpflin @ 2016-10-07 10:40 UTC (permalink / raw)
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
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Ada.Numerics, Accuracy of trigonometric functions
2016-10-07 10:13 Ada.Numerics, Accuracy of trigonometric functions Markus Schöpflin
2016-10-07 10:23 ` Markus Schöpflin
2016-10-07 10:23 ` Brian Drummond
@ 2016-10-07 14:52 ` Dennis Lee Bieber
2016-10-07 15:19 ` Markus Schöpflin
2016-10-07 16:38 ` Jeffrey R. Carter
2016-10-22 22:38 ` Robert Eachus
4 siblings, 1 reply; 10+ messages in thread
From: Dennis Lee Bieber @ 2016-10-07 14:52 UTC (permalink / raw)
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/
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Ada.Numerics, Accuracy of trigonometric functions
2016-10-07 14:52 ` Dennis Lee Bieber
@ 2016-10-07 15:19 ` Markus Schöpflin
2016-10-07 22:15 ` Dennis Lee Bieber
0 siblings, 1 reply; 10+ messages in thread
From: Markus Schöpflin @ 2016-10-07 15:19 UTC (permalink / raw)
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
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Ada.Numerics, Accuracy of trigonometric functions
2016-10-07 10:13 Ada.Numerics, Accuracy of trigonometric functions Markus Schöpflin
` (2 preceding siblings ...)
2016-10-07 14:52 ` Dennis Lee Bieber
@ 2016-10-07 16:38 ` Jeffrey R. Carter
2016-10-10 7:44 ` Markus Schöpflin
2016-10-22 22:38 ` Robert Eachus
4 siblings, 1 reply; 10+ messages in thread
From: Jeffrey R. Carter @ 2016-10-07 16:38 UTC (permalink / raw)
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
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Ada.Numerics, Accuracy of trigonometric functions
2016-10-07 15:19 ` Markus Schöpflin
@ 2016-10-07 22:15 ` Dennis Lee Bieber
0 siblings, 0 replies; 10+ messages in thread
From: Dennis Lee Bieber @ 2016-10-07 22:15 UTC (permalink / raw)
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/
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Ada.Numerics, Accuracy of trigonometric functions
2016-10-07 16:38 ` Jeffrey R. Carter
@ 2016-10-10 7:44 ` Markus Schöpflin
0 siblings, 0 replies; 10+ messages in thread
From: Markus Schöpflin @ 2016-10-10 7:44 UTC (permalink / raw)
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
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Ada.Numerics, Accuracy of trigonometric functions
2016-10-07 10:13 Ada.Numerics, Accuracy of trigonometric functions Markus Schöpflin
` (3 preceding siblings ...)
2016-10-07 16:38 ` Jeffrey R. Carter
@ 2016-10-22 22:38 ` Robert Eachus
4 siblings, 0 replies; 10+ messages in thread
From: Robert Eachus @ 2016-10-22 22:38 UTC (permalink / raw)
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;
^ permalink raw reply [flat|nested] 10+ messages in thread
end of thread, other threads:[~2016-10-22 22:38 UTC | newest]
Thread overview: 10+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2016-10-07 10:13 Ada.Numerics, Accuracy of trigonometric functions Markus Schöpflin
2016-10-07 10:23 ` Markus Schöpflin
2016-10-07 10:23 ` Brian Drummond
2016-10-07 10:40 ` Markus Schöpflin
2016-10-07 14:52 ` Dennis Lee Bieber
2016-10-07 15:19 ` Markus Schöpflin
2016-10-07 22:15 ` Dennis Lee Bieber
2016-10-07 16:38 ` Jeffrey R. Carter
2016-10-10 7:44 ` Markus Schöpflin
2016-10-22 22:38 ` Robert Eachus
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox