home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / stlpt453.zip / STLport-4.5.3 / src / complex.cpp < prev    next >
C/C++ Source or Header  |  2001-05-03  |  8KB  |  285 lines
  1. /*
  2.  * Copyright (c) 1999
  3.  * Silicon Graphics Computer Systems, Inc.
  4.  *
  5.  * Copyright (c) 1999 
  6.  * Boris Fomitchev
  7.  *
  8.  * This material is provided "as is", with absolutely no warranty expressed
  9.  * or implied. Any use is at your own risk.
  10.  *
  11.  * Permission to use or copy this software for any purpose is hereby granted 
  12.  * without fee, provided the above notices are retained on all copies.
  13.  * Permission to modify the code and to distribute modified code is granted,
  14.  * provided the above notices are retained, and a notice that the code was
  15.  * modified is included with the above copyright notice.
  16.  *
  17.  */ 
  18. # include "stlport_prefix.h"
  19. // Complex division and square roots.
  20.  
  21. #include "complex_impl.h"
  22.  
  23. _STLP_BEGIN_NAMESPACE
  24.  
  25. // Absolute value
  26. _STLP_DECLSPEC float _STLP_CALL abs(const complex<float>& __z)
  27. {
  28.   return _STLP_HYPOTF(__z._M_re, __z._M_im);
  29. }
  30.  
  31. _STLP_DECLSPEC double _STLP_CALL abs(const complex<double>& __z)
  32. {
  33.   return _STLP_HYPOT(__z._M_re, __z._M_im);
  34. }
  35.  
  36. #ifndef _STLP_NO_LONG_DOUBLE
  37. _STLP_DECLSPEC long double _STLP_CALL abs(const complex<long double>& __z)
  38. {
  39.   return _STLP_HYPOTL(__z._M_re, __z._M_im);
  40. }
  41. #endif
  42.  
  43. // Phase
  44.  
  45. _STLP_DECLSPEC float _STLP_CALL arg(const complex<float>& __z) 
  46. {
  47.   return _STLP_ATAN2F(__z._M_im, __z._M_re);
  48. }
  49.  
  50. _STLP_DECLSPEC double _STLP_CALL arg(const complex<double>& __z) 
  51. {
  52.   return _STLP_ATAN2(__z._M_im, __z._M_re);
  53. }
  54.  
  55. #ifndef _STLP_NO_LONG_DOUBLE
  56. _STLP_DECLSPEC long double _STLP_CALL arg(const complex<long double>& __z) 
  57. {
  58.   return _STLP_ATAN2L(__z._M_im, __z._M_re);
  59. }
  60. #endif
  61.  
  62. // Construct a complex number from polar representation
  63.  
  64. _STLP_DECLSPEC complex<float> _STLP_CALL polar(const float& __rho, const float& __phi) 
  65. {
  66.   return complex<float>(__rho * _STLP_COSF(__phi), __rho * _STLP_SINF(__phi));
  67. }
  68.  
  69. _STLP_DECLSPEC complex<double> _STLP_CALL polar(const double& __rho, const double& __phi) 
  70. {
  71.   return complex<double>(__rho * _STLP_COS(__phi), __rho * _STLP_SIN(__phi));
  72. }
  73.  
  74. #ifndef _STLP_NO_LONG_DOUBLE
  75. _STLP_DECLSPEC complex<long double> _STLP_CALL polar(const long double& __rho, const long double& __phi)
  76. {
  77.   return complex<long double>(__rho * _STLP_COSL(__phi), __rho * _STLP_SINL(__phi));
  78. }
  79. #endif
  80.  
  81. // Division
  82.  
  83. void  _STLP_CALL
  84. complex<float>::_div(const float& __z1_r, const float& __z1_i,
  85.              const float& __z2_r, const float& __z2_i,
  86.              float& __res_r, float& __res_i) {
  87.   float __ar = __z2_r >= 0 ? __z2_r : -__z2_r;
  88.   float __ai = __z2_i >= 0 ? __z2_i : -__z2_i;
  89.  
  90.   if (__ar <= __ai) {
  91.     float __ratio = __z2_r / __z2_i;
  92.     float __denom = __z2_i * (1 + __ratio * __ratio);
  93.     __res_r = (__z1_r * __ratio + __z1_i) / __denom;
  94.     __res_i = (__z1_i * __ratio - __z1_r) / __denom;
  95.   }
  96.   else {
  97.     float __ratio = __z2_i / __z2_r;
  98.     float __denom = __z2_r * (1 + __ratio * __ratio);
  99.     __res_r = (__z1_r + __z1_i * __ratio) / __denom;
  100.     __res_i = (__z1_i - __z1_r * __ratio) / __denom;
  101.   }
  102. }
  103.  
  104. void  _STLP_CALL
  105. complex<float>::_div(const float& __z1_r,
  106.                      const float& __z2_r, const float& __z2_i,
  107.                      float& __res_r, float& __res_i) {
  108.   float __ar = __z2_r >= 0 ? __z2_r : -__z2_r;
  109.   float __ai = __z2_i >= 0 ? __z2_i : -__z2_i;
  110.  
  111.   if (__ar <= __ai) {
  112.     float __ratio = __z2_r / __z2_i;
  113.     float __denom = __z2_i * (1 + __ratio * __ratio);
  114.     __res_r = (__z1_r * __ratio) / __denom;
  115.     __res_i = - __z1_r / __denom;
  116.   }
  117.   else {
  118.     float __ratio = __z2_i / __z2_r;
  119.     float __denom = __z2_r * (1 + __ratio * __ratio);
  120.     __res_r = __z1_r / __denom;
  121.     __res_i = - (__z1_r * __ratio) / __denom;
  122.   }
  123. }
  124.  
  125.  
  126. void  _STLP_CALL
  127. complex<double>::_div(const double& __z1_r, const double& __z1_i,
  128.                       const double& __z2_r, const double& __z2_i,
  129.                       double& __res_r, double& __res_i) {
  130.   double __ar = __z2_r >= 0 ? __z2_r : -__z2_r;
  131.   double __ai = __z2_i >= 0 ? __z2_i : -__z2_i;
  132.  
  133.   if (__ar <= __ai) {
  134.     double __ratio = __z2_r / __z2_i;
  135.     double __denom = __z2_i * (1 + __ratio * __ratio);
  136.     __res_r = (__z1_r * __ratio + __z1_i) / __denom;
  137.     __res_i = (__z1_i * __ratio - __z1_r) / __denom;
  138.   }
  139.   else {
  140.     double __ratio = __z2_i / __z2_r;
  141.     double __denom = __z2_r * (1 + __ratio * __ratio);
  142.     __res_r = (__z1_r + __z1_i * __ratio) / __denom;
  143.     __res_i = (__z1_i - __z1_r * __ratio) / __denom;
  144.   }
  145. }
  146.  
  147. void _STLP_CALL
  148. complex<double>::_div(const double& __z1_r,
  149.                       const double& __z2_r, const double& __z2_i,
  150.                       double& __res_r, double& __res_i) {
  151.   double __ar = __z2_r >= 0 ? __z2_r : -__z2_r;
  152.   double __ai = __z2_i >= 0 ? __z2_i : -__z2_i;
  153.  
  154.   if (__ar <= __ai) {
  155.     double __ratio = __z2_r / __z2_i;
  156.     double __denom = __z2_i * (1 + __ratio * __ratio);
  157.     __res_r = (__z1_r * __ratio) / __denom;
  158.     __res_i = - __z1_r / __denom;
  159.   }
  160.   else {
  161.     double __ratio = __z2_i / __z2_r;
  162.     double __denom = __z2_r * (1 + __ratio * __ratio);
  163.     __res_r = __z1_r / __denom;
  164.     __res_i = - (__z1_r * __ratio) / __denom;
  165.   }
  166. }
  167.  
  168. #ifndef _STLP_NO_LONG_DOUBLE
  169. void  _STLP_CALL
  170. complex<long double>::_div(const long double& __z1_r, const long double& __z1_i,
  171.                            const long double& __z2_r, const long double& __z2_i,
  172.                            long double& __res_r, long double& __res_i) {
  173.   long double __ar = __z2_r >= 0 ? __z2_r : -__z2_r;
  174.   long double __ai = __z2_i >= 0 ? __z2_i : -__z2_i;
  175.  
  176.   if (__ar <= __ai) {
  177.     long double __ratio = __z2_r / __z2_i;
  178.     long double __denom = __z2_i * (1 + __ratio * __ratio);
  179.     __res_r = (__z1_r * __ratio + __z1_i) / __denom;
  180.     __res_i = (__z1_i * __ratio - __z1_r) / __denom;
  181.   }
  182.   else {
  183.     long double __ratio = __z2_i / __z2_r;
  184.     long double __denom = __z2_r * (1 + __ratio * __ratio);
  185.     __res_r = (__z1_r + __z1_i * __ratio) / __denom;
  186.     __res_i = (__z1_i - __z1_r * __ratio) / __denom;
  187.   }
  188. }
  189.  
  190.  
  191. void _STLP_CALL
  192. complex<long double>::_div(const long double& __z1_r,
  193.                            const long double& __z2_r, const long double& __z2_i,
  194.                            long double& __res_r, long double& __res_i) {
  195.   long double __ar = __z2_r >= 0 ? __z2_r : -__z2_r;
  196.   long double __ai = __z2_i >= 0 ? __z2_i : -__z2_i;
  197.  
  198.   if (__ar <= __ai) {
  199.     long double __ratio = __z2_r / __z2_i;
  200.     long double __denom = __z2_i * (1 + __ratio * __ratio);
  201.     __res_r = (__z1_r * __ratio) / __denom;
  202.     __res_i = - __z1_r / __denom;
  203.   }
  204.   else {
  205.     long double __ratio = __z2_i / __z2_r;
  206.     long double __denom = __z2_r * (1 + __ratio * __ratio);
  207.     __res_r = __z1_r / __denom;
  208.     __res_i = - (__z1_r * __ratio) / __denom;
  209.   }
  210. }
  211. #endif
  212.  
  213. //----------------------------------------------------------------------
  214. // Square root
  215.  
  216.  
  217. complex<float> _STLP_CALL
  218. sqrt(const complex<float>& z) {
  219.   float re = z._M_re;
  220.   float im = z._M_im;
  221.   float mag = _STLP_HYPOTF(re, im);
  222.   complex<float> result;
  223.  
  224.   if (mag == 0.) {
  225.     result._M_re = result._M_im = 0.f;
  226.   } else if (re > 0.f) {
  227.     result._M_re = _STLP_SQRTF(0.5f * (mag + re));
  228.     result._M_im = im/result._M_re/2.f;
  229.   } else {
  230.     result._M_im = _STLP_SQRTF(0.5f * (mag - re));
  231.     if (im < 0.f)
  232.       result._M_im = - result._M_im;
  233.     result._M_re = im/result._M_im/2.f;
  234.   }
  235.   return result;
  236. }
  237.  
  238.  
  239. complex<double>  _STLP_CALL
  240. sqrt(const complex<double>& z) {
  241.   double re = z._M_re;
  242.   double im = z._M_im;
  243.   double mag = _STLP_HYPOT(re, im);
  244.   complex<double> result;
  245.  
  246.   if (mag == 0.) {
  247.     result._M_re = result._M_im = 0.;
  248.   } else if (re > 0.) {
  249.     result._M_re = _STLP_SQRT(0.5 * (mag + re));
  250.     result._M_im = im/result._M_re/2;
  251.   } else {
  252.     result._M_im = _STLP_SQRT(0.5 * (mag - re));
  253.     if (im < 0.)
  254.       result._M_im = - result._M_im;
  255.     result._M_re = im/result._M_im/2;
  256.   }
  257.   return result;
  258. }
  259.  
  260. #ifndef _STLP_NO_LONG_DOUBLE
  261. complex<long double> _STLP_CALL
  262. sqrt(const complex<long double>& z) {
  263.   long double re = z._M_re;
  264.   long double im = z._M_im;
  265.   long double mag = _STLP_HYPOTL(re, im);
  266.   complex<long double> result;
  267.  
  268.   if (mag == 0.L) {
  269.     result._M_re = result._M_im = 0.L;
  270.   } else if (re > 0.L) {
  271.     result._M_re = _STLP_SQRTL(0.5L * (mag + re));
  272.     result._M_im = (im/result._M_re) * .5L;
  273.   } else {
  274.     result._M_im = _STLP_SQRTL(0.5L * (mag - re));
  275.     if (im < 0.L)
  276.       result._M_im = - result._M_im;
  277.     result._M_re = (im/result._M_im) * .5L;
  278.   }
  279.   return result;
  280. }
  281. #endif
  282.  
  283. _STLP_END_NAMESPACE
  284.  
  285.