home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / unix / volume26 / calc / part05 / comfunc.c next >
C/C++ Source or Header  |  1992-05-09  |  10KB  |  517 lines
  1. /*
  2.  * Copyright (c) 1992 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Extended precision complex arithmetic non-primitive routines
  7.  */
  8.  
  9. #include "calc.h"
  10.  
  11.  
  12. /*
  13.  * Round a complex number to the specified number of decimal places.
  14.  * This simply means to round each of the components of the number.
  15.  * Zero decimal places means round to the nearest complex integer.
  16.  */
  17. COMPLEX *
  18. cround(c, places)
  19.     COMPLEX *c;
  20.     long places;
  21. {
  22.     COMPLEX *res;        /* result */
  23.  
  24.     res = comalloc();
  25.     res->real = qround(c->real, places);
  26.     res->imag = qround(c->imag, places);
  27.     return res;
  28. }
  29.  
  30.  
  31. /*
  32.  * Round a complex number to the specified number of binary decimal places.
  33.  * This simply means to round each of the components of the number.
  34.  * Zero binary places means round to the nearest complex integer.
  35.  */
  36. COMPLEX *
  37. cbround(c, places)
  38.     COMPLEX *c;
  39.     long places;
  40. {
  41.     COMPLEX *res;        /* result */
  42.  
  43.     res = comalloc();
  44.     res->real = qbround(c->real, places);
  45.     res->imag = qbround(c->imag, places);
  46.     return res;
  47. }
  48.  
  49.  
  50. /*
  51.  * Compute the result of raising a complex number to an integer power.
  52.  */
  53. COMPLEX *
  54. cpowi(c, q)
  55.     COMPLEX *c;        /* complex number to be raised */
  56.     NUMBER *q;        /* power to raise it to */
  57. {
  58.     COMPLEX *tmp, *res;    /* temporary values */
  59.     long power;        /* power to raise to */
  60.     unsigned long bit;    /* current bit value */
  61.     int sign;
  62.  
  63.     if (qisfrac(q))
  64.         error("Raising number to non-integral power");
  65.     if (isbig(q->num))
  66.         error("Raising number to very large power");
  67.     power = (istiny(q->num) ? z1tol(q->num) : z2tol(q->num));
  68.     if (ciszero(c) && (power == 0))
  69.         error("Raising zero to zeroth power");
  70.     sign = 1;
  71.     if (qisneg(q))
  72.         sign = -1;
  73.     /*
  74.      * Handle some low powers specially
  75.      */
  76.     if (power <= 4) {
  77.         switch ((int) (power * sign)) {
  78.             case 0:
  79.                 return clink(&_cone_);
  80.             case 1:
  81.                 return clink(c);
  82.             case -1:
  83.                 return cinv(c);
  84.             case 2:
  85.                 return csquare(c);
  86.             case -2:
  87.                 tmp = csquare(c);
  88.                 res = cinv(tmp);
  89.                 comfree(tmp);
  90.                 return res;
  91.             case 3:
  92.                 tmp = csquare(c);
  93.                 res = cmul(c, tmp);
  94.                 comfree(tmp);
  95.                 return res;
  96.             case 4:
  97.                 tmp = csquare(c);
  98.                 res = csquare(tmp);
  99.                 comfree(tmp);
  100.                 return res;
  101.         }
  102.     }
  103.     /*
  104.      * Compute the power by squaring and multiplying.
  105.      * This uses the left to right method of power raising.
  106.      */
  107.     bit = TOPFULL;
  108.     while ((bit & power) == 0)
  109.         bit >>= 1L;
  110.     bit >>= 1L;
  111.     res = csquare(c);
  112.     if (bit & power) {
  113.         tmp = cmul(res, c);
  114.         comfree(res);
  115.         res = tmp;
  116.     }
  117.     bit >>= 1L;
  118.     while (bit) {
  119.         tmp = csquare(res);
  120.         comfree(res);
  121.         res = tmp;
  122.         if (bit & power) {
  123.             tmp = cmul(res, c);
  124.             comfree(res);
  125.             res = tmp;
  126.         }
  127.         bit >>= 1L;
  128.     }
  129.     if (sign < 0) {
  130.         tmp = cinv(res);
  131.         comfree(res);
  132.         res = tmp;
  133.     }
  134.     return res;
  135. }
  136.  
  137.  
  138. /*
  139.  * Calculate the square root of a complex number, with each component
  140.  * within the specified error. This uses the formula:
  141.  *    sqrt(a+bi) = sqrt((a+sqrt(a^2+b^2))/2) + sqrt((-a+sqrt(a^2+b^2))/2)i.
  142.  */
  143. COMPLEX *
  144. csqrt(c, epsilon)
  145.     COMPLEX *c;
  146.     NUMBER *epsilon;
  147. {
  148.     COMPLEX *r;
  149.     NUMBER *hypot, *tmp1, *tmp2;
  150.  
  151.     if (ciszero(c) || cisone(c))
  152.         return clink(c);
  153.     r = comalloc();
  154.     if (cisreal(c)) {
  155.         if (!qisneg(c->real)) {
  156.             r->real = qsqrt(c->real, epsilon);
  157.             return r;
  158.         }
  159.         tmp1 = qneg(c->real);
  160.         r->imag = qsqrt(tmp1, epsilon);
  161.         qfree(tmp1);
  162.         return r;
  163.     }
  164.     hypot = qhypot(c->real, c->imag, epsilon);
  165.     tmp1 = qadd(hypot, c->real);
  166.     tmp2 = qscale(tmp1, -1L);
  167.     qfree(tmp1);
  168.     r->real = qsqrt(tmp2, epsilon);
  169.     qfree(tmp2);
  170.     tmp1 = qsub(hypot, c->real);
  171.     qfree(hypot);
  172.     tmp2 = qscale(tmp1, -1L);
  173.     qfree(tmp1);
  174.     tmp1 = qsqrt(tmp2, epsilon);
  175.     qfree(tmp2);
  176.     if (qisneg(c->imag)) {
  177.         tmp2 = qneg(tmp1);
  178.         qfree(tmp1);
  179.         tmp1 = tmp2;
  180.     }
  181.     r->imag = tmp1;
  182.     return r;
  183. }
  184.  
  185.  
  186. /*
  187.  * Take the Nth root of a complex number, where N is a positive integer.
  188.  * Each component of the result is within the specified error.
  189.  */
  190. COMPLEX *
  191. croot(c, q, epsilon)
  192.     COMPLEX *c;
  193.     NUMBER *q, *epsilon;
  194. {
  195.     COMPLEX *r;
  196.     NUMBER *a2pb2, *root, *tmp1, *tmp2, *epsilon2;
  197.  
  198.     if (qisneg(q) || qiszero(q) || qisfrac(q))
  199.         error("Taking bad root of complex number");
  200.     if (cisone(c) || qisone(q))
  201.         return clink(c);
  202.     if (qistwo(q))
  203.         return csqrt(c, epsilon);
  204.     r = comalloc();
  205.     if (cisreal(c) && !qisneg(c->real)) {
  206.         r->real = qroot(c->real, q, epsilon);
  207.         return r;
  208.     }
  209.     /*
  210.      * Calculate the root using the formula:
  211.      *    croot(a + bi, n) =
  212.      *        cpolar(qroot(a^2 + b^2, 2 * n), qatan2(b, a) / n).
  213.      */
  214.     epsilon2 = qscale(epsilon, -8L);
  215.     tmp1 = qsquare(c->real);
  216.     tmp2 = qsquare(c->imag);
  217.     a2pb2 = qadd(tmp1, tmp2);
  218.     qfree(tmp1);
  219.     qfree(tmp2);
  220.     tmp1 = qscale(q, 1L);
  221.     root = qroot(a2pb2, tmp1, epsilon2);
  222.     qfree(a2pb2);
  223.     qfree(tmp1);
  224.     tmp1 = qatan2(c->imag, c->real, epsilon2);
  225.     qfree(epsilon2);
  226.     tmp2 = qdiv(tmp1, q);
  227.     qfree(tmp1);
  228.     r = cpolar(root, tmp2, epsilon);
  229.     qfree(root);
  230.     qfree(tmp2);
  231.     return r;
  232. }
  233.  
  234.  
  235. /*
  236.  * Calculate the complex exponential function to the desired accuracy.
  237.  * We use the formula:
  238.  *    exp(a + bi) = exp(a) * (cos(b) + i * sin(b)).
  239.  */
  240. COMPLEX *
  241. cexp(c, epsilon)
  242.     COMPLEX *c;
  243.     NUMBER *epsilon;
  244. {
  245.     COMPLEX *r;
  246.     NUMBER *tmp1, *tmp2, *epsilon2;
  247.  
  248.     if (ciszero(c))
  249.         return clink(&_cone_);
  250.     r = comalloc();
  251.     if (cisreal(c)) {
  252.         r->real = qexp(c->real, epsilon);
  253.         return r;
  254.     }
  255.     epsilon2 = qscale(epsilon, -2L);
  256.     r->real = qcos(c->imag, epsilon2);
  257.     r->imag = qlegtoleg(r->real, epsilon2, _sinisneg_);
  258.     if (qiszero(c->real)) {
  259.         qfree(epsilon2);
  260.         return r;
  261.     }
  262.     tmp1 = qexp(c->real, epsilon2);
  263.     qfree(epsilon2);
  264.     tmp2 = qmul(r->real, tmp1);
  265.     qfree(r->real);
  266.     r->real = tmp2;
  267.     tmp2 = qmul(r->imag, tmp1);
  268.     qfree(r->imag);
  269.     qfree(tmp1);
  270.     r->imag = tmp2;
  271.     return r;
  272. }
  273.  
  274.  
  275. /*
  276.  * Calculate the natural logarithm of a complex number within the specified
  277.  * error.  We use the formula:
  278.  *    ln(a + bi) = ln(a^2 + b^2) / 2 + i * atan2(b, a).
  279.  */
  280. COMPLEX *
  281. cln(c, epsilon)
  282.     COMPLEX *c;
  283.     NUMBER *epsilon;
  284. {
  285.     COMPLEX *r;
  286.     NUMBER *a2b2, *tmp1, *tmp2;
  287.  
  288.     if (ciszero(c))
  289.         error("Logarithm of zero");
  290.     if (cisone(c))
  291.         return clink(&_czero_);
  292.     r = comalloc();
  293.     if (cisreal(c) && !qisneg(c->real)) {
  294.         r->real = qln(c->real, epsilon);
  295.         return r;
  296.     }
  297.     tmp1 = qsquare(c->real);
  298.     tmp2 = qsquare(c->imag);
  299.     a2b2 = qadd(tmp1, tmp2);
  300.     qfree(tmp1);
  301.     qfree(tmp2);
  302.     tmp1 = qln(a2b2, epsilon);
  303.     qfree(a2b2);
  304.     r->real = qscale(tmp1, -1L);
  305.     qfree(tmp1);
  306.     r->imag = qatan2(c->imag, c->real, epsilon);
  307.     return r;
  308. }
  309.  
  310.  
  311. /*
  312.  * Calculate the complex cosine within the specified accuracy.
  313.  * This uses the formula:
  314.  *    cos(a + bi) = cos(a) * cosh(b) - sin(a) * sinh(b) * i.
  315.  */
  316. COMPLEX *
  317. ccos(c, epsilon)
  318.     COMPLEX *c;
  319.     NUMBER *epsilon;
  320. {
  321.     COMPLEX *r;
  322.     NUMBER *cosval, *coshval, *tmp1, *tmp2, *tmp3, *epsilon2;
  323.     int negimag;
  324.  
  325.     if (ciszero(c))
  326.         return clink(&_cone_);
  327.     r = comalloc();
  328.     if (cisreal(c)) {
  329.         r->real = qcos(c->real, epsilon);
  330.         return r;
  331.     }
  332.     if (qiszero(c->real)) {
  333.         r->real = qcosh(c->imag, epsilon);
  334.         return r;
  335.     }
  336.     epsilon2 = qscale(epsilon, -2L);
  337.     coshval = qcosh(c->imag, epsilon2);
  338.     cosval = qcos(c->real, epsilon2);
  339.     negimag = !_sinisneg_;
  340.     if (qisneg(c->imag))
  341.         negimag = !negimag;
  342.     r->real = qmul(cosval, coshval);
  343.     /*
  344.      * Calculate the imaginary part using the formula:
  345.      *    sin(a) * sinh(b) = sqrt((1 - a^2) * (b^2 - 1)).
  346.      */
  347.     tmp1 = qsquare(cosval);
  348.     qfree(cosval);
  349.     tmp2 = qdec(tmp1);
  350.     qfree(tmp1);
  351.     tmp1 = qneg(tmp2);
  352.     qfree(tmp2);
  353.     tmp2 = qsquare(coshval);
  354.     qfree(coshval);
  355.     tmp3 = qdec(tmp2);
  356.     qfree(tmp2);
  357.     tmp2 = qmul(tmp1, tmp3);
  358.     qfree(tmp1);
  359.     qfree(tmp3);
  360.     r->imag = qsqrt(tmp2, epsilon2);
  361.     qfree(tmp2);
  362.     qfree(epsilon2);
  363.     if (negimag) {
  364.         tmp1 = qneg(r->imag);
  365.         qfree(r->imag);
  366.         r->imag = tmp1;
  367.     }
  368.     return r;
  369. }
  370.  
  371.  
  372. /*
  373.  * Calculate the complex sine within the specified accuracy.
  374.  * This uses the formula:
  375.  *    sin(a + bi) = sin(a) * cosh(b) + cos(a) * sinh(b) * i.
  376.  */
  377. COMPLEX *
  378. csin(c, epsilon)
  379.     COMPLEX *c;
  380.     NUMBER *epsilon;
  381. {
  382.     COMPLEX *r;
  383.  
  384.     NUMBER *cosval, *coshval, *tmp1, *tmp2, *epsilon2;
  385.  
  386.     if (ciszero(c))
  387.         return clink(&_czero_);
  388.     r = comalloc();
  389.     if (cisreal(c)) {
  390.         r->real = qsin(c->real, epsilon);
  391.         return r;
  392.     }
  393.     if (qiszero(c->real)) {
  394.         r->imag = qsinh(c->imag, epsilon);
  395.         return r;
  396.     }
  397.     epsilon2 = qscale(epsilon, -2L);
  398.     coshval = qcosh(c->imag, epsilon2);
  399.     cosval = qcos(c->real, epsilon2);
  400.     tmp1 = qlegtoleg(cosval, epsilon2, _sinisneg_);
  401.     r->real = qmul(tmp1, coshval);
  402.     qfree(tmp1);
  403.     tmp1 = qsquare(coshval);
  404.     qfree(coshval);
  405.     tmp2 = qdec(tmp1);
  406.     qfree(tmp1);
  407.     tmp1 = qsqrt(tmp2, epsilon2);
  408.     qfree(tmp2);
  409.     r->imag = qmul(tmp1, cosval);
  410.     qfree(tmp1);
  411.     qfree(cosval);
  412.     if (qisneg(c->imag)) {
  413.         tmp1 = qneg(r->imag);
  414.         qfree(r->imag);
  415.         r->imag = tmp1;
  416.     }
  417.     return r;
  418. }
  419.  
  420.  
  421. /*
  422.  * Convert a number from polar coordinates to normal complex number form
  423.  * within the specified accuracy.  This produces the value:
  424.  *    q1 * cos(q2) + q1 * sin(q2) * i.
  425.  */
  426. COMPLEX *
  427. cpolar(q1, q2, epsilon)
  428.     NUMBER *q1, *q2, *epsilon;
  429. {
  430.     COMPLEX *r;
  431.     NUMBER *tmp, *epsilon2;
  432.     long scale;
  433.  
  434.     r = comalloc();
  435.     if (qiszero(q1) || qiszero(q2)) {
  436.         r->real = qlink(q1);
  437.         return r;
  438.     }
  439.     epsilon2 = epsilon;
  440.     if (!qisunit(q1)) {
  441.         scale = zhighbit(q1->num) - zhighbit(q1->den) + 1;
  442.         if (scale > 0)
  443.             epsilon2 = qscale(epsilon, -scale);
  444.     }
  445.     r->real = qcos(q2, epsilon2);
  446.     r->imag = qlegtoleg(r->real, epsilon2, _sinisneg_);
  447.     if (epsilon2 != epsilon)
  448.         qfree(epsilon2);
  449.     if (qisone(q1))
  450.         return r;
  451.     tmp = qmul(r->real, q1);
  452.     qfree(r->real);
  453.     r->real = tmp;
  454.     tmp = qmul(r->imag, q1);
  455.     qfree(r->imag);
  456.     r->imag = tmp;
  457.     return r;
  458. }
  459.  
  460.  
  461. /*
  462.  * Raise one complex number to the power of another one to within the
  463.  * specified error.
  464.  */
  465. COMPLEX *
  466. cpower(c1, c2, epsilon)
  467.     COMPLEX *c1, *c2;
  468.     NUMBER *epsilon;
  469. {
  470.     COMPLEX *tmp1, *tmp2;
  471.     NUMBER *epsilon2;
  472.  
  473.     if (cisreal(c2) && qisint(c2->real))
  474.         return cpowi(c1, c2->real);
  475.     if (cisone(c1) || ciszero(c1))
  476.         return clink(c1);
  477.     epsilon2 = qscale(epsilon, -4L);
  478.     tmp1 = cln(c1, epsilon2);
  479.     tmp2 = cmul(tmp1, c2);
  480.     comfree(tmp1);
  481.     tmp1 = cexp(tmp2, epsilon);
  482.     comfree(tmp2);
  483.     qfree(epsilon2);
  484.     return tmp1;
  485. }
  486.  
  487.  
  488. /*
  489.  * Print a complex number.
  490.  */
  491. void
  492. comprint(c)
  493.     COMPLEX *c;
  494. {
  495.     NUMBER qtmp;
  496.  
  497.     if (_outmode_ == MODE_FRAC) {
  498.         cprintfr(c);
  499.         return;
  500.     }
  501.     if (!qiszero(c->real) || qiszero(c->imag))
  502.         qprintnum(c->real, MODE_DEFAULT);
  503.     qtmp = c->imag[0];
  504.     if (qiszero(&qtmp))
  505.         return;
  506.     if (!qiszero(c->real) && !qisneg(&qtmp))
  507.         math_chr('+');
  508.     if (qisneg(&qtmp)) {
  509.         math_chr('-');
  510.         qtmp.num.sign = 0;
  511.     }
  512.     qprintnum(&qtmp, MODE_DEFAULT);
  513.     math_chr('i');
  514. }
  515.  
  516. /* END CODE */
  517.