home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / amiga / misc / icalc.lzh / icalc / src / cmath.c < prev    next >
C/C++ Source or Header  |  1992-01-26  |  7KB  |  387 lines
  1. /*
  2. *    icalc - complex-expression parser
  3. *
  4. *    Math routines for complex-number expression parser.
  5. *    In the routines, rv is variable holding 'return value'.
  6. *
  7. *    (C) Martin W Scott, 1991.
  8. */
  9.  
  10. #include <math.h>
  11. #include "complex.h"
  12. #include "constant.h"
  13.  
  14. #define NOERRCHECK    /* domain checking doesn't work when using */
  15.             /* IEEE libraries... */
  16. #ifdef NOERRCHECK
  17. #define Log(x) log(x)
  18. #endif
  19.  
  20. #define Zero(x)    ((x) == 0.0)    /* may add epsilons at a later date... */
  21.  
  22. Complex Re(z)        /* real part of complex number */
  23.     Complex z;
  24. {
  25.     z.imag = 0.0;
  26.     return z;
  27. }
  28.  
  29. Complex Im(z)        /* imaginary part of complex number */
  30.     Complex z;
  31. {
  32.     z.real = z.imag;
  33.     z.imag = 0.0;
  34.     return z;
  35. }
  36.  
  37. #define realarg(z)    atan2((z).imag,(z).real)    /* macro arg */
  38. #define realnorm(z) ((z).real*(z).real+(z).imag*(z).imag)    /* real norm */
  39.  
  40. static double sgn(x)        /* sign of x */
  41.     double x;
  42. {
  43.     if (x > 0.0) return 1.0;
  44.     else if (x < 0.0) return -1.0;
  45.     else return 0.0;
  46. }
  47.  
  48. Complex arg(z)        /* argument of complex number, in range (-PI,PI] */
  49.     Complex z;
  50. {
  51.     z.real = realarg(z);
  52.     z.imag = 0.0;
  53.     return z;
  54. }
  55.  
  56. Complex norm(z)        /* norm of a complex number */
  57.     Complex z;
  58. {
  59.     z.real = sqr(z.real) + sqr(z.imag);
  60.     z.imag = 0.0;
  61.     return z;
  62. }
  63.  
  64. Complex cabs(z)        /* absolute value of a complex number */
  65.     Complex z;
  66. {
  67.     z.real = sqrt(realnorm(z));
  68.     z.imag = 0.0;
  69.     return z;
  70. }
  71.  
  72. Complex cadd(w,z)    /* complex addition */
  73.     Complex w,z;
  74. {
  75.     w.real = w.real + z.real;
  76.     w.imag = w.imag + z.imag;
  77.  
  78.     return w;
  79. }
  80.  
  81. Complex csub(w,z)    /* complex subtraction */
  82.     Complex w,z;
  83. {
  84.     w.real = w.real - z.real;
  85.     w.imag = w.imag - z.imag;
  86.  
  87.     return w;
  88. }
  89.  
  90. Complex cmul(w,z)    /* complex multiplication */
  91.     Complex w,z;
  92. {
  93.     Complex rv;
  94.  
  95.     rv.real = w.real*z.real - w.imag*z.imag;
  96.     rv.imag = w.real*z.imag + w.imag*z.real;
  97.  
  98.     return rv;
  99. }
  100.  
  101. Complex cdiv(w,z)    /* complex division */
  102.     Complex w,z;
  103. {
  104.     Complex rv;
  105.     double temp = sqr(z.real)+sqr(z.imag);
  106.  
  107.     if (Zero(temp))
  108.         execerror("division by zero", NULL);    
  109.  
  110.     rv.real = (w.real*z.real + w.imag*z.imag)/temp;
  111.     rv.imag = (w.imag*z.real - w.real*z.imag)/temp;
  112.  
  113.     return rv;
  114. }
  115.  
  116. Complex cneg(z)        /* complex negation */
  117.     Complex z;
  118. {
  119.     z.real = -z.real;
  120.     z.imag = -z.imag;
  121.  
  122.     return z;
  123. }
  124.  
  125. Complex csqr(z)        /* complex square, w^2 */
  126.     Complex z;
  127. {
  128.     Complex rv;
  129.  
  130.     if (Zero(z.imag))    /* small optimisation for reals */
  131.     {
  132.         z.real *= z.real;
  133.         return z;
  134.     }
  135.     rv.real = sqr(z.real) - sqr(z.imag);
  136.     rv.imag = 2*z.real*z.imag;
  137.  
  138.     return rv;
  139. }
  140.  
  141. Complex csqrt(z)    /* complex square-root */
  142.     Complex z;
  143. {
  144.     if (Zero(z.imag))    /* real */
  145.     {
  146.         if (z.real >= 0.0)
  147.             z.real = sqrt(z.real);
  148.         else
  149.         {
  150.             z.imag = sqrt(-z.real);
  151.             z.real = 0.0;
  152.         }
  153.     }
  154.     else    /* complex */
  155.     {
  156.         double absval, temp;
  157.  
  158.         absval = sqrt(realnorm(z));
  159.         temp = sqrt((absval + z.real)/2.0);
  160.         z.imag = sgn(z.imag)*sqrt((absval - z.real)/2.0);
  161.         z.real = temp;
  162.     }
  163.  
  164.     return z;
  165. }
  166.  
  167. Complex conj(z)        /* conjugate of w */
  168.     Complex z;
  169. {
  170.     z.imag = -z.imag;
  171.  
  172.     return z;
  173. }
  174.  
  175. Complex cexp(z)        /* complex exponential function e^z */
  176.     Complex z;
  177. {
  178.     double temp = exp(z.real);
  179.  
  180.     if (Zero(z.imag))    /* small optimisation for reals */
  181.         z.real = temp;
  182.     else
  183.     {
  184.         z.real = temp*cos(z.imag);
  185.         z.imag = temp*sin(z.imag);
  186.     }
  187.     return z;
  188. }
  189.  
  190. Complex clog(z)        /* complex natural logarithm */
  191.     Complex z;
  192. {
  193.     Complex rv;
  194.  
  195.     if (Zero(z.imag) && Zero(z.real))
  196.         execerror("domain error in logarithm", NULL);
  197.  
  198.     rv.real = Log(realnorm(z))*0.5;
  199.     rv.imag = realarg(z);
  200.  
  201.     return rv;
  202. }
  203.  
  204. Complex cpow(w,z)    /* complex exponential, w^z */
  205.     Complex w,z;
  206. {
  207.     return cexp(cmul(z,clog(w)));
  208. }
  209.  
  210. Complex csin(z)        /* complex sine */
  211.     Complex z;
  212. {
  213.     if (Zero(z.imag))    /* small optimisation for reals */
  214.     {
  215.         z.real = sin(z.real);
  216.         return z;
  217.     }
  218.     else
  219.     {
  220.         Complex rv;
  221.  
  222.         rv.real = sin(z.real)*cosh(z.imag);
  223.         rv.imag = cos(z.real)*sinh(z.imag);
  224.  
  225.         return rv;
  226.     }
  227. }
  228.  
  229. Complex ccos(z)        /* complex cosine */
  230.     Complex z;
  231. {
  232.     if (Zero(z.imag))    /* small optimisation for reals */
  233.     {
  234.         z.real = cos(z.real);
  235.         return z;
  236.     }
  237.     else
  238.     {
  239.         Complex rv;
  240.  
  241.         rv.real = cos(z.real)*cosh(z.imag);
  242.         rv.imag = -(sin(z.real)*sinh(z.imag));
  243.  
  244.         return rv;
  245.     }
  246. }
  247.  
  248. Complex ctan(z)        /* complex tangent */
  249.     Complex z;
  250. {
  251.     if (Zero(z.imag))    /* small optimisation for reals */
  252.     {
  253.         z.real = tan(z.real);
  254.         return z;
  255.     }
  256.     else
  257.         return cdiv(csin(z),ccos(z));
  258. }
  259.  
  260. Complex casin(z)    /* complex arcsine - lazy version */
  261.     Complex z;
  262. {
  263.     /* asin z = -ilog(iz + sqrt(1-z^2)) */
  264.  
  265.     /* PVR taken as follows:
  266.      *    u+iv = acos(z), then
  267.      *    -PI/2 <= u <= PI/2    when v >= 0,
  268.      *    -PI/2  < u  < PI/2    when v < 0.
  269.      * This seems most logical, though conventions vary...
  270.      */
  271.  
  272.     double multiplier = 1.0;   /* to get standard definition into PVR */
  273.  
  274.     if (Zero(z.imag))
  275.     {
  276.         if (abs(z.real) <= 1.0)        /* real method */
  277.         {
  278.             z.real = asin(z.real);
  279.             return z;
  280.         }
  281.         multiplier = -sgn(z.real);
  282.             /* to get into PVR */
  283.     }
  284.  
  285.     z = cmul(minuseye,clog(cadd(cmul(eye,z),csqrt(csub(one,csqr(z))))));
  286.     z.imag *= multiplier;
  287.  
  288.     return z;
  289. }
  290.  
  291. Complex cacos(z)    /* complex arccosine */
  292.     Complex z;    /* if someone knows a neater way, tell me! */
  293. {
  294.     /* acos z = -ilog(z + sqrt(z^2-1)), into PVR */
  295.  
  296.     /* PVR taken as follows:
  297.      *    u+iv = acos(z), then
  298.      *    0 <= u <= PI    when v >= 0,
  299.      *    0  < u  < PI    when v < 0.
  300.      * This seems most logical, though conventions vary...
  301.      */
  302.  
  303.     double multiplier;    /* to get standard definition into PVR */
  304.  
  305.     if (Zero(z.imag))        /* on real axis */
  306.     {
  307.         if (abs(z.real) <= 1.0)        /* real-valued */
  308.         {
  309.             z.real = acos(z.real);
  310.             return z;
  311.         }
  312.  
  313.         multiplier = -sgn(z.real);
  314.             /* because csqrt returns prinipal value */
  315.             /* (1,inf)  -> (0,inf)i */
  316.             /* (-inf,1) -> PI+(0,inf)i */
  317.  
  318.     }
  319.     else if (!Zero(z.real)) multiplier = sgn(z.real)*sgn(z.imag);
  320.     else multiplier = 1.0;
  321.             /* again to get in PVR */
  322.  
  323.     z = cmul(minuseye,clog(cadd(z,csqrt(csub(csqr(z),one)))));
  324.     z.real *= multiplier;
  325.     z.imag *= multiplier;
  326.  
  327.     return z;
  328. }
  329.  
  330.  
  331. Complex catan(z)    /* complex arctangent - not so lazy version */
  332.     Complex z;
  333. {
  334.     /* easier than other two - no ambiguity about which root to take 
  335.      * (there aren't any) and no boundaries to worry about...
  336.      *
  337.      * PVR as follows: u+iv = atan(z), then
  338.      *    -PI/2 < u < PI/2
  339.      *    -inf  < v < inf
  340.      */
  341.  
  342.     if (Zero(z.imag))    /* small optimisation for reals */
  343.         z.real = atan(z.real);
  344.     else
  345.     {
  346.         Complex ctemp;
  347.         double temp = realnorm(z);
  348.  
  349.         ctemp.real = ctemp.imag = 1.0 / (1.0 + temp - 2.0 * z.imag);
  350.         ctemp.real *= 1.0 - temp;
  351.         ctemp.imag *= -2.0*z.real;
  352.         ctemp = clog(ctemp);
  353.         z.real = -0.5*ctemp.imag;
  354.         z.imag = 0.5*ctemp.real;
  355.     }
  356.  
  357.     return z;
  358. }
  359.  
  360. Complex csinh(z)    /* complex hyperbolic sine */
  361.     Complex z;
  362. {
  363.     Complex rv;
  364.  
  365.     rv.real = cos(z.imag)*sinh(z.real);
  366.     rv.imag = sin(z.imag)*cosh(z.real);
  367.  
  368.     return rv;
  369. }
  370.  
  371. Complex ccosh(z)        /* complex hyperbolic cosine */
  372.     Complex z;
  373. {
  374.     Complex rv;
  375.  
  376.     rv.real = cos(z.imag)*cosh(z.real);
  377.     rv.imag = sin(z.imag)*sinh(z.real);
  378.  
  379.     return rv;
  380. }
  381.  
  382. Complex ctanh(z)        /* complex tangent */
  383.     Complex z;
  384. {
  385.     return cdiv(csinh(z),ccosh(z));
  386. }
  387.