1. /* sdoddmul.cpp by K.Tsuru */
  2. // function ID = 324 DRADIX, BRADIX
  3. #ifndef SN_H
  4. #include "sn.h"
  5. #endif
  6. /************************************************************
  7. SDouble class and SDecimal classes m*n
  8. It uses SLong's multiplication arithmetic because real number
  9. has much subjects for test such as Sqrt(2.0)^2-2.0 is very
  10. close to zero or not.Then it is useful in debug of SLong class.
  11. *************************************************************/
  12. static const char* const func = "SD *";
  13. SDouble operator*(const SDouble& m, const SDouble& n){
  14. //different in radix
  15. if(m.Type() != n.Type()) m.SetError(m.RADIX_ERR, func, 324);
  16. //check the sign
  17. //return zero which has the type of m. SDecimal class also calls.
  18. if( (m.Sign(324) == 0) || (n.Sign(324) == 0) ) return SDZero(m);
  19. //For a simplicity of algorithm it lets longer number be multiplicand.
  20. //To avoid recursive call in DDMult().
  21. if( m.Last() < n.Last() ) return DDMult(n, m);
  22. return DDMult(m, n);
  23. }
  24. SDouble DDMult(const SDouble& m, const SDouble& n){
  25. //For a simplicity of algorithm it lets longer number be multiplicand.
  26. if( m.Last() < n.Last() ) return DDMult(n, m);
  27. int sgn = m.Sign(324)*n.Sign(324);
  28. //check the value of exponent(overflow)
  29. double dExp = (double)m.RdxExp()+(double)n.RdxExp();
  30. if(fabs(dExp) > DRADIX_EXP_MAX){
  31. if(dExp > 0) m.SetError(m.OVERFLOW_ERR, func,324);
  32. else m.SetError(m.UNDERFLOW_ERR, func,324);
  33. }
  34. //Is it reducible into SDouble*(short integer)*DRADIX^e ?
  35. uint nf = n.Last() - n.First(); // 0.0001 2345 6000(nf = 2)=123456*10^(-9) Ok
  36. if(nf <= 2){
  37. // |m| == 1 or |n| == 1 ?
  38. int isOne;
  39. if( (isOne = n.IsOne()) > 0) return m; // n = 1
  40. if( isOne < 0) return -m; // n = -1
  41. if( (isOne = m.IsOne()) > 0) return n; // m = 1
  42. if( isOne < 0) return -n; // m = -1
  43. // n = L * 10^ne ?
  44. long L, ne;
  45. if( (m.Type() == m.REAL) && n.ConvTolongExp(&L, &ne) ){
  46. //For SDecimal objects MultPow10() cannot be used.
  47. SDouble r;
  48. L = labs(L);
  49. if( L > 1 ) r = DsMult(m, (ulong)L);
  50. else r = m; // m*1
  51. if(ne) r.MultPow10(ne); // includes Reform();
  52. // ne < 0 and in the fixed point mode it is possible r == 0.
  53. if( r.Sign(324) ) r.SetSign(sgn);
  54. return r;
  55. }
  56. }
  57. int pt_fx = m.IsPointFixed();
  58. int top; //top position of non-zero in "result" = number of head zeros + 1
  59. //check m*n==0 in fixed point mode
  60. // NetRdxExp() = rdxExp - aTail +1
  61. // 0 -3 top = 1-0-(-3)+2=6 or 5
  62. if(pt_fx){ // 0.0123e1 * 0.345e(-3) = 0.0123e1 * 0.000345 = 0.0000042435e1
  63. //It gets the positon when there is no carry.
  64. top = m.FixedExp() - m.NetRdxExp() - n.NetRdxExp()+2;
  65. if( (top > 0) && (top >= (int)m.CurrentMaxSize()) ) return SDZero(m);
  66. //The position is outside of effective figures.
  67. }
  68. uint sz, s;
  69. int rexp;
  70. if(m.Type() == m.BIN_DEC) sz = m.CurrentMaxSize();
  71. else sz = min(m.Last()+n.Last(), m.CurrentMaxSize()-1);
  72. // rexp = m.rdxExp + n.rdxExp;
  73. s = sz/2+1; //divide into "s" (upper) and "s-1" figures(lower).
  74. sz = 2*s;
  75. // work area : rs = M*N
  76. SNumber::NumberType tp = ( m.Type() == m.REAL ) ? m.DEC_INT : m.BIN_INT;
  77. SLong ml(tp, s), nl(tp, s), rL(tp, sz), temp(tp, sz);
  78. int sf, e; // e : exponent value considering the carry
  79. uint cms = m.CurrentMaxSize();
  80. // Ml = Nl = 0 m*n = Mh*Nh*R^(p+q)
  81. if( (m.Last() < cms/2) && (n.Last() < cms/2) ){
  82. m.ConvertToSLong(m.First(), m.Last(), ml);
  83. if(m == n) rL = ml*ml; // m*m ver. 2.18
  84. // if(&m == &n) rL = ml*ml;
  85. else {
  86. n.ConvertToSLong(n.First(), n.Last(), nl);
  87. rL = ml*nl; // SLong * SLong (less than "sz" figures)
  88. }
  89. e = m.rdxExp - (int)m.Last() + n.rdxExp - (int)n.Last() + (int)rL.Head()+1;
  90. }
  91. // Nl = 0 m*n = Mh*Nh*R^s + Ml*Nh (figures of n is smaller than m)
  92. else if( (m.Last() > cms/2) && (n.Last() < cms/2) ){
  93. n.ConvertToSLong(n.First(), n.Last(), nl); // N
  94. m.ConvertToSLong(s + 1, sz, ml); // Ml
  95. rL = ml*nl; // Ml * N
  96. //It cuts away lower n.Last() figures.
  97. sf = (int)n.Last();
  98. rL.MultPowRdx(-sf); //shift to lower
  99. m.ConvertToSLong(m.First(), s, ml); // Mh
  100. nl = ml*nl; // Mh*N
  101. //It adjusts the position of figures by shifting to upper.
  102. temp = nl.MultPowRdx((int)s - sf); // s - sf >= 0
  103. rL += temp;
  104. e = m.rdxExp + n.rdxExp - (int)sz + (int)rL.Head()+1;
  105. }
  106. // M*N = Mh*Nh*R^(2s) + (Mh*Nl + Ml*Nh)*R^s
  107. else {
  108. SLong mh(tp, s);
  109. sf = -(int)s;
  110. m.ConvertToSLong(m.First(), s, mh); // Mh
  111. n.ConvertToSLong(s + 1, sz, nl); // Nl
  112. rL = mh*nl; // Mh*Nl
  113. /*
  114. The speed-up for square.
  115. ver.2.18 changed from "if(&m == &n)"
  116. because overhead is very small.
  117. */
  118. bool square = (m == n); // ver. 2.18
  119. // bool square = (&m == &n);
  120. if(square){ // x*x m == n Mh*Nl == Ml*Nh
  121. rL = LsMult(rL, 2);
  122. } else {
  123. n.ConvertToSLong(n.First(), s, nl); //Nh
  124. m.ConvertToSLong(s + 1, sz, ml); //Ml
  125. temp = ml*nl; //Ml*Nh
  126. rL += temp;
  127. }
  128. rL.MultPowRdx(sf);
  129. if(square) temp = mh*mh;
  130. else temp = mh*nl; // Mh*Nh, nl=Nh
  131. rL += temp;
  132. e = m.rdxExp + n.rdxExp - (int)sz + (int)rL.Head()+1;
  133. }
  134. if(m.IsPointFixed()){
  135. //It is possible in the fixed point mode.
  136. if( !rL.Sign(324) ) return SDZero(m);
  137. rexp = m.FixedExp();
  138. top = 1 + m.FixedExp() - e;
  139. if(top <= 0 ){ //It cannot be represented with fixed exponent value.
  140. if(m.Type() == m.REAL){
  141. rexp += 1 - top; top = 1;
  142. } else if(top < 0) m.SetError(m.OVERFLOW_ERR,"DDMult", 324);
  143. } else if( (uint)top >= m.CurrentMaxSize() ) return SDZero(m);
  144. } else {
  145. top = 1; rexp = e;
  146. }
  147. //free memory
  148. ml.SizeZero(); nl.SizeZero(); temp.SizeZero();
  149. //SLong ==> SDouble
  150. SDouble result(m.Type(), 0);
  151. register uint j;
  152. const fType* rLf = rL.ReadFigures();
  153. uint rsz;
  154. s = rL.Head() + top;
  155. sz = m.CurrentMaxSize();
  156. rsz = min(s+1, sz);
  157. result.valloc(rsz, -1);
  158. sz = min( s+1, result.Size()); //really allocated size
  159. fType* rv = result.figure.Elements();
  160. //set each figure
  161. #ifndef NDEBUG
  162. result.figure(sz-1);
  163. assert( (s-top) < rL.Size() );
  164. #endif
  165. //s - j = rL.Head() ...0
  166. for(j = top; j < sz ; j++) rv[j] = rLf[s-j];
  167. if(top) result.figure.clear(0, top-1);
  168. result.figure.clear(sz);
  169. //get positions
  170. j = top;
  171. while(!rv[j] && (j < sz) ) j++;
  172. if(j == sz){
  173. result.SetZero(); return result;
  174. }
  175. result.aTail = j;
  176. j = sz-1;
  177. while(!rv[j]) j--;
  178. result.aHead = j;
  179. // Set exponent and sign.
  180. // Checking |rexp| < DRADIX_EXP_MAX
  181. result.SetRdxExp(rexp);
  182. result.SetSign(sgn);
  183. result.Reform(324);
  184. return result;
  185. }

sdoddmul.cpp : last modifiled at 2017/03/17 11:10:47(6,645 bytes)
created at 2017/10/07 10:21:14
The creation time of this html file is 2017/10/07 10:30:03 (Sat Oct 07 10:30:03 2017).