Commit 9cfd948e by Craig Burley Committed by Craig Burley

revert back to netlib versions as of f2c-19990501

From-SVN: r26739
parent d0369a2f
Mon May 3 10:52:53 1999 Craig Burley <craig@jcb-sc.com>
* libF77/c_cos.c, libF77/c_div.c, libF77/c_exp.c, libF77/c_log.c,
libF77/c_sin.c, libF77/c_sqrt.c, libF77/d_cnjg.c, libF77/pow_zi.c,
libF77/r_cnjg.c, libF77/z_cos.c, libF77/z_div.c, libF77/z_exp.c,
libF77/z_log.c, libF77/z_sin.c, libF77/z_sqrt.c: Revert back to
netlib versions as of f2c-19990501.
Sun May 2 01:38:50 1999 Craig Burley <craig@jcb-sc.com> Sun May 2 01:38:50 1999 Craig Burley <craig@jcb-sc.com>
* libU77/u77-test.f (main): Declare FTELL as intrinsic. * libU77/u77-test.f (main): Declare FTELL as intrinsic.
......
...@@ -3,19 +3,15 @@ ...@@ -3,19 +3,15 @@
#ifdef KR_headers #ifdef KR_headers
extern double sin(), cos(), sinh(), cosh(); extern double sin(), cos(), sinh(), cosh();
VOID c_cos(resx, z) complex *resx, *z; VOID c_cos(r, z) complex *r, *z;
#else #else
#undef abs #undef abs
#include <math.h> #include "math.h"
void c_cos(complex *resx, complex *z) void c_cos(complex *r, complex *z)
#endif #endif
{ {
complex res; double zr = z->r;
r->r = cos(zr) * cosh(z->i);
res.r = cos(z->r) * cosh(z->i); r->i = - sin(zr) * sinh(z->i);
res.i = - sin(z->r) * sinh(z->i); }
resx->r = res.r;
resx->i = res.i;
}
...@@ -2,39 +2,36 @@ ...@@ -2,39 +2,36 @@
#ifdef KR_headers #ifdef KR_headers
extern VOID sig_die(); extern VOID sig_die();
VOID c_div(resx, a, b) VOID c_div(c, a, b)
complex *a, *b, *resx; complex *a, *b, *c;
#else #else
extern void sig_die(char*,int); extern void sig_die(char*,int);
void c_div(complex *resx, complex *a, complex *b) void c_div(complex *c, complex *a, complex *b)
#endif #endif
{ {
double ratio, den; double ratio, den;
double abr, abi; double abr, abi, cr;
complex res;
if( (abr = b->r) < 0.) if( (abr = b->r) < 0.)
abr = - abr; abr = - abr;
if( (abi = b->i) < 0.) if( (abi = b->i) < 0.)
abi = - abi; abi = - abi;
if( abr <= abi ) if( abr <= abi )
{ {
if(abi == 0) if(abi == 0)
sig_die("complex division by zero", 1); sig_die("complex division by zero", 1);
ratio = (double)b->r / b->i ; ratio = (double)b->r / b->i ;
den = b->i * (1 + ratio*ratio); den = b->i * (1 + ratio*ratio);
res.r = (a->r*ratio + a->i) / den; cr = (a->r*ratio + a->i) / den;
res.i = (a->i*ratio - a->r) / den; c->i = (a->i*ratio - a->r) / den;
} }
else else
{ {
ratio = (double)b->i / b->r ; ratio = (double)b->i / b->r ;
den = b->r * (1 + ratio*ratio); den = b->r * (1 + ratio*ratio);
res.r = (a->r + a->i*ratio) / den; cr = (a->r + a->i*ratio) / den;
res.i = (a->i - a->r*ratio) / den; c->i = (a->i - a->r*ratio) / den;
}
c->r = cr;
} }
resx->r = res.r;
resx->i = res.i;
}
...@@ -3,21 +3,17 @@ ...@@ -3,21 +3,17 @@
#ifdef KR_headers #ifdef KR_headers
extern double exp(), cos(), sin(); extern double exp(), cos(), sin();
VOID c_exp(resx, z) complex *resx, *z; VOID c_exp(r, z) complex *r, *z;
#else #else
#undef abs #undef abs
#include <math.h> #include "math.h"
void c_exp(complex *resx, complex *z) void c_exp(complex *r, complex *z)
#endif #endif
{ {
double expx; double expx;
complex res;
expx = exp(z->r); expx = exp(z->r);
res.r = expx * cos(z->i); r->r = expx * cos(z->i);
res.i = expx * sin(z->i); r->i = expx * sin(z->i);
resx->r = res.r;
resx->i = res.i;
} }
...@@ -2,20 +2,16 @@ ...@@ -2,20 +2,16 @@
#ifdef KR_headers #ifdef KR_headers
extern double log(), f__cabs(), atan2(); extern double log(), f__cabs(), atan2();
VOID c_log(resx, z) complex *resx, *z; VOID c_log(r, z) complex *r, *z;
#else #else
#undef abs #undef abs
#include <math.h> #include "math.h"
extern double f__cabs(double, double); extern double f__cabs(double, double);
void c_log(complex *resx, complex *z) void c_log(complex *r, complex *z)
#endif #endif
{ {
complex res; double zi;
r->i = atan2(zi = z->i, z->r);
res.i = atan2(z->i, z->r); r->r = log( f__cabs(z->r, zi) );
res.r = log( f__cabs(z->r, z->i) ); }
resx->r = res.r;
resx->i = res.i;
}
...@@ -3,19 +3,15 @@ ...@@ -3,19 +3,15 @@
#ifdef KR_headers #ifdef KR_headers
extern double sin(), cos(), sinh(), cosh(); extern double sin(), cos(), sinh(), cosh();
VOID c_sin(resx, z) complex *resx, *z; VOID c_sin(r, z) complex *r, *z;
#else #else
#undef abs #undef abs
#include <math.h> #include "math.h"
void c_sin(complex *resx, complex *z) void c_sin(complex *r, complex *z)
#endif #endif
{ {
complex res; double zr = z->r;
r->r = sin(zr) * cosh(z->i);
res.r = sin(z->r) * cosh(z->i); r->i = cos(zr) * sinh(z->i);
res.i = cos(z->r) * sinh(z->i); }
resx->r = res.r;
resx->i = res.i;
}
...@@ -3,36 +3,33 @@ ...@@ -3,36 +3,33 @@
#ifdef KR_headers #ifdef KR_headers
extern double sqrt(), f__cabs(); extern double sqrt(), f__cabs();
VOID c_sqrt(resx, z) complex *resx, *z; VOID c_sqrt(r, z) complex *r, *z;
#else #else
#undef abs #undef abs
#include <math.h> #include "math.h"
extern double f__cabs(double, double); extern double f__cabs(double, double);
void c_sqrt(complex *resx, complex *z) void c_sqrt(complex *r, complex *z)
#endif #endif
{ {
double mag, t; double mag, t;
complex res; double zi = z->i, zr = z->r;
if( (mag = f__cabs(z->r, z->i)) == 0.) if( (mag = f__cabs(zr, zi)) == 0.)
res.r = res.i = 0.; r->r = r->i = 0.;
else if(z->r > 0) else if(zr > 0)
{ {
res.r = t = sqrt(0.5 * (mag + z->r) ); r->r = t = sqrt(0.5 * (mag + zr) );
t = z->i / t; t = zi / t;
res.i = 0.5 * t; r->i = 0.5 * t;
}
else
{
t = sqrt(0.5 * (mag - zr) );
if(zi < 0)
t = -t;
r->i = t;
t = zi / t;
r->r = 0.5 * t;
}
} }
else
{
t = sqrt(0.5 * (mag - z->r) );
if(z->i < 0)
t = -t;
res.i = t;
t = z->i / t;
res.r = 0.5 * t;
}
resx->r = res.r;
resx->i = res.i;
}
...@@ -2,16 +2,11 @@ ...@@ -2,16 +2,11 @@
VOID VOID
#ifdef KR_headers #ifdef KR_headers
d_cnjg(resx, z) doublecomplex *resx, *z; d_cnjg(r, z) doublecomplex *r, *z;
#else #else
d_cnjg(doublecomplex *resx, doublecomplex *z) d_cnjg(doublecomplex *r, doublecomplex *z)
#endif #endif
{ {
doublecomplex res; r->r = z->r;
r->i = - z->i;
res.r = z->r;
res.i = - z->i;
resx->r = res.r;
resx->i = res.i;
} }
#include "f2c.h" #include "f2c.h"
#ifdef KR_headers #ifdef KR_headers
VOID pow_zi(resx, a, b) /* p = a**b */ VOID pow_zi(p, a, b) /* p = a**b */
doublecomplex *resx, *a; integer *b; doublecomplex *p, *a; integer *b;
#else #else
extern void z_div(doublecomplex*, doublecomplex*, doublecomplex*); extern void z_div(doublecomplex*, doublecomplex*, doublecomplex*);
void pow_zi(doublecomplex *resx, doublecomplex *a, integer *b) /* p = a**b */ void pow_zi(doublecomplex *p, doublecomplex *a, integer *b) /* p = a**b */
#endif #endif
{ {
integer n; integer n;
unsigned long u; unsigned long u;
double t; double t;
doublecomplex x; doublecomplex q, x;
doublecomplex res; static doublecomplex one = {1.0, 0.0};
static doublecomplex one = {1.0, 0.0};
n = *b; n = *b;
q.r = 1;
q.i = 0;
if(n == 0) if(n == 0)
{ goto done;
resx->r = 1; if(n < 0)
resx->i = 0;
return;
}
res.r = 1;
res.i = 0;
if(n < 0)
{
n = -n;
z_div(&x, &one, a);
}
else
{
x.r = a->r;
x.i = a->i;
}
for(u = n; ; )
{
if(u & 01)
{ {
t = res.r * x.r - res.i * x.i; n = -n;
res.i = res.r * x.i + res.i * x.r; z_div(&x, &one, a);
res.r = t;
} }
if(u >>= 1) else
{ {
t = x.r * x.r - x.i * x.i; x.r = a->r;
x.i = 2 * x.r * x.i; x.i = a->i;
x.r = t;
} }
else
break;
}
resx->r = res.r; for(u = n; ; )
resx->i = res.i; {
} if(u & 01)
{
t = q.r * x.r - q.i * x.i;
q.i = q.r * x.i + q.i * x.r;
q.r = t;
}
if(u >>= 1)
{
t = x.r * x.r - x.i * x.i;
x.i = 2 * x.r * x.i;
x.r = t;
}
else
break;
}
done:
p->i = q.i;
p->r = q.r;
}
#include "f2c.h" #include "f2c.h"
#ifdef KR_headers #ifdef KR_headers
VOID r_cnjg(resx, z) complex *resx, *z; VOID r_cnjg(r, z) complex *r, *z;
#else #else
VOID r_cnjg(complex *resx, complex *z) VOID r_cnjg(complex *r, complex *z)
#endif #endif
{ {
complex res; r->r = z->r;
r->i = - z->i;
res.r = z->r;
res.i = - z->i;
resx->r = res.r;
resx->i = res.i;
} }
...@@ -2,18 +2,14 @@ ...@@ -2,18 +2,14 @@
#ifdef KR_headers #ifdef KR_headers
double sin(), cos(), sinh(), cosh(); double sin(), cos(), sinh(), cosh();
VOID z_cos(resx, z) doublecomplex *resx, *z; VOID z_cos(r, z) doublecomplex *r, *z;
#else #else
#undef abs #undef abs
#include <math.h> #include "math.h"
void z_cos(doublecomplex *resx, doublecomplex *z) void z_cos(doublecomplex *r, doublecomplex *z)
#endif #endif
{ {
doublecomplex res; double zr = z->r;
r->r = cos(zr) * cosh(z->i);
res.r = cos(z->r) * cosh(z->i); r->i = - sin(zr) * sinh(z->i);
res.i = - sin(z->r) * sinh(z->i); }
resx->r = res.r;
resx->i = res.i;
}
...@@ -2,38 +2,35 @@ ...@@ -2,38 +2,35 @@
#ifdef KR_headers #ifdef KR_headers
extern VOID sig_die(); extern VOID sig_die();
VOID z_div(resx, a, b) doublecomplex *a, *b, *resx; VOID z_div(c, a, b) doublecomplex *a, *b, *c;
#else #else
extern void sig_die(char*, int); extern void sig_die(char*, int);
void z_div(doublecomplex *resx, doublecomplex *a, doublecomplex *b) void z_div(doublecomplex *c, doublecomplex *a, doublecomplex *b)
#endif #endif
{ {
double ratio, den; double ratio, den;
double abr, abi; double abr, abi, cr;
doublecomplex res;
if( (abr = b->r) < 0.) if( (abr = b->r) < 0.)
abr = - abr; abr = - abr;
if( (abi = b->i) < 0.) if( (abi = b->i) < 0.)
abi = - abi; abi = - abi;
if( abr <= abi ) if( abr <= abi )
{ {
if(abi == 0) if(abi == 0)
sig_die("complex division by zero", 1); sig_die("complex division by zero", 1);
ratio = b->r / b->i ; ratio = b->r / b->i ;
den = b->i * (1 + ratio*ratio); den = b->i * (1 + ratio*ratio);
res.r = (a->r*ratio + a->i) / den; cr = (a->r*ratio + a->i) / den;
res.i = (a->i*ratio - a->r) / den; c->i = (a->i*ratio - a->r) / den;
} }
else else
{ {
ratio = b->i / b->r ; ratio = b->i / b->r ;
den = b->r * (1 + ratio*ratio); den = b->r * (1 + ratio*ratio);
res.r = (a->r + a->i*ratio) / den; cr = (a->r + a->i*ratio) / den;
res.i = (a->i - a->r*ratio) / den; c->i = (a->i - a->r*ratio) / den;
}
c->r = cr;
} }
resx->r = res.r;
resx->i = res.i;
}
...@@ -2,20 +2,16 @@ ...@@ -2,20 +2,16 @@
#ifdef KR_headers #ifdef KR_headers
double exp(), cos(), sin(); double exp(), cos(), sin();
VOID z_exp(resx, z) doublecomplex *resx, *z; VOID z_exp(r, z) doublecomplex *r, *z;
#else #else
#undef abs #undef abs
#include <math.h> #include "math.h"
void z_exp(doublecomplex *resx, doublecomplex *z) void z_exp(doublecomplex *r, doublecomplex *z)
#endif #endif
{ {
double expx; double expx;
doublecomplex res;
expx = exp(z->r); expx = exp(z->r);
res.r = expx * cos(z->i); r->r = expx * cos(z->i);
res.i = expx * sin(z->i); r->i = expx * sin(z->i);
resx->r = res.r;
resx->i = res.i;
} }
...@@ -2,19 +2,15 @@ ...@@ -2,19 +2,15 @@
#ifdef KR_headers #ifdef KR_headers
double log(), f__cabs(), atan2(); double log(), f__cabs(), atan2();
VOID z_log(resx, z) doublecomplex *resx, *z; VOID z_log(r, z) doublecomplex *r, *z;
#else #else
#undef abs #undef abs
#include <math.h> #include "math.h"
extern double f__cabs(double, double); extern double f__cabs(double, double);
void z_log(doublecomplex *resx, doublecomplex *z) void z_log(doublecomplex *r, doublecomplex *z)
#endif #endif
{ {
doublecomplex res; double zi = z->i;
r->i = atan2(zi, z->r);
res.i = atan2(z->i, z->r); r->r = log( f__cabs( z->r, zi ) );
res.r = log( f__cabs( z->r, z->i ) ); }
resx->r = res.r;
resx->i = res.i;
}
...@@ -2,18 +2,14 @@ ...@@ -2,18 +2,14 @@
#ifdef KR_headers #ifdef KR_headers
double sin(), cos(), sinh(), cosh(); double sin(), cos(), sinh(), cosh();
VOID z_sin(resx, z) doublecomplex *resx, *z; VOID z_sin(r, z) doublecomplex *r, *z;
#else #else
#undef abs #undef abs
#include <math.h> #include "math.h"
void z_sin(doublecomplex *resx, doublecomplex *z) void z_sin(doublecomplex *r, doublecomplex *z)
#endif #endif
{ {
doublecomplex res; double zr = z->r;
r->r = sin(zr) * cosh(z->i);
res.r = sin(z->r) * cosh(z->i); r->i = cos(zr) * sinh(z->i);
res.i = cos(z->r) * sinh(z->i); }
resx->r = res.r;
resx->i = res.i;
}
...@@ -2,32 +2,28 @@ ...@@ -2,32 +2,28 @@
#ifdef KR_headers #ifdef KR_headers
double sqrt(), f__cabs(); double sqrt(), f__cabs();
VOID z_sqrt(resx, z) doublecomplex *resx, *z; VOID z_sqrt(r, z) doublecomplex *r, *z;
#else #else
#undef abs #undef abs
#include <math.h> #include "math.h"
extern double f__cabs(double, double); extern double f__cabs(double, double);
void z_sqrt(doublecomplex *resx, doublecomplex *z) void z_sqrt(doublecomplex *r, doublecomplex *z)
#endif #endif
{ {
double mag; double mag, zi = z->i, zr = z->r;
doublecomplex res;
if( (mag = f__cabs(z->r, z->i)) == 0.) if( (mag = f__cabs(zr, zi)) == 0.)
res.r = res.i = 0.; r->r = r->i = 0.;
else if(z->r > 0) else if(zr > 0)
{ {
res.r = sqrt(0.5 * (mag + z->r) ); r->r = sqrt(0.5 * (mag + zr) );
res.i = z->i / res.r / 2; r->i = zi / r->r / 2;
}
else
{
r->i = sqrt(0.5 * (mag - zr) );
if(zi < 0)
r->i = - r->i;
r->r = zi / r->i / 2;
}
} }
else
{
res.i = sqrt(0.5 * (mag - z->r) );
if(z->i < 0)
res.i = - res.i;
res.r = z->i / res.i / 2;
}
resx->r = res.r;
resx->i = res.i;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment