fish/icalc/cmath.c

   1 /*
   2 *       Math routines for complex-number expression parser.
   3 *       In the routines, rv is variable holding 'return value'.
   4 *
   5 *       MWS, March 17, 1991.
   6 */
   7 #include <stdio.h>
   8 #include <stdlib.h>
   9 #include <math.h>
  10 #include "complex.h"
  11 /*#include <libraries/mathieeedp.h>*/
  12
  13 /*extern struct Library * MathIeeeDoubTransBase;*/
  14
  15 extern const Complex zero, eye;
  16
  17 const Complex   one = {1.0, 0.0},
  18                 minuseye = {0.0, -1.0};
  19
  20 int precision = 12;             /* number of decimal places to print */
  21                                 /* (when they exist) */
  22
  23 void cprin(fp, prefix, suffix, z)       /* print a complex number to file fp */
  24         FILE *fp;
  25         char *prefix, *suffix;
  26         Complex z;
  27 {
  28         fputs(prefix, fp);
  29
  30         if (z.imag == 0.0)
  31                 fprintf(fp, "%.*g", precision, z.real);
  32         else if (z.real == 0.0)
  33                 fprintf(fp, "%.*g i", precision, z.imag);
  34         else
  35                 fprintf(fp, "%.*g %c %.*g i",
  36                         precision, z.real, sign(z.imag), precision, fabs(z.imag));
  37
  38         fputs(suffix, fp);
  39 }
  40
  41 double Re(z)            /* real part of complex number */
  42         Complex z;
  43 {
  44         return z.real;
  45 }
  46
  47 double Im(z)            /* imaginary part of complex number */
  48         Complex z;
  49 {
  50         return z.imag;
  51 }
  52
  53 #define marg(z) atan((z).imag / (z).real)       /* macro arg */
  54
  55 double arg(z)           /* argument of complex number, in range (-PI,PI] */
  56         Complex z;
  57 {
  58         return atan(z.imag / z.real);
  59 }
  60
  61 double norm(z)          /* norm of a complex number */
  62         Complex z;
  63 {
  64         return sqr(z.real) + sqr(z.imag);
  65 }
  66
  67 double _cabs(z)         /* absolute value of a complex number */
  68         Complex z;
  69 {
  70         return sqrt(norm(z));
  71 }
  72
  73 Complex cadd(w,z)       /* complex addition */
  74         Complex w,z;
  75 {
  76         Complex rv;
  77
  78         rv.real = w.real + z.real;
  79         rv.imag = w.imag + z.imag;
  80
  81         return rv;
  82 }
  83
  84 Complex csub(w,z)       /* complex subtraction */
  85         Complex w,z;
  86 {
  87         Complex rv;
  88
  89         rv.real = w.real - z.real;
  90         rv.imag = w.imag - z.imag;
  91
  92         return rv;
  93 }
  94
  95 Complex cmul(w,z)       /* complex multiplication */
  96         Complex w,z;
  97 {
  98         Complex rv;
  99
 100         rv.real = w.real*z.real - w.imag*z.imag;
 101         rv.imag = w.real*z.imag + w.imag*z.real;
 102
 103         return rv;
 104 }
 105
 106 Complex cdiv(w,z)       /* complex division */
 107         Complex w,z;
 108 {
 109         Complex rv;
 110         double temp = sqr(z.real)+sqr(z.imag);
 111
 112         if (temp == 0.0)
 113                 execerror("division by zero", NULL);
 114
 115         rv.real = (w.real*z.real + w.imag*z.imag)/temp;
 116         rv.imag = (w.imag*z.real - w.real*z.imag)/temp;
 117
 118         return rv;
 119 }
 120
 121 Complex cneg(z)         /* complex negation */
 122         Complex z;
 123 {
 124         z.real = -z.real;
 125         z.imag = -z.imag;
 126
 127         return z;
 128 }
 129
 130 Complex csqr(z)         /* complex square, w^2 */
 131         Complex z;
 132 {
 133         Complex rv;
 134
 135         if (z.imag == 0.0)
 136         {
 137                 z.real *= z.real;
 138                 return z;
 139         }
 140         rv.real = sqr(z.real) - sqr(z.imag);
 141         rv.imag = 2*z.real*z.imag;
 142
 143         return rv;
 144 }
 145
 146 Complex csqrt(z)        /* complex square-root */
 147         Complex z;
 148 {
 149         Complex rv;
 150         double temp = _cabs(z);
 151
 152         rv.real = Sqrt((temp + z.real)*0.5);
 153         rv.imag = Sqrt((temp - z.real)*0.5);
 154
 155         return rv;
 156 }
 157
 158 Complex conj(z)         /* conjugate of w */
 159         Complex z;
 160 {
 161         z.imag = -z.imag;
 162
 163         return z;
 164 }
 165
 166 Complex cexp(z)         /* complex exponential function e^z */
 167         Complex z;
 168 {
 169         double temp = exp(z.real);
 170
 171         if (z.imag == 0.0)
 172                 z.real = temp;
 173         else
 174         {
 175                 z.real = temp*cos(z.imag);
 176                 z.imag = temp*sin(z.imag);
 177         }
 178         return z;
 179 }
 180
 181 Complex clog(z)         /* complex natural logarithm */
 182         Complex z;
 183 {
 184         Complex rv;
 185
 186         rv.real = Log(norm(z))*0.5;
 187         rv.imag = marg(z);
 188
 189         return rv;
 190 }
 191
 192 Complex cpow(w,z)       /* complex exponential, w^z */
 193         Complex w,z;
 194 {
 195         return cexp(cmul(z,clog(w)));
 196 }
 197
 198 Complex csin(z)         /* complex sine */
 199         Complex z;
 200 {
 201         if (z.imag == 0.0)
 202         {
 203                 z.real = sin(z.real);
 204                 return z;
 205         }
 206         else
 207         {
 208                 Complex rv;
 209
 210                 rv.real = sin(z.real)*cosh(z.imag);
 211                 rv.imag = cos(z.real)*sinh(z.imag);
 212
 213                 return rv;
 214         }
 215 }
 216
 217 Complex ccos(z)         /* complex cosine */
 218         Complex z;
 219 {
 220         if (z.imag == 0.0)
 221         {
 222                 z.real = cos(z.real);
 223                 return z;
 224         }
 225         else
 226         {
 227                 Complex rv;
 228
 229                 rv.real = cos(z.real)*cosh(z.imag);
 230                 rv.imag = -(sin(z.real)*sinh(z.imag));
 231
 232                 return rv;
 233         }
 234 }
 235
 236 Complex ctan(z)         /* complex tangent */
 237         Complex z;
 238 {
 239         if (z.imag == 0.0)
 240         {
 241                 z.real = tan(z.real);
 242                 return z;
 243         }
 244         else
 245                 return cdiv(csin(z),ccos(z));
 246 }
 247
 248 Complex casin(z)        /* complex arcsine - lazy version */
 249         Complex z;
 250 {
 251         /* asin z = -ilog(iz + sqrt(1-z^2)) */
 252
 253         if (abs(z.real) <= 1.0 && z.imag == 0.0)
 254         {
 255                 z.real = Asin(z.real);
 256                 return z;
 257         }
 258         else
 259                 return cmul(minuseye,clog(cadd(cmul(eye,z),csqrt(csub(one,csqr(z))))));
 260 }
 261
 262 Complex cacos(z)        /* complex arccosine - lazy version */
 263         Complex z;
 264 {
 265         /* acos z = -ilog(z + sqrt(z^2-1)) */
 266
 267         if (abs(z.real) <= 1.0 && z.imag == 0.0)
 268         {
 269                 z.real = Acos(z.real);
 270                 return z;
 271         }
 272         else
 273                 return cmul(minuseye,clog(cadd(z,csqrt(csub(csqr(z),one)))));
 274 }
 275
 276 Complex catan(z)        /* complex arctangent - not so lazy version */
 277         Complex z;
 278 {
 279         if (z.imag == 0.0)
 280                 z.real = atan(z.real);
 281         else
 282         {
 283                 Complex ctemp;
 284                 double temp = norm(z);
 285
 286                 ctemp.real = ctemp.imag = 1.0 / (1.0 + temp - 2.0 * z.imag);
 287                 ctemp.real *= 1.0 - temp;
 288                 ctemp.imag *= -2.0*z.real;
 289                 ctemp = clog(ctemp);
 290                 z.real = -0.5*ctemp.imag;
 291                 z.imag = 0.5*ctemp.real;
 292         }
 293
 294         return z;
 295 }
 296
 297 Complex csinh(z)        /* complex hyperbolic sine */
 298         Complex z;
 299 {
 300         Complex rv;
 301
 302         rv.real = cos(z.imag)*sinh(z.real);
 303         rv.imag = sin(z.imag)*cosh(z.real);
 304
 305         return rv;
 306 }
 307
 308 Complex ccosh(z)                /* complex hyperbolic cosine */
 309         Complex z;
 310 {
 311         Complex rv;
 312
 313         rv.real = cos(z.imag)*cosh(z.real);
 314         rv.imag = sin(z.imag)*sinh(z.real);
 315
 316         return rv;
 317 }
 318
 319 Complex ctanh(z)                /* complex tangent */
 320         Complex z;
 321 {
 322         return cdiv(csinh(z),ccosh(z));
 323 }