Wrong formulae for complex elementary functions

Hi!

The man pages for cacos, cacosf, cacosl, catan, catanf, catanl, cacosh, cacoshf, cacoshl, catanh, catanhf, and catanhl contain wrong maths.

cacos, cacosf, cacosl:
 The formula given in the man page
   cacos(z) = -i clog(z + csqrt(z * z - 1))
 gives wrong results in second and fourth quadrant of complex plain. The formula
   cacos(z) = -i clog(z + I*csqrt(1 - z * z))
 gives correct results.

catan, catanf, catanl:
 The formula given in the man page
   catan(z) = 1 / 2i clog((1 + iz) / (1 - iz))
 gives wrong results on the negative imaginary axis beginning at -I (along one of the two branch cuts). Besides, the formula is written in an ambiguous way.
 The formula
   catan(z) = (clog(1 + iz) - clog(1 - iz)) / 2i
 gives correct results.

cacosh, cacoshf, cacoshl:
 The formula given in the man page
   cacosh(z) = (0.5) * clog((1 + z) / (1 - z))
 gives wrong results everywhere in the complex plain. (The formula seems to be copied from the one for catanh, where it is sometimes correct.)
 The formula
   cacosh(z) = 2 * clog(csqrt((z + 1)/2) + csqrt((z - 1)/2)) gives correct results.

catanh, catanhf, catanhl:
 The formula given in the man page
   catanh(z) = 0.5 * clog((1 + z) / (1 - z))
 gives wrong results on the positive real axis beginning at 1 (along one of the two branch cuts).
 The formula
   catanh(z) = 0.5 * (clog(1 + z) - clog(1 - z))
 gives correct results.

I've also checked casin, casinf, casinl, casinh, casinhf, and casinhl and the formulae given there
   casin(z) = -i clog(iz + csqrt(1 - z * z))
   casinh(z) = clog(z + csqrt(z * z + 1))
are actually correct.

I suspect that some of these errors are due to somebody trying to "simplify" these formulae. Since this can ruin the behavior, I recommend to add a comment to the upstream sources that warns against attempt of simplification.

For the sake of reference, I am looking at manpages-dev 3.25-1 on Debian/squeeze.

Bye!
 -richy.
--
Richard B. Kreckel
<http://www.ginac.de/~kreckel/>

Hi Andries,

Since you are the mathematician, can you comment?

Thanks,

Michael

Quoted 61 lines · Richard B. Kreckel

Hi!

The man pages for cacos, cacosf, cacosl, catan, catanf, catanl, cacosh, cacoshf, cacoshl, catanh, catanhf, and catanhl contain wrong maths.

cacos, cacosf, cacosl:
The formula given in the man page

    cacos(z) = -i clog(z + csqrt(z * z - 1))
  gives wrong results in second and fourth quadrant of complex plain.
  The formula
    cacos(z) = -i clog(z + I*csqrt(1 - z * z))

gives correct results.

catan, catanf, catanl:
The formula given in the man page

    catan(z) = 1 / 2i clog((1 + iz) / (1 - iz))
  gives wrong results on the negative imaginary axis beginning at -I
  (along one of the two branch cuts). Besides, the formula is written
  in an ambiguous way.
  The formula
    catan(z) = (clog(1 + iz) - clog(1 - iz)) / 2i

gives correct results.

cacosh, cacoshf, cacoshl:
The formula given in the man page

    cacosh(z) = (0.5) * clog((1 + z) / (1 - z))
  gives wrong results everywhere in the complex plain. (The formula
  seems to be copied from the one for catanh, where it is sometimes
  correct.)
  The formula
    cacosh(z) = 2 * clog(csqrt((z + 1)/2) + csqrt((z - 1)/2))

gives correct results.

catanh, catanhf, catanhl:
The formula given in the man page

    catanh(z) = 0.5 * clog((1 + z) / (1 - z))
  gives wrong results on the positive real axis beginning at 1 (along
  one of the two branch cuts).
  The formula
    catanh(z) = 0.5 * (clog(1 + z) - clog(1 - z))

gives correct results.

I've also checked casin, casinf, casinl, casinh, casinhf, and casinhl and the formulae given there

    casin(z) = -i clog(iz + csqrt(1 - z * z))
    casinh(z) = clog(z + csqrt(z * z + 1))

are actually correct.

I suspect that some of these errors are due to somebody trying to "simplify" these formulae. Since this can ruin the behavior, I recommend to add a comment to the upstream sources that warns against attempt of simplification.

For the sake of reference, I am looking at manpages-dev 3.25-1 on Debian/squeeze.

Bye!
-richy.
--
Richard B. Kreckel
<http://www.ginac.de/~kreckel/>

Hi everybody,

Quoted 69 lines · Michael Kerrisk

Hi Andries,

Since you are the mathematician, can you comment?

Thanks,

Michael

Quoted 58 lines · Richard B. Kreckel

Hi!

The man pages for cacos, cacosf, cacosl, catan, catanf, catanl, cacosh, cacoshf, cacoshl, catanh, catanhf, and catanhl contain wrong maths.

cacos, cacosf, cacosl:
 The formula given in the man page
   cacos(z) = -i clog(z + csqrt(z * z - 1))
 gives wrong results in second and fourth quadrant of complex plain. The formula
   cacos(z) = -i clog(z + I*csqrt(1 - z * z))
 gives correct results.

catan, catanf, catanl:
 The formula given in the man page
   catan(z) = 1 / 2i clog((1 + iz) / (1 - iz))
 gives wrong results on the negative imaginary axis beginning at -I (along one of the two branch cuts). Besides, the formula is written in an ambiguous way.
 The formula
   catan(z) = (clog(1 + iz) - clog(1 - iz)) / 2i
 gives correct results.

cacosh, cacoshf, cacoshl:
 The formula given in the man page
   cacosh(z) = (0.5) * clog((1 + z) / (1 - z))
 gives wrong results everywhere in the complex plain. (The formula seems to be copied from the one for catanh, where it is sometimes correct.)
 The formula
   cacosh(z) = 2 * clog(csqrt((z + 1)/2) + csqrt((z - 1)/2)) gives correct results.

catanh, catanhf, catanhl:
 The formula given in the man page
   catanh(z) = 0.5 * clog((1 + z) / (1 - z))
 gives wrong results on the positive real axis beginning at 1 (along one of the two branch cuts).
 The formula
   catanh(z) = 0.5 * (clog(1 + z) - clog(1 - z))
 gives correct results.

I've also checked casin, casinf, casinl, casinh, casinhf, and casinhl and the formulae given there
   casin(z) = -i clog(iz + csqrt(1 - z * z))
   casinh(z) = clog(z + csqrt(z * z + 1))
are actually correct.

I suspect that some of these errors are due to somebody trying to "simplify"
these formulae. Since this can ruin the behavior, I recommend to add a comment to the upstream sources that warns against attempt of simplification.

For the sake of reference, I am looking at manpages-dev 3.25-1 on Debian/squeeze.

Quoted 5 lines

Quoted 5 lines

Bye!
 -richy.
--
Richard B. Kreckel
<http://www.ginac.de/~kreckel/>

This doesn't have seem to have been fixed in the Linux man pages, does it?

May I suggest a non-mathematical approach to moving forward?

Just write a little C program and compare the complex results of
a) the library implementation of catan, cacos, etc.,
b) the formula written in the man page, and
c) the formula I proposed above.
If you include a couple of values from all four quadrants, you'll see that my formula always agrees with the implementation while the one from the man page doesn't.

Is there anything else I can do to convince you that this should be fixed? If yes, pretty please, let me know! I would really appreciate getting this fixed.

The above is what you meant, right?

Hello Richard,

Quoted 72 lines · Richard B. Kreckel

Hi everybody,

Quoted 69 lines · Michael Kerrisk

Hi Andries,

Since you are the mathematician, can you comment?

Thanks,

Michael

Quoted 58 lines · Richard B. Kreckel

Hi!

The man pages for cacos, cacosf, cacosl, catan, catanf, catanl, cacosh, cacoshf, cacoshl, catanh, catanhf, and catanhl contain wrong maths.

cacos, cacosf, cacosl:
 The formula given in the man page
   cacos(z) = -i clog(z + csqrt(z * z - 1))
 gives wrong results in second and fourth quadrant of complex plain. The formula
   cacos(z) = -i clog(z + I*csqrt(1 - z * z))
 gives correct results.

catan, catanf, catanl:
 The formula given in the man page
   catan(z) = 1 / 2i clog((1 + iz) / (1 - iz))
 gives wrong results on the negative imaginary axis beginning at -I (along one of the two branch cuts). Besides, the formula is written in an ambiguous way.
 The formula
   catan(z) = (clog(1 + iz) - clog(1 - iz)) / 2i
 gives correct results.

cacosh, cacoshf, cacoshl:
 The formula given in the man page
   cacosh(z) = (0.5) * clog((1 + z) / (1 - z))
 gives wrong results everywhere in the complex plain. (The formula seems to be copied from the one for catanh, where it is sometimes correct.)
 The formula
   cacosh(z) = 2 * clog(csqrt((z + 1)/2) + csqrt((z - 1)/2)) gives correct results.

catanh, catanhf, catanhl:
 The formula given in the man page
   catanh(z) = 0.5 * clog((1 + z) / (1 - z))
 gives wrong results on the positive real axis beginning at 1 (along one of the two branch cuts).
 The formula
   catanh(z) = 0.5 * (clog(1 + z) - clog(1 - z))
 gives correct results.

I've also checked casin, casinf, casinl, casinh, casinhf, and casinhl and the formulae given there
   casin(z) = -i clog(iz + csqrt(1 - z * z))
   casinh(z) = clog(z + csqrt(z * z + 1))
are actually correct.

I suspect that some of these errors are due to somebody trying to "simplify"
these formulae. Since this can ruin the behavior, I recommend to add a comment to the upstream sources that warns against attempt of simplification.

For the sake of reference, I am looking at manpages-dev 3.25-1 on Debian/squeeze.

Quoted 7 lines

Quoted 5 lines

Quoted 5 lines

Bye!
 -richy.
--
Richard B. Kreckel
<http://www.ginac.de/~kreckel/>

This doesn't have seem to have been fixed in the Linux man pages, does it?

Sorry -- no, not yet.

Quoted 9 lines

May I suggest a non-mathematical approach to moving forward?

Just write a little C program and compare the complex results of
a) the library implementation of catan, cacos, etc.,
b) the formula written in the man page, and
c) the formula I proposed above.
If you include a couple of values from all four quadrants, you'll see that my formula always agrees with the implementation while the one from the man page doesn't.

Now there's a good idea ;-). It would have been even better if you'd sent me a program--I don't mess round with complex numbers every day of the week...

Quoted 3 lines

Is there anything else I can do to convince you that this should be fixed? If yes, pretty please, let me know! I would really appreciate getting this fixed.

So, I eventually ironed out my complex math problems and came up with the below:

===
/* Link with "-lm" */

#include <complex.h>
#include <stdlib.h>
#include <unistd.h>
#include <stdio.h>

int
main(int argc, char *argv[])
{
   double complex z, c, f1, f2;
   double complex i = I;

   if (argc != 3) {
       fprintf(stderr, "Usage: %s <real> <imag>\n"); exit(EXIT_FAILURE);
   }

   z = atof(argv[1]) + atof(argv[2]) * I;

   c = cacos(z);

   printf("cacos() = %6.3f %6.3f*i\n", creal(c), cimag(c));

   f1 = -i * clog(z + csqrt(z * z - 1));
   f2 = -i * clog(z + i * csqrt(1 - z * z));

   printf("f1      = %6.3f %6.3f*i\n", creal(f1), cimag(f1)); printf("f2      = %6.3f %6.3f*i\n", creal(f2), cimag(f2));

   exit(EXIT_SUCCESS);
}
===

Some sample runs certainly seem to prove your point for that function at least:

$ ./cacos 1 1
cacos() =  0.905 -1.061*i
f1      =  0.905 -1.061*i
f2      =  0.905 -1.061*i
$ ./cacos -1 -1
cacos() =  2.237  1.061*i
f1      =  2.237  1.061*i
f2      =  2.237  1.061*i
$ ./cacos -1 1
cacos() =  2.237 -1.061*i
f1      = -2.237  1.061*i
f2      =  2.237 -1.061*i
$ ./cacos 1 -1
cacos() =  0.905  1.061*i
f1      = -0.905 -1.061*i
f2      =  0.905  1.061*i

The above is what you meant, right?

Thanks,

Michael

For man-pages-3.33, I have applied the changes below. They should be consistent with what you proposed (and Andries confirmed). Could you please carefully check...

I think it's not so necessary to add a warning to the page source. Nowadays we have version control and logs, so the reasoning should be easy to find.

Hi!

Quoted 1 line · Michael Kerrisk

Quoted 3 lines

For man-pages-3.33, I have applied the changes below. They should be consistent with what you proposed (and Andries confirmed). Could you please carefully check...

Yes, the patch looks correct. Thanks a lot for applying it.

(I would have loved to try it but couldn't find the git repository for man-pages.)

Quoted 3 lines

I think it's not so necessary to add a warning to the page source. Nowadays we have version control and logs, so the reasoning should be easy to find.

Certainly. Out of curiousity: Does the version control of man-pages date back enough to tell when and why the formulas got written the wrong way?

On Thu, Sep 15, 2011 at 11:42 PM, Richard B. Kreckel wrote:

Quoted 2 lines

(I would have loved to try it but couldn't find the git repository for man-pages.)

Well, with the current state of affairs
(http://lwn.net/Articles/457142/), no one can find it...

Hi everybody,

Quoted 3 lines · Michael Kerrisk

Hi Andries,

Since you are the mathematician, can you comment?

Quoted 6 lines · Richard B. Kreckel

Hi everybody,

Quoted 3 lines · Michael Kerrisk

Hi Andries,

Since you are the mathematician, can you comment?

Hi everybody,

Now that I see this discussion again:
Mathematically speaking there is no unique correct definitions. The function is multiple-valued and some arbitrary agreement fixes one of the values.

Richard's point of view is that the actual behavior of the library functions should be documented. Not a bad idea, but perhaps one should document what the standard says, and then file a bug report in case the implementation does something else.

I lost most of my old docs, but just rediscovered a copy of ISO-C-FDIS.1999-04.pdf, I don't know whether that still is the appropriate source to consult.

(i)
It says in 7.3.1: cacos: branch cuts outside [-1,1] along the real axis returns arccos z where Re arccos z lies in [0,pi].
So, if z = cos w then w = arccos z, but because cos (w + 2pi) = cos w and cos (-w) = cos w, we can pick w with real part in [0,pi].

It says in 7.3.7.2: clog: branch cut along the negative real axis returns log z where Im log z lies in [-I*pi,I*pi].

It says in 7.3.8.3: csqrt: branch cut along the negative real axis returns sqrt z where Im sqrt z >= 0.

We want to check whether it is true that
cacos(z) = -I* clog(z + I*csqrt(1 - z * z))

Suppose w = cacos(z). Then w is defined by cos w = z and Re w in [0,pi]. The question is whether iw = log(cos w + i sqrt(sin^2 w)). If sin w has nonnegative imaginary part, then sqrt(sin^2 w) = sin w and the question is whether iw = log(e^{iw}). This logarithm is chosen with imaginary part in [-pi,pi] but w has imaginary part in [0,pi], so this equality holds.
If sin w has negative imaginary part, then sqrt(sin^2 w) = -sin w and the question is whether iw = log(e^{-iw}), and again this holds.

I agree with Richard's equation.

(ii)
Claim: catan(z) = (clog(1 + iz) - clog(1 - iz)) / 2i

7.3.5.3: catan: branch cuts outside [-i,i] along the imaginary axis returns atan z where Re atan z lies in [-pi/2, pi/2].

Suppose w = catan(z). Then w is defined by tan w = z and Re w in [-pi/2, pi/2]. Up to a multiple of 2pi*i the RHS is
(1/2i)log((cos w + i sin w)/(cos w - i sin w)) = (1/2i) log e^{2iw}, looks good. The two multiples of 2pi*i partially cancel, so after dividing by 2i this is in [-pi/2, pi/2].

I agree with Richard's equation.

(iii)
Claim: catanh(z) = (clog(1 + z) - clog(1 - z))/2

7.3.6.3: catanh: branch cuts outside [-1,1] along the real axis returns atanh z where Im atanh z lies in [-pi/2,pi/2].

Suppose w = catanh(z).
Then w is defined by tanh w = z and Re w in [-pi/2,pi/2]. Up to a multiple of 2pi*i the RHS is correct, and just like before we get the right value.

I agree with Richard's equation.

So, if I make no mistake, Richard is correct on all counts.