home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / oath.lha / oath / src / bigInteger.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1991-08-29  |  30.7 KB  |  1,157 lines
  1. //***************************************************************************
  2. //             OATH :: Object-oriented Abstract Type Hierarchy
  3. //***************************************************************************
  4. //
  5. //  Copyright (C) 1991  Texas Instruments Incorporated
  6. //  Permission is granted to any individual or institution
  7. //  to use, copy, modify, and distribute this software,
  8. //  provided that this complete copyright and permission notice
  9. //  is maintained, intact, in all copies and supporting documentation.
  10. //
  11. //  Texas Instruments Incorporated provides this software "as is"
  12. //  without express or implied warranty.
  13. //
  14. //***************************************************************************
  15. //  bigInteger (bigIntegerA, bigIntegerG)
  16. //
  17. //  History:
  18. //    08/91  Brian M Kennedy  signMagnitudeP; functionality refined
  19. //    07/91  Brian M Kennedy  import, export, typeRegister
  20. //    06/91  Brian M Kennedy  New macros & format; remove printDiagnostic
  21. //    02/91  Brian M Kennedy  Original
  22. //
  23. //***************************************************************************
  24.  
  25. #include "copyright.h"
  26.  
  27. #include <oath/bigInteger.h>
  28.  
  29. #include <ctype.h>
  30.  
  31. #include <limits.h>
  32.  
  33. #include <iostream.h>
  34.  
  35. /////////////////////////////////////////////////////////////////////////////
  36. // unsigned long manipulations
  37.  
  38. #define FULL 32  /* (CHAR_BIT * sizeof(unsigned long)) */
  39.  
  40. #define HALF 16  /* (FULL >> 1) */
  41.  
  42. #define LOWER (~0UL >> HALF)
  43.  
  44. #define LO(Arg) ((Arg) & LOWER)
  45.  
  46. #define HI(Arg) ((Arg) >> HALF)
  47.  
  48. #define UP(Arg) ((Arg) << HALF)
  49.  
  50.     static unsigned int        // returns bit number 31=MSB, 0=LSB of MS1
  51. log2 (unsigned long UL)
  52.    {assumed(UL, "log2 not valid on zero");
  53.     unsigned int UI = 0;
  54.     unsigned long UM = HI(UL);
  55.     if(UM)
  56.        {UL = UM;
  57.         UI += HALF;
  58.        }
  59.     while(UL)
  60.        {UL >>= 1;
  61.         ++UI;
  62.        }
  63.     return UI - 1;
  64.    }
  65.  
  66.  
  67. /////////////////////////////////////////////////////////////////////////////
  68. // common call sequences
  69.  
  70. #define SM_CALL(Arg, Fn)                                \
  71.     if(Arg->isImplementedAs(Type))                        \
  72.         Fn(((bigIntegerG*)Arg)->SAM);                        \
  73.     else if(Arg->isType(integerG::Type))                    \
  74.         Fn(((integerG*)Arg)->smMantissa());                    \
  75.     else if(Arg->isType(realG::Type))                        \
  76.        {integerA Temp = ((realG*)Arg)->makeTruncate(1);                \
  77.         Fn(Temp.guts()->smMantissa());                        \
  78.        }                                    \
  79.     else if(Arg->hasZeroImag())                            \
  80.        {integerA Temp = Arg->makeReal(1).makeTruncate(1);            \
  81.         Fn(Temp.guts()->smMantissa());                        \
  82.        }                                    \
  83.     else                                    \
  84.        {ensure(FALSE, "throw offRealAxis::rangeViolation");}            \
  85.  
  86.  
  87. /////////////////////////////////////////////////////////////////////////////
  88. // bigInteger Outlines
  89.  
  90. OUTLINES(bigInteger, integer)
  91.  
  92. // Static Data Members //////////
  93.  
  94. unsigned long  bigIntegerG::ZeroSAMMagnitude; // default init to 0
  95.  
  96. signMagnitudeP bigIntegerG::ZeroSAM (3, &bigIntegerG::ZeroSAMMagnitude);
  97.  
  98. // Constructors //////////
  99.  
  100.     bigIntegerG::
  101. bigIntegerG (const char* S, int IsConst)
  102.    :integerG(IsConst), MaxSize(1), SAM(1, 1, new unsigned long [1])
  103.    {SAM.magnitude(0) = 0;
  104.     while(*S)
  105.        {if(isdigit(*S))
  106.            {smMul10();
  107.         smIncMag(*S - '0');
  108.        }
  109.         else if (*S == '-')
  110.         neg();
  111.     ++S;
  112.        }
  113.    }
  114.  
  115.     bigIntegerG::
  116. bigIntegerG (const bigIntegerG* I, int IsConst)
  117.    :integerG(IsConst),
  118.     MaxSize((I->SAM.size() > 1) ? I->SAM.size() + 4 : 1),
  119.     SAM(I->SAM.SizeAndSign, new unsigned long [MaxSize])
  120.    {for(unsigned int UI = 0; UI < SAM.size(); ++UI)
  121.         SAM.magnitude(UI) = I->SAM.magnitude(UI);
  122.    }
  123.  
  124.     bigIntegerG::
  125. bigIntegerG (const signMagnitudeP& ISAM, int IsConst)
  126.    :integerG(IsConst),
  127.     MaxSize((ISAM.size() > 1) ? ISAM.size() + 4 : 1),
  128.     SAM(ISAM.SizeAndSign, new unsigned long [MaxSize])
  129.    {for(unsigned int UI = 0; UI < SAM.size(); ++UI)
  130.         SAM.magnitude(UI) = ISAM.magnitude(UI);
  131.    }
  132.  
  133.  
  134. // oathCore Operations //////////
  135.  
  136.     void bigIntegerG::
  137. export (exportP& X) const
  138.    {X.writeType(TypeName);
  139.     X.stream() << SAM.SizeAndSign << (isConst() ? ' ' : '\0');
  140.     for(unsigned int UI = 0; UI < SAM.size(); ++UI)
  141.         X.stream() << SAM.magnitude(UI) << ' ';
  142.    }
  143.  
  144.     objA bigIntegerG::
  145. import (importP& M)
  146.    {unsigned int SAS;
  147.     M.stream() >> SAS;
  148.     char MakeConst = M.stream().get();
  149.     unsigned int Size = SAS >> 1;
  150.     unsigned int MaxSize = (Size > 1) ? Size + 4 : 1;
  151.     unsigned long * Mag = new unsigned long [MaxSize];
  152.     for(unsigned int UI = 0; UI < Size; ++UI)
  153.         (M.stream() >> Mag[UI]).get();
  154.     return new bigIntegerG (MaxSize, SAS, Mag, MakeConst);
  155.    }
  156.  
  157.  
  158. // obj Operations //////////
  159.  
  160.     int bigIntegerG::
  161. isEqual (const objG* O) const
  162.    {return O->isType(integerG::Type) && hasEqual((integerG*)O);}
  163.  
  164.  
  165. // complex Operations //////////
  166.  
  167.     void bigIntegerG::
  168. reset ()
  169.    {SAM.size(1);
  170.     SAM.magnitude(0) = 0;
  171.    }
  172.  
  173.     void bigIntegerG::
  174. reset (const complexG* C)
  175.    {SM_CALL(C, smReset)}
  176.  
  177.     int bigIntegerG::
  178. hasEqual (const complexG* C) const
  179.    {if(C->isImplementedAs(Type))
  180.         return smHasEqual(((bigIntegerG*)C)->SAM);
  181.     else if(C->isType(integerG::Type))
  182.         return smHasEqual(((integerG*)C)->smMantissa());
  183.     else
  184.         return C->hasEqual(this);
  185.    }
  186.  
  187.     void bigIntegerG::
  188. increment (const complexG* C)
  189.    {SM_CALL(C, smInc)}
  190.  
  191.     void bigIntegerG::
  192. decrement (const complexG* C)
  193.    {SM_CALL(C, smDec)}
  194.  
  195.     void bigIntegerG::
  196. multiply (const complexG* C)
  197.    {if(C->isImplementedAs(Type))
  198.         smMul(((bigIntegerG*)C)->SAM);
  199.     else if(C->isType(integerG::Type))
  200.         smMul(((integerG*)C)->smMantissa());
  201.     else if(hasZero() || C->hasZero())
  202.         reset();
  203.     else if(C->isType(realG::Type))
  204.        {realA R = (realA&)(((realG*)C)->makeCopy(FALSE));
  205.         R.guts()->multiply(this);
  206.     reset(R.guts());
  207.        }
  208.     else if(C->hasZeroImag())
  209.        {realA R = C->makeReal(1);
  210.         R.guts()->multiply(this);
  211.     reset(R.guts());
  212.        }
  213.     else
  214.        {ensure(FALSE, "throw offRealAxis::rangeViolation");}
  215.    }
  216.  
  217.     void bigIntegerG::
  218. divide (const complexG* C)
  219.    {if(C->isImplementedAs(Type))
  220.         smDiv(((bigIntegerG*)C)->SAM);
  221.     else if(C->isType(integerG::Type))
  222.         smDiv(((integerG*)C)->smMantissa());
  223.     else if(hasZero() || C->hasZero())
  224.         reset();
  225.     else if(C->isType(realG::Type))
  226.        {realA R = (realA&)(((realG*)C)->makeCopy(FALSE));
  227.         R.inv();
  228.         R.guts()->multiply(this);
  229.     reset(R.guts());
  230.        }
  231.     else if(C->hasZeroImag())
  232.        {realA R = C->makeReal(1);
  233.         R.inv();
  234.         R.guts()->multiply(this);
  235.     reset(R.guts());
  236.        }
  237.     else
  238.        {ensure(FALSE, "throw offRealAxis::rangeViolation");}
  239.    }
  240.  
  241.     void bigIntegerG::
  242. shiftLeft (const integerG* I)
  243.    {smShiftLeft(I->smMantissa());}
  244.  
  245.     void bigIntegerG::
  246. shiftRight (const integerG* I)
  247.    {smShiftRight(I->smMantissa());}
  248.  
  249.     void bigIntegerG::
  250. sqr ()
  251.    {signMagnitudeP Copy (SAM);
  252.     Copy.magnitude() = new unsigned long [Copy.size()];
  253.     for(unsigned int UI = 0; UI < Copy.size(); ++UI)
  254.         Copy.magnitude(UI) = SAM.magnitude(UI);
  255.     smMul(Copy);
  256.     delete [] Copy.magnitude();
  257.    }
  258.  
  259.     void bigIntegerG::        // This is just Newton's method, with a 
  260. sqrt ()                // start value x >>= (log2(x) / 2)
  261.    {if(!SAM.zero())
  262.        {if(!SAM.positive())
  263.             ensure(FALSE, "throw notOnRealAxis::rangeViolation");
  264.     else
  265.        {// Temporaries
  266.         bigIntegerA  GuessA = (bigIntegerA&)makeCopy(FALSE);
  267.         bigIntegerG* Guess  = GuessA.guts();
  268.         bigIntegerA  QuotA  = (bigIntegerA&)makeCopy(FALSE);
  269.         bigIntegerG* Quot   = QuotA.guts();
  270.         // Initial guess, always high but fairly close
  271.         unsigned int Count = SAM.size() - 1;
  272.         Count *= FULL;
  273.         Count += log2(SAM.magnitude(SAM.size()-1));
  274.         Count >>= 1;
  275.         Guess->smShiftRight(Count);
  276.         // Iterative Convergence
  277.         Quot->divide(Guess);
  278.         while(Quot->hasLess(Guess))
  279.            {Guess->increment(Quot);
  280.             Guess->smShiftRight(1);
  281.         Quot->reset(this);
  282.             Quot->divide(Guess);
  283.            }
  284.         reset(Guess);
  285.        }
  286.        }
  287.    }
  288.  
  289.     void bigIntegerG::
  290. pow (const integerG* I)
  291.    {if(I->hasZero())
  292.         smReset(1);
  293.     else if(hasZero())
  294.         reset();
  295.     else
  296.        {const signMagnitudeP& SMI = I->smMantissa();
  297.         if(SMI.size() > 1 || !SMI.positive())
  298.         ensure(FALSE, "throw rangeViolation");
  299.     else
  300.        {bigIntegerA I = (bigIntegerA&)makeCopy(TRUE);
  301.         for(unsigned int Count = 0; Count < SMI.magnitude(0); ++Count)
  302.             smMul(I.guts()->SAM);
  303.        }
  304.        }
  305.    }
  306.  
  307.     void bigIntegerG::
  308. print (ostream& F, char Sep) const
  309.    {int Size = SAM.size() * 10 + 12;
  310.     if(Sep) Size += Size >> 1;
  311.     char * Buffer = new char [Size];
  312.     F << string(Buffer, Size, Sep);
  313.     delete [] Buffer;
  314.    }
  315.  
  316. // real Operations //////////
  317.  
  318.     int bigIntegerG::
  319. hasLess (const realG* R) const
  320.    {if(R->isImplementedAs(Type))
  321.         return smHasLess(((bigIntegerG*)R)->SAM);
  322.     else if(R->isType(integerG::Type))
  323.         return smHasLess(((integerG*)R)->smMantissa());
  324.     else
  325.         return R->hasMore(this);
  326.    }
  327.  
  328.     int bigIntegerG::
  329. hasMore (const realG* R) const
  330.    {if(R->isImplementedAs(Type))
  331.         return smHasMore(((bigIntegerG*)R)->SAM);
  332.     else if(R->isType(integerG::Type))
  333.         return smHasMore(((integerG*)R)->smMantissa());
  334.     else
  335.         return R->hasLess(this);
  336.    }
  337.  
  338.     int bigIntegerG::
  339. hasDouble () const
  340.    {ensure(FALSE, "Sorry, not implemented!");
  341.     return FALSE;
  342.    }
  343.  
  344.     double bigIntegerG::
  345. makeDouble () const
  346.    {ensure(FALSE, "Sorry, not implemented!");
  347.     return 0.0;
  348.    }
  349.  
  350.     void bigIntegerG::
  351. min (const realG* R)
  352.    {if(hasMore(R))
  353.         reset(R);
  354.    }
  355.  
  356.     void bigIntegerG::
  357. max (const realG* R)
  358.    {if(hasLess(R))
  359.         reset(R);
  360.    }
  361.  
  362.  
  363. // rational Operations //////////
  364.  
  365.  
  366. // integer Operations //////////
  367.  
  368.     void bigIntegerG::
  369. rem (const integerG* I)
  370.    {signMagnitudeP Quotient;
  371.     unsigned int QMaxSize = 0;
  372.     smRemDivMag(I->smMantissa(), Quotient, QMaxSize);
  373.     delete [] Quotient.magnitude();
  374.    }
  375.  
  376.     void bigIntegerG::
  377. remDivX (const integerG* Divisor, integerG* Quot)
  378.    {bigIntegerA Quotient = bigIntegerA::make();
  379.     smRemDivMag(Divisor->smMantissa(),
  380.             Quotient.guts()->SAM, Quotient.guts()->MaxSize);
  381.     Quot->reset(Quotient.guts());
  382.    }
  383.  
  384. // Common Factors
  385.  
  386.     // This is based on Knuth, Art of Computer Programming, Vol 2
  387.     // Algorithm B (binary gcd algorithm) discovered by J. Stein.
  388.     // This is an alternative to Euclidean algorithms that uses
  389.     // only subtraction and shifting (no division).
  390.  
  391.     void bigIntegerG::
  392. gcd (const integerG* I)
  393.    {bigIntegerA One = bigIntegerA::make(1);
  394.     if(hasZero() && I->hasZero())
  395.         reset();
  396.     else
  397.        {abs();
  398.         integerA V = ((integerA&)I->makeCopy(FALSE)).abs();
  399.     unsigned int K = 0;
  400.         while(hasEven() && V.hasEven())
  401.        {++K;
  402.         smShiftRight(1);
  403.         V >>= One;
  404.        }
  405.     integerA T;
  406.     if(!hasEven())
  407.             T = V.makeCopy().neg();
  408.     else
  409.             T = (integerA&)makeCopy(FALSE);
  410.     while(!T.hasZero())
  411.        {while(T.hasEven())
  412.                 T >>= One;
  413.         if(T.hasPositive())
  414.                 reset(T.guts());
  415.         else
  416.                 V.reset(T.neg());
  417.         T.guts()->reset(this);
  418.         T -= V;
  419.        }
  420.         smShiftLeft(K);
  421.        }
  422.    }
  423.  
  424.     void bigIntegerG::      // lcm(I,J) = I * J / gcd(I, J)
  425. lcm (const integerG* I)
  426.    {if(!hasZero())
  427.         return;
  428.     else if(I->hasZero())
  429.         reset();
  430.     else
  431.        {bigIntegerA GCD = (bigIntegerA&)makeCopy(FALSE);
  432.         GCD.guts()->gcd(I);
  433.     divide(GCD.guts());
  434.     multiply(I);
  435.     abs();
  436.        }
  437.    }
  438.  
  439.  
  440. // bigInteger Operations //////////
  441.  
  442.  
  443.     void bigIntegerG::
  444. smGrowMag (unsigned int NewMax)
  445.    {unsigned long * NewMag = new unsigned long [NewMax];
  446.     for(unsigned int UI = 0; UI < SAM.size(); ++UI)
  447.         NewMag[UI] = SAM.magnitude(UI);
  448.     delete [] SAM.magnitude();
  449.     SAM.magnitude() = NewMag;
  450.     MaxSize = NewMax;
  451.    }
  452.  
  453.     void bigIntegerG::
  454. smReset (long Value)
  455.    {SAM.size(1);
  456.     SAM.positive(Value >= 0);
  457.     SAM.magnitude(0) = (Value >= 0) ? Value : -Value;
  458.    }
  459.  
  460.     void bigIntegerG::
  461. smReset (const signMagnitudeP& SMI)
  462.    {if(SMI.size() > MaxSize)
  463.        {MaxSize = SMI.size() + 4;
  464.         delete [] SAM.magnitude();
  465.     SAM.magnitude() = new unsigned long [MaxSize];
  466.        }
  467.     SAM.size(SMI.size());
  468.     SAM.positive(SMI.positive());
  469.     for(unsigned int UI = 0; UI < SAM.size(); ++UI)
  470.         SAM.magnitude(UI) = SMI.magnitude(UI);
  471.    }
  472.  
  473. // Shifts
  474.  
  475.     void bigIntegerG::
  476. smShiftLeft (unsigned int Count)
  477.    {unsigned int FullCount = Count / FULL;
  478.     unsigned int Shift     = Count % FULL;
  479.     if(FullCount)
  480.        {if(SAM.size()+FullCount+1 > MaxSize)
  481.            {MaxSize = SAM.size()+FullCount+4;
  482.         unsigned long* Mag = new unsigned long [MaxSize];
  483.         for(unsigned int UI = 0; UI < FullCount; ++UI)
  484.             Mag[UI] = 0;
  485.         for(UI = 0; UI < SAM.size(); ++UI)
  486.             Mag[UI+FullCount] = SAM.magnitude(UI);
  487.         SAM.size(UI + FullCount);
  488.         delete [] SAM.magnitude();
  489.         SAM.magnitude() = Mag;
  490.        }
  491.     else
  492.        {SAM.size(SAM.size() + FullCount);
  493.         for(unsigned int UI = SAM.size()-1; UI >= FullCount; --UI)
  494.             SAM.magnitude(UI) = SAM.magnitude(UI-FullCount);
  495.         for(; UI < FullCount; --UI)
  496.             SAM.magnitude(UI) = 0;
  497.        }
  498.        }
  499.     if(Shift)
  500.        {unsigned long Carry = 0;
  501.         for(unsigned int UI = 0; UI < SAM.size(); ++UI)
  502.        {unsigned long Next = SAM.magnitude(UI) >> (FULL-Shift);
  503.         SAM.magnitude(UI) <<= Shift;
  504.         SAM.magnitude(UI) |= Carry;
  505.         Carry = Next;
  506.        }
  507.     if(Carry)
  508.        {if(SAM.size() == MaxSize)
  509.             smGrowMag(MaxSize + 4);
  510.         SAM.magnitude(SAM.size()) = Carry;
  511.         SAM.incSize();
  512.        }
  513.        }
  514.    }
  515.  
  516.     void bigIntegerG::
  517. smShiftLeft (const signMagnitudeP& SMI)
  518.    {if(!SMI.zero() && !SAM.zero())
  519.        {if(!SMI.positive() || SMI.size() != 1 || SMI.magnitude(0) > UINT_MAX)
  520.            {ensure(FALSE, "throw rangeViolation");}
  521.         else
  522.            {smShiftLeft((unsigned int)SMI.magnitude(0));}
  523.        }
  524.    }
  525.  
  526.     void bigIntegerG::
  527. smShiftRight (unsigned int Count)
  528.    {unsigned int FullCount = Count / FULL;
  529.     if(FullCount >= SAM.size())
  530.         reset();
  531.     else
  532.        {if(FullCount)
  533.            {SAM.size(SAM.size() - FullCount);
  534.         for(unsigned int UI = 0; UI < SAM.size(); ++UI)
  535.             SAM.magnitude(UI) = SAM.magnitude(UI+FullCount);
  536.        }
  537.     unsigned int Shift = Count % FULL;
  538.         if(Shift)
  539.        {unsigned long Carry = 0;
  540.         for(unsigned int UI = SAM.size()-1; UI < SAM.size(); --UI)
  541.            {unsigned long Next = SAM.magnitude(UI) << (FULL-Shift);
  542.             SAM.magnitude(UI) >>= Shift;
  543.         SAM.magnitude(UI) |= Carry;
  544.         Carry = Next;
  545.            }
  546.         if(SAM.size() > 1 && !SAM.magnitude(SAM.size()-1))
  547.             SAM.decSize();
  548.        }
  549.        }
  550.    }
  551.  
  552.     void bigIntegerG::
  553. smShiftRight (const signMagnitudeP& SMI)
  554.    {if(!SMI.zero() && !SAM.zero())
  555.        {if(!SMI.positive())
  556.            {ensure(FALSE, "throw rangeViolation");}
  557.     else if(SMI.size() > 1 || SMI.magnitude(0) > UINT_MAX)
  558.            {reset();}
  559.     else
  560.            {smShiftRight((unsigned int)SMI.magnitude(0));}
  561.        }
  562.    }
  563.  
  564. // Comparison
  565.  
  566.     int bigIntegerG::
  567. smHasEqualMag (const signMagnitudeP& SMI) const
  568.    {if(SAM.size() != SMI.size())
  569.         return FALSE;
  570.     else
  571.        {for(unsigned int UI = 0; UI < SAM.size(); ++UI)
  572.            {if(SAM.magnitude(UI) != SMI.magnitude(UI))
  573.             return FALSE;
  574.        }
  575.     return TRUE;
  576.        }
  577.    }
  578.  
  579.     int bigIntegerG::
  580. smHasEqual (const signMagnitudeP& SMI) const
  581.    {if(SAM.positive() != SMI.positive())
  582.        {if(SAM.zero() && SMI.zero())
  583.             return TRUE;
  584.     else
  585.         return FALSE;
  586.        }
  587.     else    return smHasEqualMag(SMI);
  588.    }
  589.  
  590.     int bigIntegerG::
  591. smHasLessMag (const signMagnitudeP& SMI) const
  592.    {if(SAM.size() != SMI.size())
  593.        {if(SAM.size() > SMI.size())
  594.             return FALSE;
  595.     else
  596.             return TRUE;
  597.        }
  598.     else
  599.        {for(unsigned int UI = SAM.size()-1; UI < SAM.size(); --UI)
  600.            {if(SAM.magnitude(UI) < SMI.magnitude(UI))
  601.             return TRUE;
  602.         else if(SAM.magnitude(UI) > SMI.magnitude(UI))
  603.             return FALSE;
  604.        }
  605.     return FALSE;
  606.        }
  607.    }
  608.  
  609.     int bigIntegerG::
  610. smHasMoreMag (const signMagnitudeP& SMI) const
  611.    {if(SAM.size() != SMI.size())
  612.        {if(SAM.size() > SMI.size())
  613.             return TRUE;
  614.     else
  615.             return FALSE;
  616.        }
  617.     else
  618.        {for(unsigned int UI = SAM.size()-1; UI < SAM.size(); --UI)
  619.            {if(SAM.magnitude(UI) > SMI.magnitude(UI))
  620.             return TRUE;
  621.         else if(SAM.magnitude(UI) < SMI.magnitude(UI))
  622.             return FALSE;
  623.        }
  624.     return FALSE;
  625.        }
  626.    }
  627.  
  628.     int bigIntegerG::
  629. smHasLess (const signMagnitudeP& SMI) const
  630.    {if(SAM.positive() != SMI.positive())
  631.        {if(SAM.positive())
  632.             return FALSE;
  633.         else if(SAM.zero() && SMI.zero())
  634.             return FALSE;
  635.     else
  636.         return TRUE;
  637.        }
  638.     else
  639.        {if(SAM.positive())
  640.             return smHasLessMag(SMI);
  641.     else
  642.         return smHasMoreMag(SMI);
  643.        }
  644.    }
  645.  
  646.     int bigIntegerG::
  647. smHasMore (const signMagnitudeP& SMI) const
  648.    {if(SAM.positive() != SMI.positive())
  649.        {if(!SAM.positive())
  650.             return FALSE;
  651.         else if(SAM.zero() && SMI.zero())
  652.             return FALSE;
  653.     else
  654.         return TRUE;
  655.        }
  656.     else
  657.        {if(!SAM.positive())
  658.             return smHasLessMag(SMI);
  659.     else
  660.         return smHasMoreMag(SMI);
  661.        }
  662.    }
  663.  
  664. // Increment
  665.  
  666.     void bigIntegerG::
  667. smIncExtend ()
  668.    {if(SAM.size() == MaxSize)
  669.         smGrowMag(MaxSize + 4);
  670.     SAM.magnitude(SAM.size()) = 1;
  671.     SAM.incSize();
  672.    }
  673.         
  674.     void bigIntegerG::
  675. smIncMag ()
  676.    {unsigned long Carry = 1;
  677.     for(unsigned int UI = 0; UI < SAM.size() && Carry; ++UI)
  678.         Carry = !(SAM.magnitude(UI) += 1);
  679.     if(Carry)
  680.         smIncExtend();
  681.    }
  682.  
  683.     void bigIntegerG::
  684. smIncMag (unsigned long UL)
  685.    {unsigned long Carry = SAM.magnitude(0) > ~UL;
  686.     SAM.magnitude(0) += UL;
  687.     for(unsigned int UI = 1; Carry && UI < SAM.size(); ++UI)
  688.         Carry = !(SAM.magnitude(UI) += 1);
  689.     if(Carry)
  690.         smIncExtend();
  691.    }
  692.  
  693.     void bigIntegerG::
  694. smIncMag (const signMagnitudeP& SMI)
  695.    {unsigned long Carry = 0;
  696.     unsigned int Size = (SAM.size() < SMI.size()) ? SAM.size() : SMI.size();
  697.     for(unsigned int UI = 0; UI < Size; ++UI)
  698.        {if(Carry)
  699.            {Carry = (SAM.magnitude(UI) >= ~SMI.magnitude(UI));
  700.         SAM.magnitude(UI) += SMI.magnitude(UI) + 1;
  701.        }
  702.     else
  703.        {Carry = (SAM.magnitude(UI) > ~SMI.magnitude(UI));
  704.         SAM.magnitude(UI) += SMI.magnitude(UI);
  705.        }
  706.        }
  707.     if(SMI.size() > Size)
  708.        {if(SMI.size() + Carry > MaxSize)
  709.             smGrowMag(SMI.size() + 4);
  710.     SAM.size(SMI.size());
  711.     for(; UI < SAM.size(); ++UI)
  712.        {SAM.magnitude(UI) = SMI.magnitude(UI) + Carry;
  713.         Carry = Carry && !SAM.magnitude(UI);
  714.        }
  715.     if(Carry)
  716.        {SAM.magnitude(SAM.size()) = 1;
  717.         SAM.incSize();
  718.        }
  719.        }
  720.     else
  721.        {for(; Carry && UI < SAM.size(); ++UI)
  722.             Carry = !(SAM.magnitude(UI) += 1);
  723.     if(Carry)
  724.         smIncExtend();
  725.        }
  726.    }
  727.  
  728.     void bigIntegerG::
  729. smInc (const signMagnitudeP& SMI)
  730.    {if(SAM.zero())
  731.         smReset(SMI);
  732.     else if(SAM.positive() == SMI.positive())
  733.        {if(&SAM == &SMI) smShiftLeft(1);
  734.     else         smIncMag(SMI);
  735.        }
  736.     else          smDecMag(SMI);
  737.    }
  738.  
  739. // Decrement
  740.  
  741.     void bigIntegerG::
  742. smDecMag ()
  743.    {if(SAM.zero())
  744.        {SAM.magnitude(0) = 1;
  745.         SAM.positive(0);
  746.        }
  747.     else
  748.        {unsigned long Carry = 1;
  749.         for(unsigned int UI = 0; Carry && UI < SAM.size(); ++UI)
  750.        {Carry = !SAM.magnitude(UI);
  751.         SAM.magnitude(UI) -= 1;
  752.        }
  753.     for(UI = SAM.size()-1; !SAM.magnitude(UI) && UI > 0; --UI);
  754.     SAM.size(UI + 1);
  755.        }
  756.    }
  757.  
  758.     void bigIntegerG::
  759. smDecMag (unsigned long UL)
  760.    {if(SAM.size() > 1)
  761.        {unsigned long Carry = (SAM.magnitude(0) < UL ? 1 : 0);
  762.         SAM.magnitude(0) -= UL;
  763.     if(Carry)
  764.        {for(unsigned int UI = 1; Carry && UI < SAM.size(); ++UI)
  765.            {Carry = !SAM.magnitude(UI);
  766.             SAM.magnitude(UI) -= 1;
  767.            }
  768.         for(UI = SAM.size()-1; !SAM.magnitude(UI) && UI > 0; --UI);
  769.         SAM.size(UI + 1);
  770.        }
  771.        }
  772.     else
  773.        {if(SAM.magnitude(0) >= UL)
  774.             SAM.magnitude(0) -= UL;
  775.     else
  776.        {SAM.magnitude(0) = UL - SAM.magnitude(0);
  777.         SAM.negate();
  778.        }
  779.        }
  780.    }
  781.  
  782.     void bigIntegerG::
  783. smDecMag (const signMagnitudeP& SMI)
  784.    {unsigned long Carry = 0;
  785.     if(smHasLessMag(SMI))
  786.        {SAM.negate();
  787.         if(SMI.size() > MaxSize)
  788.             smGrowMag(SMI.size() + 4);
  789.     for(unsigned int UI = 0; UI < SAM.size(); ++UI)
  790.            {if(Carry)
  791.            {Carry = SMI.magnitude(UI) <= SAM.magnitude(UI);
  792.             SAM.magnitude(UI) = SMI.magnitude(UI) - SAM.magnitude(UI) - 1;
  793.            }
  794.         else
  795.            {Carry = SMI.magnitude(UI) < SAM.magnitude(UI);
  796.             SAM.magnitude(UI) = SMI.magnitude(UI) - SAM.magnitude(UI);
  797.            }
  798.        }
  799.     SAM.size(SMI.size());
  800.     for(; UI < SMI.size(); ++UI)
  801.        {if(Carry)
  802.            {Carry = !SMI.magnitude(UI);
  803.             SAM.magnitude(UI) = SMI.magnitude(UI) - 1;
  804.            }
  805.         else
  806.             SAM.magnitude(UI) = SMI.magnitude(UI);
  807.            }
  808.     if(SAM.size() > 1 && !SAM.magnitude(SAM.size()-1))
  809.         SAM.decSize();
  810.        }
  811.     else
  812.        {for(unsigned int UI = 0; UI < SMI.size(); ++UI)
  813.            {if(Carry)
  814.            {Carry = SAM.magnitude(UI) <= SMI.magnitude(UI);
  815.             SAM.magnitude(UI) -= SMI.magnitude(UI) + 1;
  816.            }
  817.         else
  818.            {Carry = SAM.magnitude(UI) < SMI.magnitude(UI);
  819.             SAM.magnitude(UI) -= SMI.magnitude(UI);
  820.            }
  821.        }
  822.     for(; Carry; ++UI)
  823.        {Carry = !SAM.magnitude(UI);
  824.         SAM.magnitude(UI) -= 1;
  825.        }
  826.     for(UI = SAM.size()-1; !SAM.magnitude(UI) && UI > 0; --UI);
  827.     SAM.size(UI + 1);
  828.        }
  829.    }
  830.  
  831.     void bigIntegerG::
  832. smDec (const signMagnitudeP& SMI)
  833.    {if(SAM.zero())
  834.        {smReset(SMI);
  835.         SAM.negate();
  836.        }
  837.     else if(SAM.positive() != SMI.positive())
  838.        {if(&SAM == &SMI) smShiftLeft(1);
  839.         else         smIncMag(SMI);
  840.        }
  841.     else          smDecMag(SMI);
  842.    }
  843.  
  844. // Multiplication
  845.  
  846.     // This algorithm is based on hand multiplication, modified
  847.     // to prevent the need to add carries within partial products
  848.     // by producing twice the partial products by multiplying by
  849.     // every other digit:
  850.     //            5  6  7  8                     5 6  7 8
  851.     //          * 1  2  3  4                   * 1 2  3 4
  852.     //          ------------                   ----------
  853.     //             4*6   4*8                     2 4  3 2
  854.     //          4*5   4*7                     2  0 2  8
  855.     //            3*6   3*8                     1  8 2  4
  856.     //         3*5   3*7                      1 5  2 1
  857.     //       2*6   2*8                1 2  1 6
  858.     //    2*5   2*7                      1  0 1  4
  859.     //    1*6   1*8                      0  6 0  8
  860.     // 1*5   1*7                       0 5  0 7
  861.     // ---------------------           ------------------
  862.     //   7  0  0  6  6  5  2             7  0 0  6 6  5 2
  863.  
  864.     void bigIntegerG::
  865. smMulMag (const signMagnitudeP& Mplier)
  866.    {signMagnitudeP Mcand (SAM);
  867.     signMagnitudeP Part  (SAM.size(), 1, new unsigned long [SAM.size()]);
  868.     MaxSize = 4 + Mcand.size() + Mplier.size();
  869.     SAM.magnitude() = new unsigned long [MaxSize];
  870.     SAM.magnitude(0) = 0;
  871.     SAM.size(1);
  872.     for(unsigned int UI = Mplier.size()-1; UI < Mplier.size(); --UI)
  873.        {unsigned long R = HI(Mplier.magnitude(UI));
  874.         for(unsigned int UJ = 0; UJ < Mcand.size(); ++UJ)
  875.         Part.magnitude(UJ) = HI(Mcand.magnitude(UJ)) * R;
  876.         smIncMag(Part);
  877.     smShiftLeft(HALF);
  878.     for(UJ = 0; UJ < Mcand.size(); ++UJ)
  879.         Part.magnitude(UJ) = LO(Mcand.magnitude(UJ)) * R;
  880.         smIncMag(Part);
  881.     R = LO(Mplier.magnitude(UI));
  882.         for(UJ = 0; UJ < Mcand.size(); ++UJ)
  883.         Part.magnitude(UJ) = HI(Mcand.magnitude(UJ)) * R;
  884.         smIncMag(Part);
  885.     smShiftLeft(HALF);
  886.     for(UJ = 0; UJ < Mcand.size(); ++UJ)
  887.         Part.magnitude(UJ) = LO(Mcand.magnitude(UJ)) * R;
  888.         smIncMag(Part);
  889.        }
  890.     delete [] Mcand.magnitude();
  891.     delete [] Part.magnitude();
  892.    }
  893.  
  894.     void bigIntegerG::
  895. smMul (const signMagnitudeP& SMI)
  896.    {if(!SAM.zero())
  897.        {if(SMI.zero())
  898.             reset();
  899.     else
  900.        {if(!SMI.positive()) neg();
  901.         smMulMag(SMI);
  902.        }
  903.        }
  904.    }
  905.  
  906. // Division
  907.  
  908.     inline unsigned long bigIntegerG::
  909. smDivRem10 ()
  910.    {unsigned long M = 0;
  911.     for(unsigned int UI = SAM.size()-1; UI < SAM.size(); --UI)
  912.        {M <<= HALF;
  913.         M |= HI(SAM.magnitude(UI));
  914.     unsigned long H = M / 10;
  915.     M %= 10;
  916.     M <<= HALF;
  917.     M |= LO(SAM.magnitude(UI));
  918.     SAM.magnitude(UI) = (H << HALF) | (M / 10);
  919.     M %= 10;
  920.        }
  921.     if(SAM.size() > 1 && !SAM.magnitude(SAM.size()-1))
  922.         SAM.decSize();
  923.     return M;
  924.    }
  925.  
  926.     // This is based on Knuth, Art of Computer Programming, Vol 2
  927.     // Algorithm D (division of nonnegative integers).
  928.  
  929.     void bigIntegerG::
  930. smRemDivMag (const signMagnitudeP& Dsor,
  931.          signMagnitudeP& Quot, unsigned int& QuotMaxSize)
  932.    {// Special case: quot = 0, rem = this
  933.     if(smHasLessMag(Dsor))
  934.        {if(!QuotMaxSize)
  935.            {QuotMaxSize = 1;
  936.         Quot.magnitude() = new unsigned long [1];
  937.        }
  938.     Quot.magnitude(0) = 0;
  939.     Quot.size(1);
  940.         return;
  941.        }
  942.  
  943.     // Normalize (shift Dsor to have msb in bit 31 of a word)
  944.     unsigned int Shift = FULL - log2(Dsor.magnitude(Dsor.size()-1)) - 1;
  945.     smShiftLeft(Shift);
  946.     // special case: quot = 1, rem = this-Dsor
  947.     if(SAM.size() == Dsor.size())
  948.        {smDecMag(Dsor);
  949.         if(!QuotMaxSize)
  950.            {QuotMaxSize = 1;
  951.         Quot.magnitude() = new unsigned long [1];
  952.        }
  953.     Quot.magnitude(0) = 1;
  954.     Quot.size(1);
  955.         return;
  956.        }
  957.     signMagnitudeP NDsor (Dsor);
  958.     NDsor.magnitude() = new unsigned long [Dsor.size()];
  959.     NDsor.magnitude(0) = Dsor.magnitude(0) << Shift;
  960.     unsigned long Carry = Dsor.magnitude(0) >> (FULL-Shift);
  961.     for(unsigned int UI = 1; UI < Dsor.size(); ++UI)
  962.        {NDsor.magnitude(UI) = (Dsor.magnitude(UI) << Shift) | Carry;
  963.     Carry = Dsor.magnitude(UI) >> (FULL-Shift);
  964.        }
  965.     assumed(!Carry, "internal error in normalization for division");
  966.     unsigned long Divisor    = HI(NDsor.magnitude(NDsor.size()-1));
  967.     unsigned long SubDivisor = LO(NDsor.magnitude(NDsor.size()-1));
  968.  
  969.     // Initialize Quotient
  970.     unsigned int QuotMax = SAM.size() - NDsor.size() + 1;
  971.     Quot.size(QuotMax);
  972.     if(QuotMax > QuotMaxSize)
  973.        {QuotMaxSize = QuotMax;
  974.         delete [] Quot.magnitude();
  975.     Quot.magnitude() = new unsigned long [QuotMax];
  976.        }
  977.  
  978.     // Initialize Product
  979.     signMagnitudeP Prod (SAM.size(), 1, new unsigned long [SAM.size()]);
  980.     for(UI = 0; UI < QuotMax; ++UI)
  981.         Prod.magnitude(UI) = 0;
  982.  
  983.     // First Iteration of Divide Loop (Quotient must be 0 or 1 due to norm.)
  984.     for(UI = 0; UI < NDsor.size(); ++UI)
  985.         Prod.magnitude(UI + QuotMax - 1) = NDsor.magnitude(UI);
  986.     if(smHasLessMag(Prod))
  987.        {Quot.magnitude(QuotMax-1) = 0;}
  988.     else
  989.        {smDecMag(Prod);
  990.         Quot.magnitude(QuotMax-1) = 1;
  991.        }
  992.  
  993.     // Subsequent Iterations of Divide Loop
  994.     unsigned int Halves = (QuotMax-1) << 1;
  995.     for(UI = Halves-1; UI < Halves; --UI)
  996.        {unsigned int QIndex = UI >> 1;
  997.         unsigned int QHalf  = UI & 1U;
  998.  
  999.         // Estimate Quotient
  1000.     unsigned long Dividend;
  1001.     unsigned long SubDividend;
  1002.     if(QHalf)
  1003.        {Dividend    =    SAM  .magnitude(QIndex + NDsor.size());
  1004.         SubDividend = HI(SAM  .magnitude(QIndex + NDsor.size() - 1));
  1005.        }
  1006.     else
  1007.        {Dividend    =   UP(SAM.magnitude(QIndex + NDsor.size()))
  1008.                   | HI(SAM.magnitude(QIndex + NDsor.size() - 1));
  1009.         SubDividend = LO(SAM  .magnitude(QIndex + NDsor.size() - 1));
  1010.        }
  1011.     unsigned long Quotient = Dividend / Divisor;
  1012.     if(Quotient > LOWER)
  1013.         Quotient = LOWER;
  1014.  
  1015.     // Heuristic to make missed guesses more rare
  1016.     // Knuth's algorithm calls for this test to be repeated,
  1017.     // but that leads to wrong answers!
  1018.     if((Quotient*SubDivisor) > (UP(Dividend-Quotient*Divisor)|SubDividend))
  1019.        {--Quotient;}
  1020.  
  1021.     if(Quotient)
  1022.        {// Multiply
  1023.         Carry = 0;
  1024.         for(unsigned int UJ = 0; UJ < NDsor.size(); ++UJ)
  1025.            {unsigned long ProdLo = Quotient * LO(NDsor.magnitude(UJ))
  1026.                        + Carry;
  1027.             unsigned long ProdHi = Quotient * HI(NDsor.magnitude(UJ)) 
  1028.                        + HI(ProdLo);
  1029.         Carry = HI(ProdHi);
  1030.         if(QHalf)
  1031.            {Prod.magnitude(UJ+QIndex)  |= UP(ProdLo);
  1032.             Prod.magnitude(UJ+QIndex+1) = LO(ProdHi);
  1033.            }
  1034.         else
  1035.            {Prod.magnitude(UJ+QIndex) = UP(ProdHi) | LO(ProdLo);}
  1036.            }
  1037.         if(QHalf)
  1038.             Prod.magnitude(NDsor.size() + QIndex) |= UP(Carry);
  1039.         else
  1040.             Prod.magnitude(NDsor.size() + QIndex)  = LO(Carry);
  1041.         Prod.size(NDsor.size()+QIndex+1);
  1042.         for(UJ = Prod.size()-1; !Prod.magnitude(UJ) && UJ > 0; --UJ);
  1043.         Prod.size(UJ + 1);
  1044.  
  1045.         // Compare (rare: this is true 1/32676 of the time, per Knuth)
  1046.         while(smHasLessMag(Prod))
  1047.            {--Quotient;
  1048.             // reduce Prod by NDsor
  1049.         Carry = 0;
  1050.         for(UJ = 0; UJ < NDsor.size(); ++UJ)
  1051.            {if(QHalf)
  1052.                    {unsigned long Temp = UP(Prod.magnitude(UJ+QIndex+1))
  1053.                        | HI(Prod.magnitude(UJ+QIndex));
  1054.                 if(Carry)
  1055.                {Carry = Temp <= NDsor.magnitude(UJ);
  1056.                 Temp -= NDsor.magnitude(UJ) + 1;
  1057.                }
  1058.             else
  1059.                {Carry = Temp <  NDsor.magnitude(UJ);
  1060.                 Temp -= NDsor.magnitude(UJ);
  1061.                }
  1062.             Prod.magnitude(UJ+QIndex) &= LOWER;
  1063.             Prod.magnitude(UJ+QIndex) |= UP(Temp);
  1064.             Prod.magnitude(UJ+QIndex+1) &= ~LOWER;
  1065.             Prod.magnitude(UJ+QIndex+1) |= HI(Temp);
  1066.                }
  1067.             else
  1068.                {if(Carry)
  1069.                    {Carry = Prod.magnitude(UJ + QIndex)
  1070.                     <= NDsor.magnitude(UJ);
  1071.                 Prod.magnitude(UJ+QIndex) -= NDsor.magnitude(UJ)+1;
  1072.                }
  1073.             else
  1074.                {Carry = Prod.magnitude(UJ+QIndex)
  1075.                     < NDsor.magnitude(UJ);
  1076.                 Prod.magnitude(UJ+QIndex) -= NDsor.magnitude(UJ);
  1077.                }
  1078.                }
  1079.            }
  1080.         if(Carry)
  1081.            {if(QHalf)
  1082.                     Prod.magnitude(NDsor.size() + QIndex) -= UP(Carry);
  1083.             else
  1084.             Prod.magnitude(NDsor.size() + QIndex) -= LO(Carry);
  1085.            }
  1086.         for(UJ = Prod.size()-1; !Prod.magnitude(UJ) && UJ > 0; --UJ);
  1087.         Prod.size(UJ + 1);
  1088.            }
  1089.  
  1090.         // Subtract
  1091.         smDecMag(Prod);
  1092.        }
  1093.  
  1094.     // Record Quotient
  1095.     if(QHalf)
  1096.        {Quot.magnitude(QIndex) = UP(Quotient);}
  1097.     else
  1098.        {Quot.magnitude(QIndex) |= Quotient;}
  1099.        }
  1100.  
  1101.     for(UI = Quot.size()-1; !Quot.magnitude(UI) && UI > 0; --UI);
  1102.     Quot.size(UI + 1);
  1103.  
  1104.     // Unnormalize
  1105.     if(Shift)
  1106.         smShiftRight(Shift);
  1107.     delete [] Prod.magnitude();
  1108.     delete [] NDsor.magnitude();
  1109.    }
  1110.  
  1111.     void bigIntegerG::
  1112. smDiv (const signMagnitudeP& SMI)
  1113.    {signMagnitudeP Quotient;
  1114.     unsigned int QMaxSize = 0;
  1115.     smRemDivMag(SMI, Quotient, QMaxSize);
  1116.     if(!SMI.positive()) SAM.negate();
  1117.     delete [] SAM.magnitude();
  1118.     SAM.magnitude() = Quotient.magnitude();
  1119.     SAM.size(Quotient.size());
  1120.     MaxSize = QMaxSize;
  1121.    }
  1122.  
  1123. // Output
  1124.  
  1125.     char * bigIntegerG::
  1126. string (char* Buffer, int N, char Sep) const
  1127.    {assumed(N > 2, "string() called with N <= 2");
  1128.     Buffer[--N] = '\0';
  1129.     if(hasZero())
  1130.         Buffer[--N] = '0';
  1131.     else
  1132.        {bigIntegerA I = (bigIntegerA&)makeCopy(FALSE);
  1133.         if(Sep)
  1134.        {int Count = 3;
  1135.         while(!I.hasZero() && N)
  1136.            {if(!Count)
  1137.                {Buffer[--N] = Sep;
  1138.             Count = 3;
  1139.            }
  1140.         Buffer[--N] = '0' + char(I.guts()->smDivRem10());
  1141.         --Count;
  1142.            }
  1143.         if(hasNegative() && N)
  1144.             Buffer[--N] = '-';
  1145.        }
  1146.     else
  1147.        {while(!I.hasZero() && N)
  1148.             Buffer[--N] = '0' + char(I.guts()->smDivRem10());
  1149.         if(hasNegative() && N)
  1150.             Buffer[--N] = '-';
  1151.        }
  1152.        }
  1153.     return &Buffer[N];
  1154.    }
  1155.  
  1156. //***************************************************************************
  1157.