Question on Inverse and Hyperbolic Trigonometric functions

[Bart]

Hi,

I need someone with more knowledge of mathematics tham me to look at this.

I have a proposed patch that implements the SecH (hyperbolic secant), CscH (hyperbolic cosecant) and CotH (hyperbolic cotangent) functions.
Based upon https://en.wikipedia.org/wiki/Hyperbolic_functions#Definitions.
Hope I did not screw that one up.

Now, I'm trying to implement inverse secant (ArcSec), the inverse cosecant (ArcCsc) and inverse cotangent (ArcCot) functions for fpc's math unit.

My knowledge of trigonometry stops with the things I learned until I was 18 years old: Sin, Cos, Tan functions.

From googling, this is what I have come to understand:

ArcSec(X) = ArcCos(1/X)
ArcCsc(X) = ArcSin(1/x)
ArcCot(X) = pi/2 - ArcCot(1/X) = ArcTan(1/X), if X > 0
ArcCot(X) = 3*pi/2 - ArcTan(1/X) = pi + ArcTan(1/X), if X < 0

Since Math unit already implements ArcSin(), ArcCos() and ArcTan(), I could implement these functions as described above, if and only if my understanding is indeed correct.

It also leaves the question what the result of ArcCot(0) should be, or better, what Delphi returns in that case.


And once the above is solved, we also need implementations for
-  inverse hyperbolic secant (ArcSecH)
-  inverse hyperbolic cosecant  (ArcCscH)
-  inverse hyperbolic cotangent (ArcCotH)

Any help would be appreciated.

Bart

[Peter H]

According to Wolfram Alpha arccot(0)=pi/2.
Wolfram Alpha also shows various ways to calculate this.
https://www.wolframalpha.com/input?i=arccot%280%29

While it seems easy to get the formulas, the question is if a naive implementation is numerically stable and accurate.
There are often many ways to calculate this and all have advantages and disadvantages.
I think if stuff like this is implemented, it must be done by specialists for numerical mathematics.
Not by programmers. ;-)

E.g. ArcSec(x) can be calculated in various ways:https://www.wolframalpha.com/input?i=simplify+arcsec%28x%29

[wp]

I personally would not use these arccot formulas because they are not continuous and have that unnecessary discontinuity at x = 0.

There is another formula, arccot(x) = pi/2 - arctan(x), which appears very natural to me. See screenshot and see the similarity with arctan(x). It is used by Excel. The arccot values here have values between 0 and pi while in the "ugly" formula they have values between -pi/2 and pi/2.

Unfortunately the "ugly" formula is used by Delphi... (and Wolfram Alpha).

It's a matter of definition: The tan(x) and cot(x) are periodic functions. Since the inverse function must be unique one period must be selected for the domain of the function values. Should it run between 0 and pi, or between -pi/2 and pi/2 (or any other interval with width pi)?

BTW: The split formula ArcCot(X) = ArcTan(1/X) for X > 0 and ArcCot(X) = pi + ArcTan(1/X) for X < 0, and ArcCot(0) = pi/2 results in the simple formula again (due to the symmetry of the curve). In order to really get the discontinuous curve the X>0 branch must be used also for negative values:

Code: Pascal  [Select][+][-]

1. function ArcCot(x: Double): Double;
2. // Results are in range -pi/2 and +pi/2
3. begin
4.   if x = 0 then
5.     Result := pi/2
6.   else
7.     Result := arctan(1/x);
8. end;

In the math library of the late Wolfgang Ehrhard (https://github.com/JulStrat/MPArith) there are two functions, ArcCot as written above, and an Arccotc ("c" for "complementary", a notation which I know from the "error function", too) defined like

Code: Pascal  [Select][+][-]

1. function ArcCotC(x: Double): Double;
2. // Results are in range between 0 and pi
3. begin
4.   Result := pi/2 - arctan(x);
5. end;

This is, maybe, a good compromise.

[Thaddy]

OpenAI chatGPT gave me these:
1. inverse hyperbolic secant

Code: Pascal  [Select][+][-]

1. Function InverseSecH(x : Real):Real;
2. Begin
3.   Result := arccosh(1/x);
4. End;

2.inverse cosecant

Code: Pascal  [Select][+][-]

1. Function InverseCosecant (X : Real) : Real;
2. Var
3.   Result : Real;
4. Begin
5.   Result := ArcCos(1/X);
6.   InverseCosecant := Result;
7. End;

3.inverse hyperbolic cosecant

Code: Pascal  [Select][+][-]

1. Function Arcsch (x : Real) : Real;
2. Begin
3.   Arcsch := ln((1/x) + sqrt(1/(x*x) + 1));
4. End;

4. inverse hyperbolic cotangent

Code: Pascal  [Select][+][-]

1. function arcoth(x : Real) : Real;
2. begin
3.   arcoth := ln((x + 1) / (x - 1)) / 2;
4. end;

Model Input was:"write a Pascal function for   inverse hyperbolic cotangent " etc.
The functions seem correct but I have no idea if they are stable..
If they are correct and stable, it is proof that OpenAI is useful.
Completely clean roomed, hope the model knew what it was doing...
There are a couple I never used myself and, like Bart, been a long time ago(47 years).

What also is notable is the different programming styles the model uses,  see e.g. 1 and 2.
Sourcecode is here: https://forum.lazarus.freepascal.org/index.php/topic,61668.msg465090.html#msg465090
Some of these functions seem to be taken from standard C and translated by the model to Pascal.
The model may gie different answers for the same request when run multiple times...

[wp]

Regarding the high error rate, I would not put too much trust into code generated by ChatGPT. At least it did not consider the divisions by zero...

[Thaddy]

Regarding the high error rate, I would not put too much trust into code generated by ChatGPT. At least it did not consider the divisions by zero...

Well, that is RTL compliant  especially in the math unit there are multiple functions that do not consider division by zero. Will take days to report all those as bugs.
In that light it is a bit of a non-argument.
The rtl assumes you are smart enough to protect your code with try/except and EDivByZero.
As always: garbage in, garbage out.
Btw: I have tested these functions against their C counter parts and the results are the same.

Also, you can ask for multiple solutions, like so:

Code: Bash  [Select][+][-]

1. What is your request?
2. write 3 functions that return the inverse hyperbolic cosecant in different ways in pascal
3.  
4. Please wait for the response, this can take some time:
5.  
6.  
7. {The inverse hyperbolic cosecant (or arcsinh) is the inverse of the hyperbolic cosecant.}
8.  
9. {Function 1}
10. Function InverseHyperbolicCosecant(x : Real): Real;
11. Begin
12.   Result := Ln(1/x + Sqrt((1/x)^2 + 1));
13. End;
14.  
15. {Function 2}
16. Function InverseHyperbolicCosecant(x : Real): Real;
17. Begin
18.   Result := Ln(x + Sqrt(x^2 + 1));
19. End;
20.  
21. {Function 3}
22. Function InverseHyperbolicCosecant(x : Real): Real;
23. Begin
24.   Result := ArcTanh(1/x);
25. End;

Where did it get that comment from?

[Bart]

@Peter H
> According to Wolfram Alpha arccot(0)=pi/2.
OK.

@WP
> Unfortunately the "ugly" formula is used by Delphi... (and Wolfram Alpha).
We must follow Delphi here, not how it's implemented, but the output should be the same.
I gather that Delpi outputs values between -pi/2 and pi/2 then.


The proposal for now would then be something like this.
(Note: for Deplhi compatibilty there muste be Single, Double and Extended versions in the final implementation, which requires other trig functions to be present in those 3 flavours as well...)

Code: Pascal  [Select][+][-]

1. function ArcSec(X: float): float;
2. begin
3.   ArcSec:=ArcCos(1/X);
4. end;
5.  
6. function ArcCsc(X: float): float;
7. begin
8.   ArcCsc:=ArcSin(1/X);
9. end;
10.  
11. function ArcCot(X: float): float;
12. begin
13.   if X>0 then
14.     ArcCot:=ArcTan(1/X)
15.   else if X<0 then
16.     ArcCot:=pi + ArcTan(1/X)
17.   else
18.     ArcCot := pi/2;
19. end;
20.  
21. function ArcSecH(X : float):float;
22. begin
23.   ArcSecH:=ArcCosH(1/X);
24. end;
25.  
26. function ArcCscH(X: float): float;
27. begin
28.   ArcCscH:=ln((1/X) + sqrt(1/(x*x) + 1));
29. end;
30.  
31. function ArcCotH(X: float): float;
32. begin
33.   ArcCotH:=ln((x + 1) / (x - 1)) / 2;
34. end;

Bart

[wp]

Thaddy, don't trust this "intelligence". Putting the AI-generated code into a simple project:

Code: Pascal  [Select][+][-]

1. program Project1;
2. uses math;
3.  
4. {Function 1}
5. Function InverseHyperbolicCosecant1(x : Real): Real;
6. Begin
7.   Result := Ln(1/x + Sqrt(sqr(1/x) + 1));
8. End;
9.  
10. {Function 2}
11. Function InverseHyperbolicCosecant2(x : Real): Real;
12. Begin
13.   Result := Ln(x + Sqrt(x*x + 1));
14. End;
15.  
16. {Function 3}
17. Function InverseHyperbolicCosecant3(x : Real): Real;
18. Begin
19.   Result := ArcTanh(1/x);
20. End;
21.  
22. const
23.   w = 15;
24. var
25.   x: Double;
26. begin
27.   WriteLn('x':w, 'Function 1':w, 'Function 2':w, 'Function 3':w);
28.   x := 2;
29.   WriteLn(x:w:5, InverseHyperbolicCosecant1(x):w:5, InverseHyperbolicCosecant2(x):w:5, InverseHyperbolicCosecant3(x):w:5);
30.   x := -2;
31.   WriteLn(x:w:5, InverseHyperbolicCosecant1(x):w:5, InverseHyperbolicCosecant2(x):w:5, InverseHyperbolicCosecant3(x):w:5);
32.   x := 1.5;
33.   WriteLn(x:w:5, InverseHyperbolicCosecant1(x):w:5, InverseHyperbolicCosecant2(x):w:5, InverseHyperbolicCosecant3(x):w:5);
34.   x := -1.01;
35.   WriteLn(x:w:5, InverseHyperbolicCosecant1(x):w:5, InverseHyperbolicCosecant2(x):w:5, InverseHyperbolicCosecant3(x):w:5);
36.  
37.   ReadLn;
38.  
39. end.
40.  

Output:

              x     Function 1     Function 2     Function 3
        2.00000        0.48121        1.44364        0.54931
       -2.00000       -0.48121       -1.44364       -0.54931
        1.50000        0.62515        1.19476        0.80472
       -1.01000       -0.87436       -0.88843       -2.65165

As can be seen all three methods yield different results...

[Thaddy]

I am not trusting it... but it is a nice toy...

[Bart]

Based upons https://en.wikipedia.org/wiki/Inverse_hyperbolic_functions:

Code: Pascal  [Select][+][-]

1. function ArcSecH(X : float):float;
2. begin
3.   //instead of ArcSecH:=ArcCosH(1/X);
4.   ArcSecH:=ln((1+(sqrt(1.0-sqr(X))))/X);  
5.   //equals: ArcSecH:=ln(1+(sqrt(1-sqr(X))))-ln(X);  not sure wich one is faster
6. end;

Bart

[avra]

Now, I'm trying to implement inverse secant (ArcSec), the inverse cosecant (ArcCsc) and inverse cotangent (ArcCot) functions for fpc's math unit.

My knowledge of trigonometry stops with the things I learned until I was 18 years old: Sin, Cos, Tan functions.

If you know what to look for, that can be sufficient.

From googling, this is what I have come to understand:

ArcSec(X) = ArcCos(1/X)
ArcCsc(X) = ArcSin(1/x)

Correct.

ArcCot(X) = pi/2 - ArcCot(1/X) = ArcTan(1/X), if X > 0
ArcCot(X) = 3*pi/2 - ArcTan(1/X) = pi + ArcTan(1/X), if X < 0

Wrong.
arctan(x) + arccot(x) = PI/2 (https://www.math-only-math.com/arctan-x-plus-arccot-x-equals-pi-by-2.html), so:
arccot(x) = PI/2 - arctan(x)

I did not look at math.pas so I do not know how it handles division by zero and defined range in general, rounding, correction, infinity and quadrant finding, but fixed point math library I made 13 years ago for ATMEL AVR microcontrollers used to deal with such things via defines, so for example you could choose if you wanted to get cos(90) as 0.000000000 instead of 0.000000012 which lib normally calculates.

It even had an optional cache for last few calculated sine values so it was super fast in complex trigonometry without need to recalculate the same sine again and again, and it also had a define to optionally sacrifice 2 bits of precision for 50% speed gain in some cases.

That lib enabled me to position client's telescope and accurately track and record nice smooth videos of arbitrary sky objects with just a simple 8-bit AVR microcontroller. Standard solution before that was to use double with ARMs (double is too slow on AVR for tracking, and single does not have enough accuracy so tracking had jerky moves).

[Bart]

First of all, thanks for your input.

ArcCot(X) = pi/2 - ArcCot(1/X) = ArcTan(1/X), if X > 0
ArcCot(X) = 3*pi/2 - ArcTan(1/X) = pi + ArcTan(1/X), if X < 0

Wrong.
arctan(x) + arccot(x) = PI/2 (https://www.math-only-math.com/arctan-x-plus-arccot-x-equals-pi-by-2.html), so:
arccot(x) = PI/2 - arctan(x)

OK, I derived my solution from https://en.wikipedia.org/wiki/Inverse_trigonometric_functions, section: reciprocal arguments:

Code: Pascal  [Select][+][-]

1.   ArcCot(1/X) = pi/2 - ArcCot(X) = ArcTan(X), if X>0
2.   ArcCot(1/X) = 3/2*pi - ArcCot(X) = pi  + ArcTan(X), if X>0
3.   //substitute Y=1/X and you get the version I had

I missed the solution you give, which is on the same page.

However, as WP in reply #2 pointed out:

There is another formula, arccot(x) = pi/2 - arctan(x), which appears very natural to me. See screenshot and see the similarity with arctan(x). It is used by Excel. The arccot values here have values between 0 and pi while in the "ugly" formula they have values between -pi/2 and pi/2.

Unfortunately the "ugly" formula is used by Delphi... (and Wolfram Alpha).

Our output values must follow Delphi, or people will complain...
 

[avra]

Our output values must follow Delphi, or people will complain...

Subtraction is lighter then division (especially when there is no FPU), but ok... let it be...

[Peter H]

I think the output of Wolfram Alpha should be most reliable. Delphi or C++ are not a reference.
It does not necessarily use these symbolic formulas which are displayed, but uses the numeric and symbolic engine of mathematica.
It will probably not display digits which have a numeric error but will display "computing time exceeded".

As a test I calculated arcsec(10) and then used the numeric result to calculate arcsec^-1(result):

[Bart]

I think the output of Wolfram Alpha should be most reliable. Delphi or C++ are not a reference.

Alas, Delphi compatibility is a must, unless the behaviour of Delphi is plain wrong.
If that's the case, then the "arccot(x) = pi/2 - arctan(x)" should perhaps be used.

I'm currently not in a position to argue either of those.

As to my referal to the wiki page (section: reciprocal arguments): is that page wrong?
If not, is the formula I derived from it wrong?

Bart