Mercurial > projects > ldc
diff druntime/src/compiler/dmd/complex.c @ 1458:e0b2d67cfe7c
Added druntime (this should be removed once it works).
author | Robert Clipsham <robert@octarineparrot.com> |
---|---|
date | Tue, 02 Jun 2009 17:43:06 +0100 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/druntime/src/compiler/dmd/complex.c Tue Jun 02 17:43:06 2009 +0100 @@ -0,0 +1,112 @@ +/** + * Implementation of complex number support routines. + * + * Copyright: Copyright Digital Mars 2000 - 2009. + * License: <a href="http://www.boost.org/LICENSE_1_0.txt>Boost License 1.0</a>. + * Authors: Walter Bright + * + * Copyright Digital Mars 2000 - 2009. + * Distributed under the Boost Software License, Version 1.0. + * (See accompanying file LICENSE_1_0.txt or copy at + * http://www.boost.org/LICENSE_1_0.txt) + */ +#include <math.h> + +typedef struct Complex +{ + long double re; + long double im; +} Complex; + +Complex _complex_div(Complex x, Complex y) +{ + Complex q; + long double r; + long double den; + + if (fabs(y.re) < fabs(y.im)) + { + r = y.re / y.im; + den = y.im + r * y.re; + q.re = (x.re * r + x.im) / den; + q.im = (x.im * r - x.re) / den; + } + else + { + r = y.im / y.re; + den = y.re + r * y.im; + q.re = (x.re + r * x.im) / den; + q.im = (x.im - r * x.re) / den; + } + return q; +} + +Complex _complex_mul(Complex x, Complex y) +{ + Complex p; + + p.re = x.re * y.re - x.im * y.im; + p.im = x.im * y.re + x.re * y.im; + return p; +} + +long double _complex_abs(Complex z) +{ + long double x,y,ans,temp; + + x = fabs(z.re); + y = fabs(z.im); + if (x == 0) + ans = y; + else if (y == 0) + ans = x; + else if (x > y) + { + temp = y / x; + ans = x * sqrt(1 + temp * temp); + } + else + { + temp = x / y; + ans = y * sqrt(1 + temp * temp); + } + return ans; +} + +Complex _complex_sqrt(Complex z) +{ + Complex c; + long double x,y,w,r; + + if (z.re == 0 && z.im == 0) + { + c.re = 0; + c.im = 0; + } + else + { + x = fabs(z.re); + y = fabs(z.im); + if (x >= y) + { + r = y / x; + w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r))); + } + else + { + r = x / y; + w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r))); + } + if (z.re >= 0) + { + c.re = w; + c.im = z.im / (w + w); + } + else + { + c.im = (z.im >= 0) ? w : -w; + c.re = z.im / (c.im + c.im); + } + } + return c; +}