Return to log.c CVS log

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

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

1.1     ! root        1: 
        !             2: #include "complex.h"
        !             3: #include "const.h"
        !             4: 
        !             5: #define        LOGWILD -1000
        !             6: #define        LOGDANGER       1e18
        !             7: #define PERIL(t) (t > LOGDANGER || (t < 1/LOGDANGER && t != 0) )
        !             8: 
        !             9: complex log(complex z)
        !            10: /*
        !            11:        The complex natural logarithm of "z".
        !            12:        If z = 0, then the answer LOGWILD + 0*i is returned.
        !            13: 
        !            14:        Stu Feldman says that the peril tests for the following function
        !            15:        are "acceptable for now", but certain things like
        !            16:        complex variables outside the over/underflow range 
        !            17:        will cause floating exceptions.
        !            18: */
        !            19: {
        !            20:        complex answer; 
        !            21:        double  partial;
        !            22: 
        !            23:        if ( z.re == 0 && z.im == 0) {
        !            24:                complex_error(C_LOG_0,0);
        !            25:                answer.re = LOGWILD;
        !            26:        }
        !            27:        else {
        !            28:        /*
        !            29:                 Check for (over/under)flow, and fixup if necessary.
        !            30:        */
        !            31:                double x = ABS(z.re);
        !            32:                double y = ABS(z.im);
        !            33:     
        !            34:                if ( x>y && PERIL(x) ) {
        !            35:                        z.im /=x;
        !            36:                        z.re /= x;  /* z.re is replaced by 1 or -1 */
        !            37:                        partial = log(x);
        !            38:                }
        !            39:                else if (PERIL(y)) {
        !            40:                        z.im /= y;  /* roles of re, im reversed from previous */
        !            41:                        z.re /= y;
        !            42:                        partial = log(y);
        !            43:                        }
        !            44:                else partial = 0;
        !            45:     
        !            46:        /*
        !            47:                z.re*z.re and z.im*z.im should not cause problems now.
        !            48:        */
        !            49:     
        !            50:                answer.im = atan2(z.im,z.re); 
        !            51:                answer.re = log(z.re*z.re + z.im*z.im)/2 + partial;
        !            52:        }
        !            53:        return answer;
        !            54: }

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