- /* sdoddmul.cpp by K.Tsuru */
- // function ID = 324 DRADIX, BRADIX
- #ifndef SN_H
- #include "sn.h"
- #endif
- /************************************************************
- SDouble class and SDecimal classes m*n
- It uses SLong's multiplication arithmetic because real number
- has much subjects for test such as Sqrt(2.0)^2-2.0 is very
- close to zero or not.Then it is useful in debug of SLong class.
- *************************************************************/
- static const char* const func = "SD *";
- SDouble operator*(const SDouble& m, const SDouble& n){
- //different in radix
- if(m.Type() != n.Type()) m.SetError(m.RADIX_ERR, func, 324);
- //check the sign
- //return zero which has the type of m. SDecimal class also calls.
- if( (m.Sign(324) == 0) || (n.Sign(324) == 0) ) return SDZero(m);
- //For a simplicity of algorithm it lets longer number be multiplicand.
- //To avoid recursive call in DDMult().
- if( m.Last() < n.Last() ) return DDMult(n, m);
- return DDMult(m, n);
- }
-
- SDouble DDMult(const SDouble& m, const SDouble& n){
- //For a simplicity of algorithm it lets longer number be multiplicand.
- if( m.Last() < n.Last() ) return DDMult(n, m);
-
- int sgn = m.Sign(324)*n.Sign(324);
- //check the value of exponent(overflow)
- double dExp = (double)m.RdxExp()+(double)n.RdxExp();
- if(fabs(dExp) > DRADIX_EXP_MAX){
- if(dExp > 0) m.SetError(m.OVERFLOW_ERR, func,324);
- else m.SetError(m.UNDERFLOW_ERR, func,324);
- }
-
- //Is it reducible into SDouble*(short integer)*DRADIX^e ?
- uint nf = n.Last() - n.First(); // 0.0001 2345 6000(nf = 2)=123456*10^(-9) Ok
- if(nf <= 2){
- // |m| == 1 or |n| == 1 ?
- int isOne;
- if( (isOne = n.IsOne()) > 0) return m; // n = 1
- if( isOne < 0) return -m; // n = -1
- if( (isOne = m.IsOne()) > 0) return n; // m = 1
- if( isOne < 0) return -n; // m = -1
- // n = L * 10^ne ?
- long L, ne;
- if( (m.Type() == m.REAL) && n.ConvTolongExp(&L, &ne) ){
- //For SDecimal objects MultPow10() cannot be used.
- SDouble r;
- L = labs(L);
- if( L > 1 ) r = DsMult(m, (ulong)L);
- else r = m; // m*1
- if(ne) r.MultPow10(ne); // includes Reform();
- // ne < 0 and in the fixed point mode it is possible r == 0.
- if( r.Sign(324) ) r.SetSign(sgn);
- return r;
- }
- }
-
- int pt_fx = m.IsPointFixed();
- int top; //top position of non-zero in "result" = number of head zeros + 1
-
- //check m*n==0 in fixed point mode
- // NetRdxExp() = rdxExp - aTail +1
- // 0 -3 top = 1-0-(-3)+2=6 or 5
- if(pt_fx){ // 0.0123e1 * 0.345e(-3) = 0.0123e1 * 0.000345 = 0.0000042435e1
- //It gets the positon when there is no carry.
- top = m.FixedExp() - m.NetRdxExp() - n.NetRdxExp()+2;
- if( (top > 0) && (top >= (int)m.CurrentMaxSize()) ) return SDZero(m);
- //The position is outside of effective figures.
- }
- uint sz, s;
- int rexp;
- if(m.Type() == m.BIN_DEC) sz = m.CurrentMaxSize();
- else sz = min(m.Last()+n.Last(), m.CurrentMaxSize()-1);
-
- // rexp = m.rdxExp + n.rdxExp;
- s = sz/2+1; //divide into "s" (upper) and "s-1" figures(lower).
- sz = 2*s;
- // work area : rs = M*N
- SNumber::NumberType tp = ( m.Type() == m.REAL ) ? m.DEC_INT : m.BIN_INT;
- SLong ml(tp, s), nl(tp, s), rL(tp, sz), temp(tp, sz);
- int sf, e; // e : exponent value considering the carry
- uint cms = m.CurrentMaxSize();
- // Ml = Nl = 0 m*n = Mh*Nh*R^(p+q)
- if( (m.Last() < cms/2) && (n.Last() < cms/2) ){
- m.ConvertToSLong(m.First(), m.Last(), ml);
- if(m == n) rL = ml*ml; // m*m ver. 2.18
- // if(&m == &n) rL = ml*ml;
- else {
- n.ConvertToSLong(n.First(), n.Last(), nl);
- rL = ml*nl; // SLong * SLong (less than "sz" figures)
- }
- e = m.rdxExp - (int)m.Last() + n.rdxExp - (int)n.Last() + (int)rL.Head()+1;
- }
- // Nl = 0 m*n = Mh*Nh*R^s + Ml*Nh (figures of n is smaller than m)
- else if( (m.Last() > cms/2) && (n.Last() < cms/2) ){
- n.ConvertToSLong(n.First(), n.Last(), nl); // N
- m.ConvertToSLong(s + 1, sz, ml); // Ml
- rL = ml*nl; // Ml * N
- //It cuts away lower n.Last() figures.
- sf = (int)n.Last();
- rL.MultPowRdx(-sf); //shift to lower
- m.ConvertToSLong(m.First(), s, ml); // Mh
- nl = ml*nl; // Mh*N
- //It adjusts the position of figures by shifting to upper.
- temp = nl.MultPowRdx((int)s - sf); // s - sf >= 0
- rL += temp;
- e = m.rdxExp + n.rdxExp - (int)sz + (int)rL.Head()+1;
- }
- // M*N = Mh*Nh*R^(2s) + (Mh*Nl + Ml*Nh)*R^s
- else {
- SLong mh(tp, s);
- sf = -(int)s;
- m.ConvertToSLong(m.First(), s, mh); // Mh
- n.ConvertToSLong(s + 1, sz, nl); // Nl
- rL = mh*nl; // Mh*Nl
- /*
- The speed-up for square.
- ver.2.18 changed from "if(&m == &n)"
- because overhead is very small.
- */
- bool square = (m == n); // ver. 2.18
- // bool square = (&m == &n);
- if(square){ // x*x m == n Mh*Nl == Ml*Nh
- rL = LsMult(rL, 2);
- } else {
- n.ConvertToSLong(n.First(), s, nl); //Nh
- m.ConvertToSLong(s + 1, sz, ml); //Ml
- temp = ml*nl; //Ml*Nh
- rL += temp;
- }
-
- rL.MultPowRdx(sf);
-
- if(square) temp = mh*mh;
- else temp = mh*nl; // Mh*Nh, nl=Nh
-
- rL += temp;
- e = m.rdxExp + n.rdxExp - (int)sz + (int)rL.Head()+1;
- }
-
- if(m.IsPointFixed()){
- //It is possible in the fixed point mode.
- if( !rL.Sign(324) ) return SDZero(m);
- rexp = m.FixedExp();
- top = 1 + m.FixedExp() - e;
- if(top <= 0 ){ //It cannot be represented with fixed exponent value.
- if(m.Type() == m.REAL){
- rexp += 1 - top; top = 1;
- } else if(top < 0) m.SetError(m.OVERFLOW_ERR,"DDMult", 324);
- } else if( (uint)top >= m.CurrentMaxSize() ) return SDZero(m);
- } else {
- top = 1; rexp = e;
- }
- //free memory
- ml.SizeZero(); nl.SizeZero(); temp.SizeZero();
- //SLong ==> SDouble
- SDouble result(m.Type(), 0);
-
- register uint j;
- const fType* rLf = rL.ReadFigures();
- uint rsz;
- s = rL.Head() + top;
- sz = m.CurrentMaxSize();
- rsz = min(s+1, sz);
- result.valloc(rsz, -1);
- sz = min( s+1, result.Size()); //really allocated size
-
- fType* rv = result.figure.Elements();
- //set each figure
- #ifndef NDEBUG
- result.figure(sz-1);
- assert( (s-top) < rL.Size() );
- #endif
- //s - j = rL.Head() ...0
- for(j = top; j < sz ; j++) rv[j] = rLf[s-j];
- if(top) result.figure.clear(0, top-1);
- result.figure.clear(sz);
- //get positions
- j = top;
- while(!rv[j] && (j < sz) ) j++;
-
- if(j == sz){
- result.SetZero(); return result;
- }
- result.aTail = j;
-
- j = sz-1;
- while(!rv[j]) j--;
- result.aHead = j;
- // Set exponent and sign.
- // Checking |rexp| < DRADIX_EXP_MAX
- result.SetRdxExp(rexp);
- result.SetSign(sgn);
- result.Reform(324);
- return result;
- }
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).