Use of C99 extra long double math functions after r236148
Stephen Montgomery-Smith
stephen at missouri.edu
Sun Aug 12 23:02:44 UTC 2012
I just realized that catan(z) = reverse( catanh(reverse (z))), just like
casin relates to casinh (remember reverse(x+I*y) = y+I*x). This is a
consequence of catan and catanh being odd functions, as well as the
standard relation catan(z) = -I*catanh(I*z).
So I would modify Peter's code by taking out the minus signs.
Maybe it would make a difference if the answers involved -0.
On 07/22/12 07:12, Peter Jeremy wrote:
> /*
> * Arc-tangent of a complex argument z = x + I*y.
> *
> * catan(z) = reverse( catanh(reverse (z)))
> *
> */
> double complex
> catan(double complex z)
> {
> double complex r;
>
> r = catanh(cpack(cimag(z), creal(z)));
> return (cpack(cimag(r), creal(r)));
> }
>
I am attaching code that computes all six arc-trig-hyp functions. Peter
made a remark that catanh(z) would be hard to compute when |z|=1.
Fortunately that is not the case, because when |z|=1, the imaginary part
of catanh(z) is plus or minus PI/4. So we won't face the same problems
that clog(z) has.
I feel that I am done with these functions for now. I tried to change
my comments to conform to the style given to me by Bruce. However
spacing inside mathematical expressions is something where I am
inconsistent.
The functions still need a lot of work to handle -0, infs and NaNs
correctly. I will leave that to you guys, because you seem so much
better at it than me. I still don't understand why the proper test is
"if (x!=x) return(x+x)" rather than "if (isnan(x)) return(NAN)".
However, I think you might find that a lot of the handling when the
input or output is infinity works without any changes. For example,
catanh(1) seems to produce the correct answer, and maybe even
catanh(1-0*I) will work better than expected.
I am also attaching the test code. I run it like this, so that only
ULPs greater than 3 appear:
./test3 | perl -lne '@a=split " ",$_;print if $a[2]>3 || $a[3]>3'
and get outputs like this:
21987 atanh 3.13643 0.299452 0.443282 0.0665108 0.473296 0.0824781
67013 atan 0.377411 3.03922 -0.0170315 -0.442191 -0.0211662 -0.474753
70044 acosh 0.883474 3.23108 1.06353 0.107343 0.433696 0.242278
70044 acos 3.23108 0.883474 1.06353 0.107343 0.242278 -0.433696
96124 atan 0.509279 3.21631 -0.0346851 -0.461121 -0.0440054 -0.497841
The first example is the count of which example I am trying. The third
and fourth entries are the ULPs of the real and imaginary parts of the
answer. The 5th and 6th entries are the real and imaginary parts of the
input, and the 7th and 8th entries are the real and imaginary parts of
the answer.
I use the unuran port to generate data where the x and y values of the
inputs are normally distributed N(0,1). As you can see from the parts I
commented out, I tried many variations (mostly to check edge cases:
close to zero, very large, or close to branch cuts). In particular, if
you want random data on the unit disk, use h=hypot(x,y); x/=h; y/=h;
My next project will be to get clog(z) to work well when |z|=1.
Stephen
-------------- next part --------------
#include <complex.h>
#include <float.h>
#include <math.h>
#include "math_private.h"
/*
* gcc doesn't implement complex multiplication or division correctly,
* so we need to handle infinities specially. We turn on this pragma to
* notify conforming c99 compilers that the fast-but-incorrect code that
* gcc generates is acceptable, since the special cases have already been
* handled.
*/
#pragma STDC CX_LIMITED_RANGE ON
complex double clog(complex double z);
static const double
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
huge= 1.00000000000000000000e+300;
/*
* Testing indicates that all these functions are accurate up to 4 ULP.
*/
/*
* The algorithm is very close to that in "Implementing the complex arcsine
* and arccosine functions using exception handling" by T. E. Hull,
* Thomas F. Fairgrieve, and Ping Tak Peter Tang, published in ACM
* Transactions on Mathematical Software, Volume 23 Issue 3, 1997, Pages
* 299-335, http://dl.acm.org/citation.cfm?id=275324
*
* casinh(x+iy) = sign(x)*log(A+sqrt(A*A-1)) + sign(y)*I*asin(B)
* where
* A = 0.5(|z+I| + |z-I|) = f(x,1+y) + f(x,1-y) + 1
* B = 0.5(|z+I| - |z-I|)
* z = x+I*y
* f(x,y) = 0.5*(hypot(x,y)-y)
* We also use
* asin(B) = atan2(sqrt(A*A-y*y),y)
* A-y = f(x,y+1) + f(x,y-1).
*
* Much of the difficulty comes because computing f(x,y) may produce
* underflows.
*/
/*
* Returns 0.5*(hypot(x,y)-y). It assumes x is positive, and that y does
* not satisfy the inequalities 0 < fabs(y) < 1e-20.
* If reporting the answer risks an underflow, the underflow flag is set,
*and it returns 0.5*(hypot(x,y)-y)/x/x.
*/
static double f(double x, double y, int *underflow) {
if (x==0) {
*underflow = 0;
if (y > 0)
return 0;
return -y;
}
if (y==0) {
*underflow = 0;
return 0.5*x;
}
if (x < 1e-100 && x < y) {
*underflow = 1;
return 0.5/(hypot(x,y)+y);
}
if (x < y) {
*underflow = 0;
return 0.5*x*x/(hypot(x,y)+y);
}
*underflow = 0;
return 0.5*(hypot(x,y)-y);
}
/*
* All the hard work is contained in this function.
* Upon return:
* rx = Re(casinh(x+I*y))
* B_good is set to 1 if the value of B is usable.
* If B_good is set to 0, A2my2 = A*A-y*y.
*/
static void do_hard_work(double x, double y, double *rx, int *B_good, double *B, double *A2my2)
{
double R, S, A, fp, fm;
int fpuf, fmuf;
R = hypot(x,y+1);
S = hypot(x,y-1);
A = 0.5*(R + S);
if (A < 10) {
fp = f(x,1+y,&fpuf);
fm = f(x,1-y,&fmuf);
if (fpuf == 1 && fmuf == 1) {
if (huge+x>one) /* set inexact flag. */
*rx = log1p(x*sqrt((fp+fm)*(A+1)));
} else if (fmuf == 1) {
/* Overflow not possible because fp < 1e50 and x > 1e-100.
Underflow not possible because either fm=0 or fm
approximately bigger than 1e-200. */
if (huge+x>one) /* set inexact flag. */
*rx = log1p(fp+sqrt(x)*sqrt((fp/x+fm*x)*(A+1)));
} else if (fpuf == 1) {
/* Similar arguments against over/underflow. */
if (huge+x>one) /* set inexact flag. */
*rx = log1p(fm+sqrt(x)*sqrt((fm/x+fp*x)*(A+1)));
} else {
*rx = log1p(fp + fm + sqrt((fp+fm)*(A+1)));
}
} else
*rx = log(A + sqrt(A*A-1));
*B = y/A; /* = 0.5*(R - S) */
*B_good = 1;
if (*B > 0.5) {
*B_good = 0;
fp = f(x,y+1,&fpuf);
fm = f(x,y-1,&fmuf);
if (fpuf == 1 && fmuf == 1)
*A2my2 =x*sqrt((A+y)*(fp+fm));
else if (fmuf == 1)
/* Overflow not possible because fp < 1e50 and x > 1e-100.
Underflow not possible because either fm=0 or fm
approximately bigger than 1e-200. */
*A2my2 = sqrt(x)*sqrt((A+y)*(fp/x+fm*x));
else if (fpuf == 1)
/* Similar arguments against over/underflow. */
*A2my2 = sqrt(x)*sqrt((A+y)*(fm/x+fp*x));
else
*A2my2 = sqrt((A+y)*(fp+fm));
}
}
double complex
casinh(double complex z)
{
double x, y, rx, ry, B, A2my2;
int sx, sy;
int B_good;
x = creal(z);
y = cimag(z);
sx = signbit(x);
sy = signbit(y);
x = fabs(x);
y = fabs(y);
if (cabs(z) > 1e20) {
if (huge+x>one) { /* set inexact flag. */
if (sx == 0) return clog(2*z);
if (sx == 1) return -clog(-2*z);
}
}
if (cabs(z) < 1e-20)
if (huge+x>one) /* set inexact flag. */
return z;
do_hard_work(x, y, &rx, &B_good, &B, &A2my2);
if (B_good)
ry = asin(B);
else
ry = atan2(y,A2my2);
if (sx == 1) rx = -rx;
if (sy == 1) ry = -ry;
return cpack(rx,ry);
}
/*
* casin(z) = reverse(casinh(reverse(z)))
* where reverse(x+I*y) = y+x*I = I*conj(x+I*y).
*/
double complex
casin(double complex z)
{
complex result;
result = casinh(cpack(cimag(z),creal(z)));
return cpack(cimag(result),creal(result));
}
/*
* cacos(z) = PI/2 - casin(z)
* but do the computation carefully so cacos(z) is accurate when z is close to 1.
*/
double complex
cacos(double complex z)
{
double x, y, rx, ry, B, A2my2;
int sx, sy;
int B_good;
complex w;
x = creal(z);
y = cimag(z);
sx = signbit(x);
sy = signbit(y);
x = fabs(x);
y = fabs(y);
if (cabs(z) > 1e20) {
if (huge+x>one) { /* set inexact flag. */
w = clog(2*z);
if (signbit(cimag(w)) == 0)
return cpack(cimag(w),-creal(w));
return cpack(-cimag(w),creal(w));
}
}
if (cabs(z) < 1e-10)
if (huge+x>one) /* set inexact flag. */
return cpack(M_PI_2-creal(z),-cimag(z));
do_hard_work(y, x, &ry, &B_good, &B, &A2my2);
if (B_good) {
if (sx==0)
rx = acos(B);
else
rx = acos(-B);
} else {
if (sx==0)
rx = atan2(A2my2,x);
else
rx = atan2(A2my2,-x);
}
if (sy==0) ry = -ry;
return cpack(rx,ry);
}
/*
* cacosh(z) = I*cacos(z) or -I*cacos(z)
* where the sign is chosen so Re(cacosh(z)) >= 0 .
*/
double complex
cacosh(double complex z)
{
complex double w;
w = cacos(z);
if (signbit(cimag(w)) == 0)
return cpack(cimag(w),-creal(w));
else
return cpack(-cimag(w),creal(w));
}
/*
* catanh(z) = 0.25 * log((z+1)/(z-1))
* = 0.25 * log(|z+1|/|z-1|) + 0.5 * I * atan2(2y/(1-x*x-y*y))
*/
double complex
catanh(double complex z)
{
double x, y, rx, ry, hp, hm;
x = creal(z);
y = cimag(z);
if (cabs(z) < 1e-20)
if (huge+x>one) /* set inexact flag. */
return z;
if (cabs(z) > 1e20)
if (huge+x>one) { /* set inexact flag. */
if (signbit(x) == 0)
return cpack(0,M_PI_2);
return cpack(0,-M_PI_2);
}
if (fabs(y) < 1e-100) {
hp = (x+1)*(x+1);
hm = (x-1)*(x-1);
} else {
hp = (x+1)*(x+1)+y*y; /* |z+1| */
hm = (x-1)*(x-1)+y*y; /* |z-1| */
}
if (hp < 0.5 || hm < 0.5)
rx = 0.25*(log(hp/hm));
else if (x > 0)
rx = 0.25*log1p(4*x/hm);
else
rx = -0.25*log1p(-4*x/hp);
if (x==1 || x==-1) {
if (signbit(y) == 0)
ry = atan2(2, -y)/2;
else
ry = atan2(-2, y)/2;
} else if (fabs(y) < 1e-100) {
if (huge+x>one) /* set inexact flag. */
ry = atan2(2*y, (1-x)*(1+x))/2;
} else
ry = atan2(2*y, (1-x)*(1+x)-y*y)/2;
return cpack(rx,ry);
}
/*
* catan(z) = reverse(catanh(reverse(z)))
* where reverse(x+I*y) = y+x*I = I*conj(x+I*y).
*/
double complex
catan(double complex z)
{
complex result;
result = catanh(cpack(cimag(z),creal(z)));
return cpack(cimag(result),creal(result));
}
-------------- next part --------------
#include <stdio.h>
#include <stdlib.h>
#include <sys/cdefs.h>
#include <float.h>
#include <fenv.h>
#include "complex.h"
#include "math.h"
#include "math_private.h"
#include "mpfr.h"
#include "mpc.h"
#include <unuran.h>
complex double casinh(complex double z);
complex double casin(complex double z);
complex double cacosh(complex double z);
complex double cacos(complex double z);
complex double catanh(complex double z);
complex double catan(complex double z);
mpc_t zz, rr;
mpfr_t rxx, ryy, exx, eyy;
void eval(double x,double y,complex double (*f)(complex double),void (*mpc_f)(mpc_t,mpc_t,mpc_rnd_t),double *rx,double *ry,double *ex,double *ey) {
complex double result;
result = f(cpack(x,y));
*rx = creal(result);
*ry = cimag(result);
mpc_set_d_d(zz, x, y, MPC_RNDNN);
mpc_f(rr, zz, MPC_RNDNN);
/*
mpc_out_str(stdout, 10, 100, zz, MPC_RNDNN);
puts("");
mpc_out_str(stdout, 10, 100, rr, MPC_RNDNN);
puts("");
*/
mpc_real(rxx, rr, MPFR_RNDN);
mpfr_sub_d(exx,rxx,*rx,MPFR_RNDN);
mpfr_abs(exx,exx,MPFR_RNDN);
mpfr_mul_2si(exx,exx,- mpfr_get_exp(rxx)+DBL_MANT_DIG,MPFR_RNDN);
*ex = mpfr_get_d(exx,MPFR_RNDN);
mpc_imag(ryy, rr, MPFR_RNDN);
mpfr_sub_d(eyy,ryy,*ry,MPFR_RNDN);
mpfr_abs(eyy,eyy,MPFR_RNDN);
mpfr_mul_2si(eyy,eyy,- mpfr_get_exp(ryy)+DBL_MANT_DIG,MPFR_RNDN);
*ey = mpfr_get_d(eyy,MPFR_RNDN);
}
int main(int argc, char **argv) {
double x,y,rx,ry,ex,ey,h;
UNUR_DISTR *distr;
UNUR_PAR *par;
UNUR_GEN *gen;
int count = 0;
distr = unur_distr_normal(NULL,0);
par = unur_auto_new(distr);
gen = unur_init(par);
mpc_init2(zz,300);
mpc_init2(rr,300);
mpfr_set_default_prec(300);
mpfr_init(rxx);
mpfr_init(ryy);
mpfr_init(exx);
mpfr_init(eyy);
while (1) {
x = unur_sample_cont(gen);
y = unur_sample_cont(gen);
/*
h = hypot(x,y);
x = x/h;
y = y/h;
x = 0;
y = 1;
x += 1e-200*unur_sample_cont(gen);
y += 1e-200*unur_sample_cont(gen);
*/
count++;
eval(x,y,casinh,mpc_asinh,&rx,&ry,&ex,&ey);
printf("%d asinh %g %g %g %g %g %g\n",count,ex,ey,x,y,rx,ry);
eval(x,y,cacosh,mpc_acosh,&rx,&ry,&ex,&ey);
printf("%d acosh %g %g %g %g %g %g\n",count,ex,ey,x,y,rx,ry);
eval(x,y,catanh,mpc_atanh,&rx,&ry,&ex,&ey);
printf("%d atanh %g %g %g %g %g %g\n",count,ex,ey,x,y,rx,ry);
eval(x,y,casin,mpc_asin,&rx,&ry,&ex,&ey);
printf("%d asin %g %g %g %g %g %g\n",count,ex,ey,x,y,rx,ry);
eval(x,y,cacos,mpc_acos,&rx,&ry,&ex,&ey);
printf("%d acos %g %g %g %g %g %g\n",count,ex,ey,x,y,rx,ry);
eval(x,y,catan,mpc_atan,&rx,&ry,&ex,&ey);
printf("%d atan %g %g %g %g %g %g\n",count,ex,ey,x,y,rx,ry);
}
}
More information about the freebsd-numerics
mailing list