Return to pow.c CVS log

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

Annotation of researchv8dc/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:        register 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: //     register complex x, p = 1;
                     30: register complex x; 
                     31: register complex p = 1;
                     32:        if (n == 0) return p;
                     33: 
                     34:        if (n < 0) {
                     35:                n = -n;
                     36:                x = 1/a;
                     37:        }
                     38:        else    x = a;
                     39: 
                     40:        for( ; ; ) {
                     41:                if(n & 01) {
                     42:                        register double t = p.re * x.re - p.im * x.im;
                     43:                        p.im = p.re * x.im + p.im * x.re;
                     44:                        p.re = t;
                     45:                }
                     46:                if(n >>= 1) {
                     47:                        register double t = x.re * x.re - x.im * x.im;
                     48:                        x.im = 2 * x.re * x.im;
                     49:                        x.re = t;
                     50:                }
                     51:                else    break;
                     52:        }
                     53:        return p;
                     54: }
                     55: 
                     56: complex pow(complex a, double b)
                     57: /*
                     58:        complex to real power: a**b.
                     59: */
                     60: {
                     61:        register double logr = log( abs(a) );
                     62:        register double logi = atan2(a.im, a.re);
                     63:        register double x = exp( b*logr );
                     64:        register double y = b * logi;
                     65:        return complex(x*cos(y), x*sin(y));
                     66: }
                     67: 
                     68: 
                     69: complex pow(complex base, complex sup)
                     70: /*
                     71:        complex to complex power: base**sup.
                     72: */
                     73: {
                     74:        complex result;
                     75:        register double logr, logi;
                     76:        register double xx, yy;
                     77:        double a = abs(base);
                     78: 
                     79:        if (a == 0) return result;
                     80: 
                     81:        logr = log( a );
                     82:        logi = atan2(base.im, base.re);
                     83: 
                     84:        xx = exp( logr * sup.re - logi * sup.im );
                     85:        yy = logr * sup.im + logi * sup.re;
                     86: 
                     87:        result.re = xx * cos(yy);
                     88:        result.im = xx * sin(yy);
                     89: 
                     90:        return result;
                     91: }

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