z_log.c

Go to the documentation of this file.

#include "f2c.h"

#ifdef KR_headers
double log(), f__cabs(), atan2();
#define ANSI(x) ()
#else
#define ANSI(x) x
#undef abs
#include "math.h"
#ifdef __cplusplus
extern "C" {
#endif
extern double f__cabs(double, double);
#endif

#ifndef NO_DOUBLE_EXTENDED
#ifndef GCC_COMPARE_BUG_FIXED
#ifndef Pre20000310
#ifdef Comment
Some versions of gcc, such as 2.95.3 and 3.0.4, are buggy under -O2 or -O3:
on IA32 (Intel 80x87) systems, they may do comparisons on values computed
in extended-precision registers.  This can lead to the test "s > s0" that
was used below being carried out incorrectly.  The fix below cannot be
spoiled by overzealous optimization, since the compiler cannot know
whether gcc_bug_bypass_diff_F2C will be nonzero.  (We expect it always
to be zero.  The weird name is unlikely to collide with anything.)

An example (provided by Ulrich Jakobus) where the bug fix matters is

        double complex a, b
        a = (.1099557428756427618354862829619, .9857360542953131909982289471372)
        b = log(a)

An alternative to the fix below would be to use 53-bit rounding precision,
but the means of specifying this 80x87 feature are highly unportable.
#endif /*Comment*/
#define BYPASS_GCC_COMPARE_BUG
double (*gcc_bug_bypass_diff_F2C) ANSI((double*,double*));
 static double
#ifdef KR_headers
diff1(a,b) double *a, *b;
#else
diff1(double *a, double *b)
#endif
{ return *a - *b; }
#endif /*Pre20000310*/
#endif /*GCC_COMPARE_BUG_FIXED*/
#endif /*NO_DOUBLE_EXTENDED*/

#ifdef KR_headers
VOID z_log(r, z) doublecomplex *r, *z;
#else
void z_log(doublecomplex *r, doublecomplex *z)
#endif
{
        double s, s0, t, t2, u, v;
        double zi = z->i, zr = z->r;
#ifdef BYPASS_GCC_COMPARE_BUG
        double (*diff) ANSI((double*,double*));
#endif

        r->i = atan2(zi, zr);
#ifdef Pre20000310
        r->r = log( f__cabs( zr, zi ) );
#else
        if (zi < 0)
                zi = -zi;
        if (zr < 0)
                zr = -zr;
        if (zr < zi) {
                t = zi;
                zi = zr;
                zr = t;
                }
        t = zi/zr;
        s = zr * sqrt(1 + t*t);
        /* now s = f__cabs(zi,zr), and zr = |zr| >= |zi| = zi */
        if ((t = s - 1) < 0)
                t = -t;
        if (t > .01)
                r->r = log(s);
        else {

#ifdef Comment

        log(1+x) = x - x^2/2 + x^3/3 - x^4/4 + - ...

                 = x(1 - x/2 + x^2/3 -+...)

        [sqrt(y^2 + z^2) - 1] * [sqrt(y^2 + z^2) + 1] = y^2 + z^2 - 1, so

        sqrt(y^2 + z^2) - 1 = (y^2 + z^2 - 1) / [sqrt(y^2 + z^2) + 1]

#endif /*Comment*/

#ifdef BYPASS_GCC_COMPARE_BUG
                if (!(diff = gcc_bug_bypass_diff_F2C))
                        diff = diff1;
#endif
                t = ((zr*zr - 1.) + zi*zi) / (s + 1);
                t2 = t*t;
                s = 1. - 0.5*t;
                u = v = 1;
                do {
                        s0 = s;
                        u *= t2;
                        v += 2;
                        s += u/v - t*u/(v+1);
                        }
#ifdef BYPASS_GCC_COMPARE_BUG
                        while(s - s0 > 1e-18 || (*diff)(&s,&s0) > 0.);
#else
                        while(s > s0);
#endif
                r->r = s*t;
                }
#endif
        }
#ifdef __cplusplus
}
#endif