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

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.