### Bookstore

 Computer Math and Games in Pascal (preview) Lazarus Handbook

### Author Topic: Strange error with arccos()  (Read 3719 times)

#### Paolo

• Sr. Member
• Posts: 333
##### Re: Strange error with arccos()
« Reply #45 on: August 09, 2022, 02:47:38 pm »
Done some testing, some improvements are observed even if very little.

#### Josh

• Hero Member
• Posts: 1024
##### Re: Strange error with arccos()
« Reply #46 on: August 09, 2022, 03:23:08 pm »
Hi

I think it needs to be fixed, as the adding pi and subtracting 2pi, is causing a round error.

compact code; which is equally as fast but highly doubt if would be accepted
Code: Pascal  [Select][+][-]
1. function JoshArcTan2(Y,X:Float):Float;
2. function SignGreaterEqualZero(Const AValue:Float):Integer;inline;
3. begin
4.   If AValue<0.0 Then SignGreaterEqualZero:=-1 Else SignGreaterEqualZero:=+1;
5. end;
6. begin
7.   If X>0.0 Then JoshArcTan2:=Arctan(Y/X)
8.   Else
9.     If X<0.0 then JoshArcTan2:=Arctan(Y/X)+SignGreaterEqualZero(Y)*pi
10.     else JoshArcTan2:=Sign(Y)*pi/2;
11. end;
12.

or just correct existing code to stop the rounding issue

Code: Pascal  [Select][+][-]
1. function arctan2(y,x : float) : float;
2. begin
3.   if (x=0) then
4.     begin
5.       if y=0 then
6.         arctan2:=0.0
7.       else if y>0 then
8.         arctan2:=pi/2
9.       else if y<0 then
10.         arctan2:=-pi/2;
11.     end
12.   else
13.   if x<0.0 then
14.     begin
15.       if y<0.0 then
16.         arctan2:=arctan(y/x)-pi
17.       else
18.         arctan2:= arctan(y/x)+pi;
19.     end
20.   else arctan2:=ArcTan(y/x);
21. end;
22.

I started investigating FPC tanh function, on initial test it is flawed, as the range check is not correct,
try calling it with a value of -356, you get a FT error, according to the range check it is in range..
Not too sure about the formula used, if mem serves it should be using cosh,sinh for calculation, but there is probably other equations to solve it, I just have a gut feeling something is wrong; but at the moment i do not have a good time slot to check it out.
« Last Edit: August 09, 2022, 03:27:29 pm by Josh »
The best way to get accurate information on the forum is to post something wrong and wait for corrections.

#### Paolo

• Sr. Member
• Posts: 333
##### Re: Strange error with arccos()
« Reply #47 on: August 09, 2022, 03:43:18 pm »
I'll submit the second version, for the same reasons you exposed. I'll check tanh !
« Last Edit: August 09, 2022, 03:47:21 pm by Paolo »

#### Paolo

• Sr. Member
• Posts: 333
##### Re: Strange error with arccos()
« Reply #48 on: August 09, 2022, 04:43:20 pm »
for tanh the suspect is:

Code: Pascal  [Select][+][-]
1. x:=-359;
2. ...
3. y:=sinh(x)/cosh(x);   //works as expected
4. y:=(Exp(x)-Exp(-x))/(Exp(X)+Exp(-x); //works as above
5. y:=Tanh(x); //crash !
6. y:=(1-exp(-2x))/(1+Exp(-2*x)); //that is could be as implemented tanh and crash again,
7. //and you see way : -2*X is > 0 and Exp itself crash ! too large value
8.
9.

if the hypotesis above are right and the computation of tanh over extended range
should be (using the tanh(-x)=-tanh(x), x > 0)

Code: Pascal  [Select][+][-]
1. funcyion NewTanh(x:double):double;
2. if X < 0 then
3.   result:=-(1-exp(2x))/(1+Exp(2*x)); //ie -tanh(-x) (here X < 0, tanh(x)=-tanh(-x)
4. else
5.   result:=(1-exp(-2x))/(1+Exp(-2*x)); //ie tanh(x)
6.

now working for x:=-359;

do you agree ?
« Last Edit: August 09, 2022, 04:49:47 pm by Paolo »

#### Paolo

• Sr. Member
• Posts: 333
##### Re: Strange error with arccos()
« Reply #49 on: August 10, 2022, 12:42:41 pm »
but in any case Tanh(x) for x approaching 18 already gives 1, so tanh with X above such values is = 1, non need to really compute it.

#### Josh

• Hero Member
• Posts: 1024
##### Re: Strange error with arccos()
« Reply #50 on: August 10, 2022, 02:18:58 pm »
hi

i agree, when going to 50 digits, the min/max is only 60 before it is +-1/

below code that test varius routine, note joshtanh is fast and just as accurate and going in the range +-60(so has an advantage)

Code: Pascal  [Select][+][-]
1. program Project1;
2. uses sysutils,math,crt;
3. const n=20000000;
4.       numberofroutines=4;
5.
6. Type TTanHFunc = function(const x: Float): Float;
7.
8. var i:longint;
9.     v:float;
10.     SampleData,Differences:array[1..N] of Float;
11.     results:array[1..n,1..numberofroutines] of Float;
12.     verify:boolean=true;
13.     BaseTime:double;
14.     TimeInRoutine:Array[1..numberofroutines] of double;
15.     RoutineNames:Array[1..numberofroutines] of String;
16.     mc:integer=0;
17.
18. const  MaxTanh = 5678.22249441322; // Ln(MaxExtended)/2
19.        TestMaxRange=350;
20.
21. function mathtanh(const x : float) : float;
22. var Temp : float;
23. begin
24.   if x>MaxTanh then exit(1.0)
25.   else if x<-MaxTanh then exit (-1.0);
26.   temp:=exp(-2*x);
27.   mathtanh:=(1-temp)/(1+temp)
28. end;
29.
30.
31. function NewTanhIdea1(const x:float):float;
32. var SignX:Integer;
33. begin
34.   SignX:=sign(x);
35.   if SignX=0 then result:=0
36.   else result:=SignX*(
37.                      (1-exp(((SignX*-1)*2)*x))
38.                           /
39.                      (1+exp(((SignX*-1)*2)*x)));
40. end;
41.
42. function NewTanh2(const x:float):float;
43. var SignX:Integer;
44. begin
45.   if X < 0 then
46.   result:=-(1-exp(2*x))/(1+Exp(2*x)) //ie -tanh(-x) (here X < 0, tanh(x)=-tanh(-x)
47. else
48.   result:=(1-exp(-2*x))/(1+Exp(-2*x)); //ie tanh(x)
49. end;
50.
51.
52.
53.
54. function JoshNewTanh(const X: float): float;
55. const xrange=60;
56. begin
57.   if X>xrange then Result := 1.0
58.   else
59.     if X<-xrange then Result := -1.0
60.     else
61.       begin
62.         Result := Exp(X);
63.         Result := Result*Result;
64.         Result := (Result-1.0)/(Result+1.0);
65.       end;
66. end;
67.
68. var t:tdatetime;
69.     d:float;
70.     tol:float=0.00000000000001;
71.     dx:integer=48;
72.     dy:integer=1;
73.
74. procedure PrepareData;
75. var j,k:longint;
76. begin
77.   WriteLn('Preparing Dataset');
78.   for j:=1 to n do
79.   begin
80.     SampleData[j]:=random-0.5*(2*TestMaxRange);
81.     for k:=1 to numberofroutines do results[j,k]:=0;
82.   end;
83.   for k:=1 to numberofroutines do differences[k]:=0;
84. end;
85.
86. Function Test(Descr: String; TanHFunc: TTanHFunc; r:integer):Double;
87. var
88.   t: TDateTime;
89.   i: Integer;
90.   ans:float;
91. begin
92.   Write('  Testing ' + Descr + '...');
93.   t:=Now;
94.   for i := 1 to N do
95.   begin
96.     ans:=TanHFunc(SampleData[i]);
97.     results[i,r]:=ans;
98.     if verify then if abs(results[i,r]-results[i,1])>tol then differences[r]:=differences[r]+abs(results[i,r]-results[i,1]);
99.   end;
100.   t := Now - t;
101.   if r=1 then BaseTime:=t;
102.   WriteLn(' done: ', FormatDateTime('s.zzz', t), ' sec');
103.   if verify then writeLn('     Differences :',differences[r]:8:8);
104.   WriteLn;
105.   TimeInRoutine[r]:=TimeInRoutine[r]+t;
106.   result:=t;
107. end;
108.
109.
110.
111. begin
112.   RoutineNames[1]:='MathTanH....';
113.   RoutineNames[2]:='NewTanh2....';
114.   RoutineNames[3]:='NewTanhIdea1';
115.   RoutineNames[4]:='JoshNewTanh.';
116.   randomize;
117.   randseed:=0;
118.   mc:=0;
119.   repeat
120.     inc(mc);
121.     PrepareData;
122.     t:=test(RoutineNames[1],@mathtanh,1);
123.     t:=test(RoutineNames[2],@NewTanh2,2);
124.     t:=test(RoutineNames[3],@NewTanhIdea1,3);
125.     t:=test(RoutineNames[4],@JoshNewTanh,4);
126.     TextBackground(blue);
127.     gotoxy(dx,dy);write('        Total Times   Count '+mc.ToString); clreol;
128.     gotoxy(dx,dy+1);write('        hr mn se ms     Avg');  clreol;
129.     if mc>0 then
130.     begin
131.       for i:=1 to numberofroutines do
132.       begin;
133.         gotoxy(dx,dy+1+i);
134.         write(' '+RoutineNames[i]+':', FormatDateTime('hh.nn.ss.zzz', TimeInRoutine[i]), '  ',FormatDateTime('s.zzz',TimeInRoutine[i]/mc), ' sec'); clreol;
135.       end;
136.     end;
137.     TextBackground(black);
138.     gotoxy(1,1);
139.   until false;
140. end.
141.

What u think?
« Last Edit: August 10, 2022, 03:36:37 pm by Josh »
The best way to get accurate information on the forum is to post something wrong and wait for corrections.

#### wp

• Hero Member
• Posts: 10291
##### Re: Strange error with arccos()
« Reply #51 on: August 10, 2022, 03:49:08 pm »
Code: Pascal  [Select][+][-]
1. procedure PrepareData;
2. [...]
3.     SampleData[j]:=random-0.5*(2*TestMaxRange);
4.
What's this? Why don't you cancel the 2? Or do you mean this, with brackets around "random - 0.5"?
Code: Pascal  [Select][+][-]
1.     SampleData[j]:=(random-0.5)*(2*TestMaxRange);
2.

#### Paolo

• Sr. Member
• Posts: 333
##### Re: Strange error with arccos()
« Reply #52 on: August 10, 2022, 04:07:00 pm »
I like newtanh2, because in any case it computes exp function in the right direction (ie toward 0) so no overflow can occur and it should work for any x without check input value range
« Last Edit: August 10, 2022, 05:08:45 pm by Paolo »

#### Josh

• Hero Member
• Posts: 1024
##### Re: Strange error with arccos()
« Reply #53 on: August 10, 2022, 05:43:06 pm »
@wp
yes it should be
Code: [Select]
`SampleData[j]:=(random-0.5)*(2*TestMaxRange)`Was not thinking straight, initially i was checking in small range of values to verify results across routines. then later expanded to larger values.

@Paolo
I like the newtan2 function because it does not burp on large numbers,

The fact the the Math routine checks for a +- range of 5678.22249441322, which itself cannot handle and causes a crash is a problem and needs to be fixed, a range check of +-26 should be adequate for a precision of 22 decimal places, which i think is the accuracy of the extended type.

The best way to get accurate information on the forum is to post something wrong and wait for corrections.

#### Paolo

• Sr. Member
• Posts: 333
##### Re: Strange error with arccos()
« Reply #54 on: August 10, 2022, 06:21:56 pm »
My suspect is the fact is that the code below

Quote
Const MaxTanh = 5678.22249441322; // Ln(MaxExtended)/2

function tanh(x : float) : float;
var Temp : float;
begin
if x>MaxTanh then
exit(1.0)
else
if x<-MaxTanh then
exit (-1.0);
temp:=exp(-2*x);
tanh:=(1-temp)/(1+temp)
end;

is a problem for:
1- it does not consider correctly the large negaive numbers (as we aleady discussed)
2- as in the comment "//Ln(MaxExtended)/2" for win 64 there is not "Extended" (it is an alias of double), so the max imput is wrongly assumed, not checked ma it should be

const MAXTANH = somevalue;  //ln(Maxdouble)/2, something just less 365

#### Paolo

• Sr. Member
• Posts: 333
##### Re: Strange error with arccos()
« Reply #55 on: August 10, 2022, 06:40:02 pm »

#### wp

• Hero Member
• Posts: 10291
##### Re: Strange error with arccos()
« Reply #56 on: August 10, 2022, 06:47:27 pm »
for the arctan2

https://gitlab.com/freepascal.org/fpc/source/-/issues/39861
and fixed 20 minutes later - that's the power of Open Source.

#### Paolo

• Sr. Member
• Posts: 333
##### Re: Strange error with arccos()
« Reply #57 on: August 10, 2022, 07:00:19 pm »
@wp fixed even before I copy paste here the link, unbeliveable  !!!!! faster than the arctan

#### Paolo

• Sr. Member
• Posts: 333
##### Re: Strange error with arccos()
« Reply #58 on: August 18, 2022, 12:06:46 am »
For who is interested on tanh issue topic

https://gitlab.com/freepascal.org/fpc/source/-/issues/39867