PageRenderTime 49ms CodeModel.GetById 11ms RepoModel.GetById 0ms app.codeStats 0ms

/cln-1.3.2/src/complex/algebraic/cl_LF_hypot.cc

#
C++ | 77 lines | 44 code | 10 blank | 23 comment | 9 complexity | 160c2a32322aea089c4ffcafc442bc29 MD5 | raw file
Possible License(s): GPL-2.0
  1. // cl_hypot().
  2. // General includes.
  3. #include "base/cl_sysdep.h"
  4. // Specification.
  5. #include "complex/cl_C.h"
  6. // Implementation.
  7. #include "cln/lfloat.h"
  8. #include "float/lfloat/cl_LF.h"
  9. #include "float/lfloat/cl_LF_impl.h"
  10. /* For inline version of minusp */
  11. #include "base/cl_inline.h"
  12. #include "float/lfloat/elem/cl_LF_minusp.cc"
  13. namespace cln {
  14. ALL_cl_LF_OPERATIONS_SAME_PRECISION()
  15. const cl_LF cl_hypot (const cl_LF& a, const cl_LF& b)
  16. {
  17. // Zuerst a und b auf gleiche Länge bringen: den längeren runden.
  18. // a=0.0 -> liefere abs(b).
  19. // b=0.0 -> liefere abs(a).
  20. // e:=max(exponent(a),exponent(b)).
  21. // a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
  22. // oder beim Quadrieren a'*a': 2*(e-exponent(a))>exp_mid-exp_low-1
  23. // d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
  24. // b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
  25. // oder beim Quadrieren b'*b': 2*(e-exponent(b))>exp_mid-exp_low-1
  26. // d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
  27. // c':=a'*a'+b'*b', c':=sqrt(c'), liefere 2^e*c'.
  28. { Mutable(cl_LF,a);
  29. Mutable(cl_LF,b);
  30. {
  31. var uintC a_len = TheLfloat(a)->len;
  32. var uintC b_len = TheLfloat(b)->len;
  33. if (!(a_len == b_len)) {
  34. if (a_len < b_len)
  35. b = shorten(b,a_len);
  36. else
  37. a = shorten(a,b_len);
  38. }
  39. }
  40. var sintE a_exp;
  41. var sintE b_exp;
  42. {
  43. // Exponenten von a holen:
  44. var uintE uexp = TheLfloat(a)->expo;
  45. if (uexp == 0)
  46. // a=0.0 -> liefere (abs b) :
  47. return (minusp_inline(b) ? -b : b);
  48. a_exp = (sintE)(uexp - LF_exp_mid);
  49. }
  50. {
  51. // Exponenten von b holen:
  52. var uintE uexp = TheLfloat(b)->expo;
  53. if (uexp == 0)
  54. // b=0.0 -> liefere (abs a) :
  55. return (minusp_inline(a) ? -a : a);
  56. b_exp = (sintE)(uexp - LF_exp_mid);
  57. }
  58. // Nun a_exp = float_exponent(a), b_exp = float_exponent(b).
  59. var sintE e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten
  60. // a und b durch 2^e dividieren:
  61. var cl_LF na = ((b_exp > a_exp) && ((uintE)(b_exp-a_exp) > (uintE)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(a)->len) : scale_float(a,-e));
  62. var cl_LF nb = ((a_exp > b_exp) && ((uintE)(a_exp-b_exp) > (uintE)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(b)->len) : scale_float(b,-e));
  63. // c' := a'*a'+b'*b' berechnen:
  64. var cl_LF nc = square(na) + square(nb);
  65. return scale_float(sqrt(nc),e); // c' := sqrt(c'), 2^e*c'
  66. }}
  67. } // namespace cln