/src/rt/complex.c

http://github.com/AlexeyProkhin/druntime · 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. /* Copyright Digital Mars 2000 - 2010.
  9. * Distributed under the Boost Software License, Version 1.0.
  10. * (See accompanying file LICENSE or copy at
  11. * http://www.boost.org/LICENSE_1_0.txt)
  12. */
  13. #include <math.h>
  14. typedef struct Complex
  15. {
  16. long double re;
  17. long double im;
  18. } Complex;
  19. Complex _complex_div(Complex x, Complex y)
  20. {
  21. Complex q;
  22. long double r;
  23. long double den;
  24. if (fabs(y.re) < fabs(y.im))
  25. {
  26. r = y.re / y.im;
  27. den = y.im + r * y.re;
  28. q.re = (x.re * r + x.im) / den;
  29. q.im = (x.im * r - x.re) / den;
  30. }
  31. else
  32. {
  33. r = y.im / y.re;
  34. den = y.re + r * y.im;
  35. q.re = (x.re + r * x.im) / den;
  36. q.im = (x.im - r * x.re) / den;
  37. }
  38. return q;
  39. }
  40. Complex _complex_mul(Complex x, Complex y)
  41. {
  42. Complex p;
  43. p.re = x.re * y.re - x.im * y.im;
  44. p.im = x.im * y.re + x.re * y.im;
  45. return p;
  46. }
  47. long double _complex_abs(Complex z)
  48. {
  49. long double x,y,ans,temp;
  50. x = fabs(z.re);
  51. y = fabs(z.im);
  52. if (x == 0)
  53. ans = y;
  54. else if (y == 0)
  55. ans = x;
  56. else if (x > y)
  57. {
  58. temp = y / x;
  59. ans = x * sqrt(1 + temp * temp);
  60. }
  61. else
  62. {
  63. temp = x / y;
  64. ans = y * sqrt(1 + temp * temp);
  65. }
  66. return ans;
  67. }
  68. Complex _complex_sqrt(Complex z)
  69. {
  70. Complex c;
  71. long double x,y,w,r;
  72. if (z.re == 0 && z.im == 0)
  73. {
  74. c.re = 0;
  75. c.im = 0;
  76. }
  77. else
  78. {
  79. x = fabs(z.re);
  80. y = fabs(z.im);
  81. if (x >= y)
  82. {
  83. r = y / x;
  84. w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
  85. }
  86. else
  87. {
  88. r = x / y;
  89. w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
  90. }
  91. if (z.re >= 0)
  92. {
  93. c.re = w;
  94. c.im = z.im / (w + w);
  95. }
  96. else
  97. {
  98. c.im = (z.im >= 0) ? w : -w;
  99. c.re = z.im / (c.im + c.im);
  100. }
  101. }
  102. return c;
  103. }