Return to complex.c CVS log

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

Annotation of researchv8dc/cmd/map/libmap/complex.c, revision 1.1

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

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