PageRenderTime 42ms CodeModel.GetById 19ms RepoModel.GetById 0ms app.codeStats 0ms

/cln-1.3.2/src/float/lfloat/elem/cl_LF_I_mul.cc

#
C++ | 91 lines | 62 code | 8 blank | 21 comment | 18 complexity | da163b1ff652e6e23f65b7a1fe28b48a MD5 | raw file
Possible License(s): GPL-2.0
  1. // cl_LF_I_mul().
  2. // General includes.
  3. #include "base/cl_sysdep.h"
  4. // Specification.
  5. #include "float/lfloat/cl_LF.h"
  6. // Implementation.
  7. #include "float/lfloat/cl_LF_impl.h"
  8. #include "cln/integer.h"
  9. #include "integer/cl_I.h"
  10. #include "base/digitseq/cl_DS.h"
  11. #include "float/cl_F.h"
  12. namespace cln {
  13. const cl_R cl_LF_I_mul (const cl_LF& x, const cl_I& y)
  14. {
  15. // Method:
  16. // If y=0, return 0.
  17. // If x=0.0, return x.
  18. // If y is longer than x, convert y to a float and multiply.
  19. // Else multiply the mantissa of x with the absolute value of y, then round.
  20. if (eq(y,0)) { return 0; }
  21. if (TheLfloat(x)->expo == 0) { return x; }
  22. var cl_signean sign = -(cl_signean)minusp(y); // Vorzeichen von y
  23. var cl_I abs_y = (sign==0 ? y : -y);
  24. var uintC y_exp = integer_length(abs_y);
  25. var uintC len = TheLfloat(x)->len;
  26. #ifndef CL_LF_PEDANTIC
  27. if (ceiling(y_exp,intDsize) > len)
  28. return x * cl_I_to_LF(y,len);
  29. #endif
  30. // x länger als y, direkt multiplizieren.
  31. CL_ALLOCA_STACK;
  32. var const uintD* y_MSDptr;
  33. var uintC y_len;
  34. var const uintD* y_LSDptr;
  35. I_to_NDS_nocopy(abs_y, y_MSDptr=,y_len=,y_LSDptr=,false,); // NDS zu y bilden, y_len>0
  36. if (mspref(y_MSDptr,0)==0) y_len--; // NUDS zu y bilden, y_len>0
  37. // Multiplizieren.
  38. var uintD* prodMSDptr;
  39. var uintC prodlen;
  40. UDS_UDS_mul_UDS(len,arrayLSDptr(TheLfloat(x)->data,len),
  41. y_len,y_LSDptr,
  42. prodMSDptr=,prodlen=,);
  43. // x fing mit 0 Nullbits an, y mit maximal intDsize-1 Nullbits,
  44. // daher fängt das Produkt mit maximal intDsize Nullbits an.
  45. var uintL shiftcount;
  46. if (mspref(prodMSDptr,0)==0) {
  47. shiftcount = intDsize;
  48. msshrink(prodMSDptr); prodlen--;
  49. } else {
  50. integerlengthD(mspref(prodMSDptr,0), shiftcount = intDsize -);
  51. if (shiftcount > 0)
  52. shiftleft_loop_lsp(prodMSDptr mspop (len+1),len+1,shiftcount,0);
  53. }
  54. // Produkt ist nun normalisiert: höchstes Bit =1.
  55. // exponent := exponent(x) + intDsize*y_len - shiftcount
  56. var uintE uexp = TheLfloat(x)->expo;
  57. var uintE iexp = intDsize*y_len - shiftcount; // >= 0 !
  58. uexp = uexp + iexp;
  59. if ((uexp < iexp) || (uexp > LF_exp_high))
  60. throw floating_point_overflow_exception();
  61. // Runden:
  62. var uintD* midptr = prodMSDptr mspop len;
  63. var uintC restlen = prodlen - len;
  64. if ( (restlen==0)
  65. || ((sintD)mspref(midptr,0) >= 0) // nächstes Bit =0 -> abrunden
  66. || ( ((mspref(midptr,0) & ((uintD)bit(intDsize-1)-1)) ==0) // Bit =1, weitere Bits >0 -> aufrunden
  67. && !test_loop_msp(midptr mspop 1,restlen-1)
  68. // round-to-even
  69. && ((lspref(midptr,0) & bit(0)) ==0)
  70. ) )
  71. // abrunden
  72. {}
  73. else
  74. // aufrunden
  75. { if ( inc_loop_lsp(midptr,len) )
  76. // Übertrag durchs Aufrunden
  77. { mspref(prodMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
  78. if (++uexp == LF_exp_high+1) { throw floating_point_overflow_exception(); }
  79. } }
  80. return encode_LFu(TheLfloat(x)->sign ^ sign, uexp, prodMSDptr, len);
  81. }
  82. // Bit complexity (N = max(length(x),length(y))): O(M(N)).
  83. } // namespace cln