/src/rt/complex.c
C | 113 lines | 91 code | 10 blank | 12 comment | 14 complexity | d1047f96e46231aacc71b3f7587b5902 MD5 | raw file
1/** 2 * Implementation of complex number support routines. 3 * 4 * Copyright: Copyright Digital Mars 2000 - 2010. 5 * License: <a href="http://www.boost.org/LICENSE_1_0.txt">Boost License 1.0</a>. 6 * Authors: Walter Bright 7 */ 8 9/* Copyright Digital Mars 2000 - 2010. 10 * Distributed under the Boost Software License, Version 1.0. 11 * (See accompanying file LICENSE or copy at 12 * http://www.boost.org/LICENSE_1_0.txt) 13 */ 14#include <math.h> 15 16typedef struct Complex 17{ 18 long double re; 19 long double im; 20} Complex; 21 22Complex _complex_div(Complex x, Complex y) 23{ 24 Complex q; 25 long double r; 26 long double den; 27 28 if (fabs(y.re) < fabs(y.im)) 29 { 30 r = y.re / y.im; 31 den = y.im + r * y.re; 32 q.re = (x.re * r + x.im) / den; 33 q.im = (x.im * r - x.re) / den; 34 } 35 else 36 { 37 r = y.im / y.re; 38 den = y.re + r * y.im; 39 q.re = (x.re + r * x.im) / den; 40 q.im = (x.im - r * x.re) / den; 41 } 42 return q; 43} 44 45Complex _complex_mul(Complex x, Complex y) 46{ 47 Complex p; 48 49 p.re = x.re * y.re - x.im * y.im; 50 p.im = x.im * y.re + x.re * y.im; 51 return p; 52} 53 54long double _complex_abs(Complex z) 55{ 56 long double x,y,ans,temp; 57 58 x = fabs(z.re); 59 y = fabs(z.im); 60 if (x == 0) 61 ans = y; 62 else if (y == 0) 63 ans = x; 64 else if (x > y) 65 { 66 temp = y / x; 67 ans = x * sqrt(1 + temp * temp); 68 } 69 else 70 { 71 temp = x / y; 72 ans = y * sqrt(1 + temp * temp); 73 } 74 return ans; 75} 76 77Complex _complex_sqrt(Complex z) 78{ 79 Complex c; 80 long double x,y,w,r; 81 82 if (z.re == 0 && z.im == 0) 83 { 84 c.re = 0; 85 c.im = 0; 86 } 87 else 88 { 89 x = fabs(z.re); 90 y = fabs(z.im); 91 if (x >= y) 92 { 93 r = y / x; 94 w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r))); 95 } 96 else 97 { 98 r = x / y; 99 w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r))); 100 } 101 if (z.re >= 0) 102 { 103 c.re = w; 104 c.im = z.im / (w + w); 105 } 106 else 107 { 108 c.im = (z.im >= 0) ? w : -w; 109 c.re = z.im / (c.im + c.im); 110 } 111 } 112 return c; 113}