/src/rt/complex.c

  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}