Return to pow.c CVS log

Up to [Research Unix] / researchv9 / cmd / cfront / libC / complex

Annotation of researchv9/cmd/cfront/libC/complex/pow.c, revision 1.1

1.1     ! root        1: 
        !             2: #include "complex.h"
        !             3: 
        !             4: complex pow(double base, complex z)
        !             5: /*
        !             6:        real to complex power: base**z.
        !             7: */
        !             8: {
        !             9:        complex y;
        !            10:        
        !            11:        if (base == 0) return y;        /* even for singularity */
        !            12: 
        !            13:        if (0 < base) {
        !            14:                double lb = log(base);
        !            15:                y.re = z.re * lb;
        !            16:                y.im = z.im * lb;
        !            17:                return exp(y);
        !            18:        }
        !            19:        
        !            20:        return pow(complex(base), z);   /* use complex power fct */
        !            21: }
        !            22: 
        !            23: 
        !            24: complex  pow(complex a, int n)
        !            25: /*
        !            26:        complex to integer power: a**n.
        !            27: */
        !            28: {
        !            29:        complex x; 
        !            30:        complex p = 1;
        !            31: 
        !            32:        if (n == 0) return p;
        !            33: 
        !            34:        if (n < 0) {
        !            35:                n = -n;
        !            36:                x = 1/a;
        !            37:        }
        !            38:        else
        !            39:                x = a;
        !            40: 
        !            41:        for( ; ; ) {
        !            42:                if(n & 01) {
        !            43:                        register double t = p.re * x.re - p.im * x.im;
        !            44:                        p.im = p.re * x.im + p.im * x.re;
        !            45:                        p.re = t;
        !            46:                }
        !            47:                if(n >>= 1) {
        !            48:                        register double t = x.re * x.re - x.im * x.im;
        !            49:                        x.im = 2 * x.re * x.im;
        !            50:                        x.re = t;
        !            51:                }
        !            52:                else    break;
        !            53:        }
        !            54:        return p;
        !            55: }
        !            56: 
        !            57: complex pow(complex a, double b)
        !            58: /*
        !            59:        complex to real power: a**b.
        !            60: */
        !            61: {
        !            62:        register double logr = log( abs(a) );
        !            63:        register double logi = atan2(a.im, a.re);
        !            64:        register double x = exp( b*logr );
        !            65:        register double y = b * logi;
        !            66:        return complex(x*cos(y), x*sin(y));
        !            67: }
        !            68: 
        !            69: 
        !            70: complex pow(complex base, complex sup)
        !            71: /*
        !            72:        complex to complex power: base**sup.
        !            73: */
        !            74: {
        !            75:        complex result;
        !            76:        register double logr, logi;
        !            77:        register double xx, yy;
        !            78:        double a = abs(base);
        !            79: 
        !            80:        if (a == 0) return result;
        !            81: 
        !            82:        logr = log( a );
        !            83:        logi = atan2(base.im, base.re);
        !            84: 
        !            85:        xx = exp( logr * sup.re - logi * sup.im );
        !            86:        yy = logr * sup.im + logi * sup.re;
        !            87: 
        !            88:        result.re = xx * cos(yy);
        !            89:        result.im = xx * sin(yy);
        !            90: 
        !            91:        return result;
        !            92: }

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