* 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 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 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 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