Return to complex.c CVS log

Up to [Research Unix] / researchv10no / cmd / map / libmap

Annotation of researchv10no/cmd/map/libmap/complex.c, revision 1.1.1.1

1.1       root        1: #include "map.h"
                      2: 
                      3: /*complex divide, defensive against overflow from
                      4:  *     * and /, but not from + and -
                      5:  *     assumes underflow yields 0.0
                      6:  *     uses identities:
                      7:  *     (a + bi)/(c + di) = ((a + bd/c) + (b - ad/c)i)/(c + dd/c)
                      8:  *     (a + bi)/(c + di) = (b - ai)/(d - ci)
                      9: */
                     10: void
                     11: cdiv(double a, double b, double c, double d, double *u, double *v)
                     12: {
                     13:        double r,t;
                     14:        if(fabs(c)<fabs(d)) {
                     15:                t = -c; c = d; d = t;
                     16:                t = -a; a = b; b = t;
                     17:        }
                     18:        r = d/c;
                     19:        t = c + r*d;
                     20:        *u = (a + r*b)/t;
                     21:        *v = (b - r*a)/t;
                     22: }
                     23: 
                     24: void
                     25: cmul(double c1, double c2, double d1, double d2, double *e1, double *e2)
                     26: {
                     27:        *e1 = c1*d1 - c2*d2;
                     28:        *e2 = c1*d2 + c2*d1;
                     29: }
                     30: 
                     31: void
                     32: csq(double c1, double c2, double *e1, double *e2)
                     33: {
                     34:        *e1 = c1*c1 - c2*c2;
                     35:        *e2 = c1*c2*2;
                     36: }
                     37: 
                     38: /* complex square root
                     39:  *     assumes underflow yields 0.0
                     40:  *     uses these identities:
                     41:  *     sqrt(x+_iy) = sqrt(r(cos(t)+_isin(t))
                     42:  *                = sqrt(r)(cos(t/2)+_isin(t/2))
                     43:  *     cos(t/2) = sin(t)/2sin(t/2) = sqrt((1+cos(t)/2)
                     44:  *     sin(t/2) = sin(t)/2cos(t/2) = sqrt((1-cos(t)/2)
                     45: */
                     46: void
                     47: csqrt(double c1, double c2, double *e1, double *e2)
                     48: {
                     49:        double r,s;
                     50:        double x,y;
                     51:        x = fabs(c1);
                     52:        y = fabs(c2);
                     53:        if(x>=y) {
                     54:                if(x==0) {
                     55:                        *e1 = *e2 = 0;
                     56:                        return;
                     57:                }
                     58:                r = x;
                     59:                s = y/x;
                     60:        } else {
                     61:                r = y;
                     62:                s = x/y;
                     63:        }
                     64:        r *= sqrt(1+ s*s);
                     65:        if(c1>0) {
                     66:                *e1 = sqrt((r+c1)/2);
                     67:                *e2 = c2/(2* *e1);
                     68:        } else {
                     69:                *e2 = sqrt((r-c1)/2);
                     70:                if(c2<0)
                     71:                        *e2 = -*e2;
                     72:                *e1 = c2/(2* *e2);
                     73:        }
                     74: }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.