Return to log.c CVS log

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

Annotation of researchv10no/cmd/cfront/libC/complex/log.c, revision 1.1.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.