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