aboutsummaryrefslogtreecommitdiff
path: root/libf2c/libF77/z_log.c
diff options
context:
space:
mode:
authorToon Moene <toon@moene.indiv.nluug.nl>2000-12-09 16:34:53 +0100
committerToon Moene <toon@gcc.gnu.org>2000-12-09 15:34:53 +0000
commita40bb4d3457a605b5fdb63da04aae6bfa02938e5 (patch)
treee96069addc998cc95cee16e821cf01a3d8ae6b6f /libf2c/libF77/z_log.c
parent976654912051f4b6680973c01da45c3697645f24 (diff)
downloadgcc-a40bb4d3457a605b5fdb63da04aae6bfa02938e5.zip
gcc-a40bb4d3457a605b5fdb63da04aae6bfa02938e5.tar.gz
gcc-a40bb4d3457a605b5fdb63da04aae6bfa02938e5.tar.bz2
Update to Netlib version 20001205.
2000-12-09 Toon Moene <toon@moene.indiv.nluug.nl> Update to Netlib version 20001205. Thanks go to David M. Gay for these updates. * libF77/Version.c: Update version information. * libF77/z_log.c: Improve accuracy of real(log(z)) for z near (+-1,eps) with |eps| small. * libF77/s_cat.c: Adjust call when ftnint and ftnlen are of different size. * libF77/dtime_.c, libF77/etime_.c: Use floating point divide. * libI77/Version.c: Update version information. * libI77/rsne.c, libI77/xwsne.c: Adjust code for when ftnint and ftnlen differ in size. * libI77/lread.c: Fix reading of namelist logical values followed by <name>= where <name> starts with T or F. From-SVN: r38152
Diffstat (limited to 'libf2c/libF77/z_log.c')
-rw-r--r--libf2c/libF77/z_log.c47
1 files changed, 47 insertions, 0 deletions
diff --git a/libf2c/libF77/z_log.c b/libf2c/libF77/z_log.c
index 34c56d4..9dcc7f7 100644
--- a/libf2c/libF77/z_log.c
+++ b/libf2c/libF77/z_log.c
@@ -10,7 +10,54 @@ extern double f__cabs(double, double);
void z_log(doublecomplex *r, doublecomplex *z)
#endif
{
+ double s, s0, t, t2, u, v;
double zi = z->i, zr = z->r;
+
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*/
+
+ 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);
+ } while(s > s0);
+ r->r = s*t;
+ }
+#endif
}