home *** CD-ROM | disk | FTP | other *** search
/ Collection of Education / collectionofeducationcarat1997.iso / SCIENCE / PIW131.ZIP / MULTIFDW.CPP next >
Encoding:
C/C++ Source or Header  |  1994-07-29  |  36.5 KB  |  1,407 lines
  1. // ------- MultiFDW.Cpp  Multiple-precision floating decimal algorithms
  2. /*
  3.   Version : Version 3.11, last revised: 1994-07-29, 0600 hours
  4.   Author  : Copyright (c) 1981-1994 by author: Harry J. Smith,
  5.   Address : 19628 Via Monte Dr., Saratoga, CA 95070.  All rights reserved.
  6. */
  7.  
  8. #include "MultiFDW.h" // Multiple-precision Floating decimal algorithms module
  9. #include <limits.h>
  10.  
  11. // ------- Convert +Double to a String integer ?????? is this needed ??????
  12. void StrDI(char* St, double D);
  13.  
  14. // Developed in Turbo Pascal 5.5.  Converted to Borland C++ 3.1 for Windows.
  15.  
  16. // ------- The following are set when InitMultiF is called
  17. Mu1Digit FRound;   // Floating point rounding constant, Base / 2 or 0
  18. long     FMC;      // Current max # of super digits in a result
  19. long     FTn;      // Current # of super digits to truncate display
  20. long     FDn;      // Current max decimal digits to display
  21. Bool     ScieN;    // Force scientific notation on flag
  22. Bool     Expanded; // True if running in expanded precision
  23. long     FMCB;     // Base FMC when running expanded precision
  24.  
  25. // ------- Implementation
  26.  
  27. // ------- Common hidden routines
  28. //         none
  29. // ------- End of common hidden routines
  30.  
  31. // ------- MultiF's method implementation
  32.  
  33. // ------- this = X
  34. void MultiF::SetTo( MultiF& X)
  35. {
  36.   long   I, J, K;
  37.   double D;
  38.   Bool   RoundIt;
  39.   long   L;
  40.  
  41.   K = MinL(M, FMC);
  42.   if (K >= X.N) {
  43.     MultiSI::SetTo(X);
  44.     C = X.C;  return;
  45.   }
  46.   N = K;  S = X.S;
  47.   J = X.N - K;
  48.   for ( I = 0; I < K; I++)  V[I] = X.V[I + J];
  49.   D = X.C + J;
  50.   if (D > LONG_MAX) {
  51.     Clear();  MuWriteErr("Floating SetTo overflow, continuing...");
  52.   }
  53.   else {
  54.     C = X.C + J;
  55.     if ((FRound > 0) && (X.V[J - 1] >= FRound)) {
  56.       RoundIt = True;
  57.       if ((X.V[J - 1] == FRound) && (J == 1)) {
  58.     L = X.V[J];
  59.     if (!(L & 1))  RoundIt = False; // if not odd
  60.       }
  61.       if (RoundIt) {
  62.         MultiF T(1);
  63.     T.V[0] = 1;
  64.     T.S = S;
  65.     T.C = C;
  66.         RAdd(T);
  67.       }
  68.     }
  69.     Norm();
  70.   }
  71. } // --end-- MultiF::SetTo
  72.  
  73. // ------- Normalize
  74. void MultiF::Norm()
  75. {
  76.   long     I, J;
  77.   Mu1Digit SaveFR;
  78.   long     SaveFMC;
  79.   double   D;
  80.   long     L;
  81.   Bool     RoundIt;
  82.  
  83.   MultiSI::Norm();
  84.   if (ZTest()) { C = 0;  return; }    
  85.   I = 0;
  86.   if (N > FMC) {
  87.     I = N - FMC;
  88.     D = C + I;
  89.     if (D > LONG_MAX) {
  90.       MuWriteErr("Floating Norm overflow, continuing...");
  91.       Clear();  return;
  92.     }
  93.     if ((FRound > 0) && (V[I - 1] >= FRound)) {
  94.       RoundIt = True;
  95.       if (V[I - 1] == FRound) {
  96.     L = V[I];
  97.     if (!(L & 1)) { // if not odd
  98.       RoundIt = False;
  99.       J = I - 2;
  100.       while (J >= 0) {
  101.         if (V[J] != 0) {
  102.           J = 0;  RoundIt = True;
  103.         }
  104.             J--;
  105.       }
  106.     }
  107.       }
  108.       if (RoundIt) {
  109.     SaveFR = FRound;
  110.     FRound = 0;
  111.     SaveFMC = FMC;
  112.     FMC = M;
  113.     MultiF T(1);
  114.     T.V[0] = 1;
  115.     T.S = S;
  116.     T.C = C + I;
  117.     if (V[0] == 0)  V[0] = 1; // Prevent Normalizing early
  118.     RAdd(T);
  119.     FRound = SaveFR;
  120.     FMC = SaveFMC;
  121.       }
  122.     }
  123.     N -= I;  C += I;
  124.   }
  125.   while (V[I] == 0) {
  126.     I++;  // I = # of trailing zero super digits
  127.     N = N - 1;
  128.     D = C + 1.0;
  129.     if (D > LONG_MAX) {
  130.       MuWriteErr("Floating Norm overflow, continuing...");
  131.       Clear();  return;
  132.     }
  133.     C = C + 1;
  134.   }
  135.   if (I > 0) {
  136.     for (J = 0; J < N; J++) {
  137.       V[J] = V[I];
  138.       I++;
  139.     }
  140.   }
  141. } // --end-- MultiF::Norm
  142.  
  143. // ------- this = this, R "digits" lost, unnormalized
  144. void MultiF::ShiftR( long R)
  145. // R < 0 => Shift left, Sets MuErr True if overflow
  146. {
  147.   long I;
  148.  
  149.   if (R == 0)  return;
  150.   if (R  < 0) { ShiftL(-R);  return; }
  151.   if (R >= N) { Clear();  return; }
  152.   if ((FRound > 0) && (!V[R] || (V[R] == FRound))) {
  153.     I = 0;
  154.     while (I < R) {
  155.       if (V[I] != 0) {
  156.     I = R;  V[R] = V[R] + 1; // Add sticky bit to meet IEEE standard
  157.       }
  158.       I++;
  159.     }
  160.   }
  161.   for (I = R; I < N; I++)  V[I - R] = V[I];
  162.   N = N - R;
  163.   if (C > (LONG_MAX - R)) {
  164.     Clear();  MuWriteErr("Floating shift right underflow, continuing...");
  165.   }
  166.   else
  167.     C = C + R;
  168. } // --end-- MultiF::ShiftR
  169.  
  170. // ------- this = this,L "digits" more, unnormalized
  171. void MultiF::ShiftL( long L)
  172. // R < 0 => Shift right, Sets MuErr True if overflow
  173. {
  174.   long I;
  175.  
  176.   if (L == 0)  return;
  177.   if (L < 0) { ShiftR(-L);  return; }
  178.   if ((N + L) > M) {
  179.     MuWriteErr("Floating shift left overflow, continuing...");
  180.     Clear();  return;
  181.   }
  182.   for (I = (N - 1); I >= 0; I--)  V[I + L] = V[I];
  183.   for (I = 0; I < L; I++)  V[I] = 0;
  184.   N = N + L;
  185.   if (C < (L - LONG_MAX))  Clear();
  186.   else
  187.     C = C - L;
  188. } // --end-- MultiF::ShiftL
  189.  
  190. // ------- this = this + X
  191. void MultiF::RAdd( MultiF& X)
  192. // Sets MuErr True if overflow
  193. {
  194.   Add( *this, X);
  195. } // --end-- MultiF::RAdd
  196.  
  197. // ------- this = this - X
  198. void MultiF::RSub( MultiF& X)
  199. // Sets MuErr True if overflow
  200. {
  201.   if (V == X.V)
  202.     Clear();
  203.   else {
  204.     X.S = !X.S;
  205.     RAdd(X);
  206.     X.S = !X.S;
  207.   }
  208. } // --end-- MultiF::RSub
  209.  
  210. // ------- this = X + Y
  211. void MultiF::Add( MultiF& X, MultiF& Y)
  212. // Sets MuErr True if overflow
  213. {
  214.   double   XF,  YF,  TF; // Length of mantissa
  215.   MultiF  *XL, *YL, *TL; // Local pointers
  216.   long    MaxS;
  217.  
  218.   if (X.ZTest()) { SetTo(Y);  return; }
  219.   if (Y.ZTest()) { SetTo(X);  return; }
  220.   XL = &X;  XF = (double)(X.C) + X.N;
  221.   YL = &Y;  YF = (double)(Y.C) + Y.N;
  222.   if (YF > XF) {
  223.     TL = XL;  XL = YL;  YL = TL;
  224.     TF = XF;  XF = YF;  YF = TF;
  225.   }
  226.   if ((XF - YF - 3) >= M)  SetTo( *XL);
  227.   else {
  228.     XF = XF - MinL( X.C, Y.C);
  229.     MaxS = MinD( XF, M) + 1 + 3; // 3 guard + 1 overflow "Digits"
  230.     MultiF XT( MaxS);  XT.SetTo( *XL);
  231.     MultiF YT( MaxS);  YT.SetTo( *YL); // XT, YT are Temp on heap
  232.     XT.ShiftL( XT.M - XT.N - 1);
  233.     YT.ShiftR( XT.C - YT.C);
  234.     XT.MultiSI::RAdd( YT);
  235.     XT.Norm();
  236.     SetTo( XT);
  237.   }
  238. } // --end-- MultiF::Add
  239.  
  240. // ------- this = X - Y
  241. void MultiF::Sub( MultiF& X, MultiF& Y)
  242. // Sets MuErr True if overflow
  243. {
  244.   if (&X == &Y)  Clear();
  245.   else if (V == X.V) 
  246.     RSub(Y);
  247.   else if (V == Y.V) {
  248.     RSub(X);  ChS();
  249.   }
  250.   else {
  251.     Y.S = !Y.S;
  252.     Add(X, Y);
  253.     Y.S = !Y.S;
  254.   }
  255. } // --end-- MultiF::Sub
  256.  
  257. // -------
  258. void MultiF::Value( char* St, int& I)
  259. // Convert String to this, returns I = # of character used
  260. // Allows St to be a literal, MuErr set True if overflow
  261. {
  262.   int      J, K, L;
  263.   Mu1Digit D1;
  264.   double   Db, DI;
  265.   Bool     Sign;
  266.  
  267.   Sign = (St[0] == '-');
  268.   MultiSI::Value( St, J);
  269.   I= J;
  270.   C = 0;
  271.   Norm();
  272.   if (J) strcpy( St, St+J); // Delete first J characters
  273.   if (strlen( St) == 0)  return;
  274.   if ((strlen( St) >= 2) && (St[0] == '.') &&
  275.      (St[1] >= '0') && (St[1] <= '9')) {
  276.     strcpy( St, St+1); // Delete first character
  277.     MultiF X(1 + strlen( St) / MuDMax); // Temp on heap
  278.     X.MultiSI::Value( St, J);
  279.     I += J + 1;
  280.     K = 0;
  281.     for (L = 0; L < J; L++)  if ((St[L] >= '0') && (St[L] <= '9'))  K++;
  282.     X.C = -((K - 1) / MuDMax) - 1;
  283.     K = MuDMax - ((K - 1) % MuDMax) - 1;
  284.     D1 = 1;
  285.     for (L = 0; L < K; L++)  D1 *= 10;
  286.     X.MultiSI::RMul1( D1);
  287.     X.S = Sign;
  288.     RAdd(X);
  289.     if (J) strcpy( St, St+J); // Delete first J characters
  290.   }
  291.   if (strlen( St) < 2)  return;
  292.  
  293.   if ( ((St[0] == 'E') || (St[0] == 'e')) &&    (
  294.      ((St[1] >= '0') && (St[1] <= '9')) ||
  295.      ((strlen(St) >= 3) &&
  296.       ((St[1] == '+') || (St[1] == '-')) &&
  297.       (St[2] >= '0') && (St[2] <= '9'))     ) ) {
  298.  
  299.     strcpy( St, St+1); // Delete first character
  300.     MultiF X(1 + strlen( St) / MuDMax); // Temp on heap
  301.     X.MultiSI::Value( St, J);
  302.     I += J + 1;
  303.     X.MultiSI::GetD( Db);
  304.     modf( Db / MuDMax, &DI); // DI = Integer part
  305.     if (Db < 0.0)  DI--;
  306.     if ((DI + C) > LONG_MAX) {
  307.       MuWriteErr("Floating Value overflow, continuing...");
  308.       Clear();  return;
  309.     }
  310.     if ((DI + C) < -LONG_MAX) {
  311.       Clear();  return;
  312.     }
  313.     C = DI + C;
  314.     K = Db - DI * MuDMax;
  315.     D1 = 1;
  316.     for (L = 1; L <= K; L++)  D1 *= 10;
  317.     RMul1( D1);
  318.   }
  319.   Norm();
  320. } // --end-- MultiF::Value
  321.  
  322. // ------- this = D1 Mod Base
  323. void MultiF::SetTo1( Mu1Digit D1)
  324. // Allows a literal D1, D1 an integer, |D1| < Base, D1 not changed
  325. {
  326.   MultiSI::SetTo1( D1);
  327.   C = 0;
  328. } // --end-- MultiF::SetTo1
  329.  
  330. // ------- this = Db Mod Base**4
  331. void MultiF::SetToD( double Db)
  332. // Allows a literal Db, Db not changed, MuErr set True if overflow
  333. {
  334.   MultiSI::SetToD( Db);
  335.   C = 0;
  336.   Norm();
  337. } // --end-- MultiF::SetToD
  338.  
  339. // ------- D1 = this Mod Base
  340. void MultiF::Get1( Mu1Digit& D1)
  341. // this is not changed
  342. {
  343.   Bool    SaveB;
  344.  
  345.   SaveB = MuErrRep;  MuErrRep = False; // Suppress error messages
  346.   MultiF X(5);  X.SetTo( *this);       // X is a temp on the heap
  347.   if (X.C > 4) {
  348.     MuErrRep = SaveB;
  349.     X.Clear();
  350.     MuWriteErr("Get1 overflow, continuing...");
  351.   }
  352.   else  X.ShiftL(X.C);
  353.   X.MultiSI::Get1( D1);
  354.   MuErrRep = SaveB;
  355. } // --end-- MultiF::Get1
  356.  
  357. // ------- Db = this Mod Base**4
  358. void MultiF::GetD( double& Db)
  359. // this is not changed
  360. {
  361.   Bool    SaveB;
  362.  
  363.   SaveB = MuErrRep;  MuErrRep = False; // Suppress error messages
  364.   MultiF X(5);  X.SetTo( *this);
  365.   if (X.C > 4) {
  366.     MuErrRep = SaveB;
  367.     X.Clear();  MuWriteErr("GetD overflow, continuing...");
  368.   }
  369.   else  X.ShiftL(X.C);
  370.   X.MultiSI::GetD( Db);
  371.   MuErrRep = SaveB;
  372. } // --end-- MultiF::GetD
  373.  
  374. // ------- Db = this Mod Base**4, in a range, + or -
  375. void MultiF::GetDIn( double& Db, double Lo, double Hi)
  376. {
  377.   Db = N;  Db = Db + C;
  378.   if (Db > 4.0) {
  379.     if (S)  Db = Lo;
  380.     else  Db = Hi;
  381.     return;
  382.   }
  383.   GetD( Db);
  384.   if (Db < Lo)  Db = Lo;
  385.   if (Db > Hi)  Db = Hi;
  386. } // --end-- MultiF::GrtDIn
  387.  
  388. // ------- this = this + D1
  389. void MultiF::RAdd1( Mu1Digit D1)
  390. // Allows literal D1, D1 is not changed, sets MuErr True if overflow
  391. {
  392.   MultiF X(1);  X.SetTo1( D1);
  393.   RAdd(X);
  394. } // --end-- MultiF::RAdd1
  395.  
  396. // ------- this = X + D1
  397. void MultiF::Add1( MultiF& X, Mu1Digit D1)
  398. // Allows literal D1, D1 is not changed, sets MuErr True if overflow
  399. {
  400.   SetTo(X);
  401.   RAdd1( D1);
  402. } // --end-- MultiF::Add1
  403.  
  404. // ------- Comp = sign( this - X)
  405. int MultiF::Comp( MultiF& X)
  406. // Sign = -1 if this < X; Sign = 0 if this = X; Sign = +1 if this > X
  407. // this.S and X.S used but not changed;  this, X not changed
  408. {
  409.   MultiF Dif(1 + MaxL( N, X.N));
  410.   Dif.Sub( *this, X);
  411.   if      (Dif.ZTest())  return 0;
  412.   else if (Dif.S)        return -1;
  413.   else  return +1;
  414. } // --end-- MultiF::Comp
  415.  
  416. // ------- Output String and line feed every 80, Tot = characters on line so far
  417. void Write80( ostream& Out, char* St, long& Tot)
  418. {
  419.   if (&Out == &cout)  Out << St;
  420.   else
  421.     for ( int i = 0; St[i] != '\0'; i++) {
  422.       Out << St[i];
  423.       Tot++;
  424.       if (Tot % 80 == 0)  Out << nL;
  425.     }
  426. } // --end-- Write80
  427.  
  428. // ------- Output String and line feed every 80, Tot = characters on line so far
  429. void Write80Ln( ostream& Out, char* St, long& Tot)
  430. {
  431.   Write80( Out, St, Tot);  Out << nL;
  432. } // --end-- Write80Ln
  433.  
  434. // ------- Convert +Double to a String integer
  435. void StrDI(char* St, double D);
  436.  
  437. // ------- Convert +Double to a String integer
  438. void StrDI(char* St, double D)
  439.    // Kludge around problem with Windows 3.0 386 Mode, Str( Single, St)
  440.    //  or Write( Single) or Double or Real can hang system }
  441. {
  442.   long I1, I2;
  443.   char St1[ 11];
  444.   char St2[ 11];
  445.   double DI;
  446.  
  447.   D = D + 0.5;
  448.   modf( D / 1.0E9, &DI);  I1 = DI; // integer part
  449.   if (!I1) {
  450.     modf( D, &DI);  I2 = DI; // integer part
  451.     sprintf( St, "%ld", I2);
  452.   }
  453.   else {
  454.     modf(D - 1.0E9 * I1, &DI);  I2 = DI + 1000000000L;
  455.     sprintf( St1, "%ld", I1);
  456.     sprintf( St2, "%ld", I2);
  457.     St = strcat( St1, St2+1);
  458.   }
  459. } // --end--  StrDI
  460.  
  461.  
  462. // -------
  463. void Write1( ostream& Out, long& Tot, int J, int K, long& Did,
  464.              long DidZ, char* St) // Part of WriteFixed
  465. {
  466.   char* StL = " ";
  467.   while (J < K) {
  468.     J++;
  469.     if ((Did != DidZ) && (MuDpG) && !(Did % MuDpG))
  470.        Write80( Out, ",", Tot);
  471.     StL[0] = St[J];
  472.     Write80( Out, StL, Tot);
  473.     Did++;
  474.   }
  475. } // --end--  Write1
  476.  
  477. // -------
  478. void WriteFixed( MultiF& T, ostream& Out, long& Tot)
  479. {
  480.   int J, K;
  481.   long Did, DidZ;
  482.   char St[ 15];
  483.  
  484.   Did = 0;
  485.   if (T.S)  Write80( Out, "-", Tot);
  486.   if (T.N + T.C > 0) {
  487.     sprintf( St, "%.0f", Base + T.V[ T.N - 1]);
  488.     J = 0;
  489.     do
  490.       J++;
  491.     while ((St[J] == '0') && (J != MuDMax));
  492.     Did = J - 1 - MuDMax;
  493.     if (T.N + T.C > 1)  Did -= MuDMax;
  494.     if (MuDpG != 0)     Did  = Did % MuDpG;
  495.     if (Did < 0)        Did += MuDpG;
  496.     DidZ = Did;
  497.     K = MuDMax;  J--;
  498.     Write1( Out, Tot, J, K, Did, DidZ, St);
  499.     if ((T.N + T.C) > 1) {
  500.       if (T.N > 1)
  501.     sprintf( St, "%.0f", Base + T.V[ T.N - 2]);
  502.       else
  503.     sprintf( St, "%.0f", Base);
  504.       K = MuDMax;  J = 0;
  505.       Write1( Out, Tot, J, K, Did, DidZ, St);
  506.     }
  507.   }
  508.   else
  509.     Write80( Out, "0", Tot);
  510.   Write80( Out, ".", Tot);
  511.   DidZ = Did;
  512.   if (-T.C == 2) {
  513.     sprintf( St, "%.0f", Base + T.V[1]);
  514.     K = MuDMax;  J = 0;
  515.     Write1( Out, Tot, J, K, Did, DidZ, St);
  516.   }
  517.   if (-T.C > 0) {
  518.     sprintf( St, "%.0f", Base + T.V[0]);
  519.     J = MuDMax +1;
  520.     do
  521.       J--;
  522.     while (St[J] == '0');
  523.     K = J;  J = 0;
  524.     Write1( Out, Tot, J, K, Did, DidZ, St);
  525.   }
  526.   else
  527.     Write80( Out, "0", Tot);
  528. } // --end--  WriteFixed
  529.  
  530. // ------- Output this as a line of text, Tot = characters on line so far
  531. void MultiF::Writ( ostream& Out, long& Tot)
  532. // this is not modified
  533. {
  534.   MultiF T; // Temp on heap
  535.   long   LZ;
  536.   char*  Ch = " ";
  537.   int    J, K;
  538.   long   I, Did;
  539.   char   St[ 15];
  540.   char   St2[ 15];
  541.   char*  StL = " ";
  542.   char*  StP = "0.";
  543.   char*  Sep = ",";
  544.   double DE;
  545.   Bool   NewT;
  546.   Bool   Fixed;
  547.  
  548.   if (ZTest())  { Write80( Out, "0.0 (False)", Tot);  return; }
  549.   if (!C && !S && EQ1(1)) {
  550.     Write80( Out, "1.0 (True) ", Tot);   return;
  551.   }
  552.   Fixed = (!ScieN) && (N + C <= 2) && (-C <= 2);
  553.   if (Fixed) {
  554.     WriteFixed( *this, Out, Tot);  // format <= 99999999999999.99999999999999
  555.     return;
  556.   }
  557.   T = *this;
  558.   NewT = False;
  559.   if ((FTn > 0) && (FTn < FMC) && (N > FMC - FTn)) {
  560.     MultiF TL( FMC - FTn);  NewT = True;
  561.     TL.SetTo( *this);  T = TL;  TL.V = NULL;
  562.   }
  563.   if (T.S)  Write80( Out, "-", Tot);
  564. //Str( Base + T.V[ T.N - 1] :0:0, St);
  565.   sprintf( St, "%.0f", Base + T.V[ T.N - 1]);
  566.   J = 0;
  567.   do
  568.     J++;
  569.   while ((St[J] == '0') && (J != MuDMax));
  570.   if (MuDOnly) {
  571.     Sep = " ";  StP = "0 ";
  572.     if (!MuDpG)  StP = "0";
  573.   }
  574.   StP[0] = St[J];
  575.   Write80( Out, StP, Tot);
  576.   LZ = J - 1;
  577.   Did = 0;
  578.   K = MuDMax;
  579.   if (T.N == 1)  while (St[K] == '0')  K--;
  580.   if ((J >= K) && (T.N == 1 ))  Write80( Out, "0", Tot);
  581.   while (J < K) {
  582.     J++;
  583.     if (Did && MuDpG && !(Did % MuDpG))
  584.        Write80( Out, Sep, Tot);
  585.     StL[0] = St[J];
  586.     Write80( Out, StL, Tot);
  587.     Did++;
  588.     if (Did >= FDn-1)  K = 0;
  589.   }
  590.   if (Did < FDn-1)  I = 2;  else  I = T.N + 1;
  591.   while (I <= T.N) {
  592. //  Str( Base + T.V[ T.N - I] :0:0, St);
  593.     sprintf( St, "%.0f", Base + T.V[ T.N - I]);
  594.     K = MuDMax;
  595.     if (I == T.N)  while (St[K] == '0')  K--;
  596.     J = 1;
  597.     while (J <= K) {
  598.       if (Did && MuDpG && !(Did % MuDpG))
  599.         Write80( Out, Sep, Tot);
  600.       StL[0] = St[J];
  601.       Write80( Out, StL, Tot);
  602.       Did++;
  603.       if (Did >= FDn-1) { K = 0;  I = T.N; }
  604.       J++;
  605.     }
  606.     MuInterrupt;
  607.     if (MuAbort)  I = T.N;
  608.     I++;
  609.   }
  610.   DE = T.C;
  611.   DE = MuDMax * (DE + T.N) - LZ - 1.0;
  612.   if (DE >= 0)  Ch = "+";
  613.   else  Ch = "-";
  614.   if (!MuDOnly) {
  615.     Write80( Out, " E", Tot);
  616.     Write80( Out,   Ch, Tot);
  617. //Str( Abs( DE) :0:0, St);
  618. //StrDi( St, DE);
  619.     sprintf( St, "%.0f", fabs( DE));
  620.     K = strlen( St);
  621.     if ((K > MuDpG) && MuDpG) {
  622.       J = K + 1;
  623.       for (I = 1; I <= ((K-1) / MuDpG); I++) {
  624.         J -= MuDpG;
  625.         strcpy( St2, St+J-1);  St[J-1] = ',';  strcpy(St+J, St2);
  626.       }
  627.     }
  628.     Write80( Out, St, Tot);
  629.     if (Did >= MuLenD) {
  630.       sprintf( St, "%ld", Did + 1);
  631.       Write80( Out, " (", Tot);
  632.       Write80( Out,   St, Tot);
  633.       Write80( Out,  ")", Tot);
  634.     }
  635.     sprintf( St, "%ld", N * (long)( MuDMax));
  636.     Write80( Out, " [", Tot);
  637.     Write80( Out,   St, Tot);
  638.     Write80( Out,  "]", Tot);
  639.   }
  640.   if (!NewT) T.V = NULL;
  641. } // --end-- MultiF::Writ
  642.  
  643. // ------- Output this as a line of text, Tot = characters on line so far
  644. void MultiF::WritLn( ostream& Out, long& Tot)
  645. {
  646.   Writ( Out, Tot);  Out << nL;
  647. } // --end-- MultiF::WritLn
  648.  
  649. // ------- this = this * D1
  650. void MultiF::RMul1( Mu1Digit D1)
  651. // Multi-precision one super digit multiply; D1 is not changed
  652. // MuErr set true if this overflows
  653. {
  654.   MultiSI::RMul1( D1);
  655.   Norm();
  656. } // --end-- MultiF::RMul1
  657.  
  658. // ------- this = X * D1
  659. void MultiF::Mul1( MultiF& X, Mu1Digit D1)
  660. // Multi-precision one super digit multiply; X and D1 are not changed
  661. // MuErr set true if this overflows
  662. {
  663.   SetTo(X);
  664.   RMul1( D1);
  665. } // --end-- MultiF::Mul1
  666.  
  667. // ------- this = X * Y
  668. void MultiF::MulSlow( MultiF& X, MultiF& Y)
  669. // X or Y or both may have the same address in memory as this
  670. // X, Y are not changed unless they share memory with this
  671. // Sets MuErr True if overflow
  672. {
  673.   long     MM;
  674.   long     I, J, II, IJ;
  675.   long     Bias;
  676.   Mu1Digit K;
  677.   Mu2Digit T;
  678.   double   D;
  679.  
  680.   MultiF XL = X; // Local pointer to X (or Y) so XL.N >= YL.N
  681.   MultiF YL = Y; // Local pointer to Y (or X)
  682.   if (XL.N < YL.N) {
  683.     XL = Y;  YL = X;
  684.   }
  685.   I = 1;
  686.   if (!Expanded)  I = 4; // Four guard "digits", no sticky bit
  687.   MM = MinL( XL.N + YL.N, M + I);
  688.   MM = MinL( MM, FMC + I);
  689.   MM = MinL( MM, MuNMax);
  690.   MultiF Z(MM); // Temp on heap
  691.   Z.S = (XL.S && !YL.S) || (!XL.S && YL.S); // XOR
  692.   Bias = XL.N + YL.N - Z.M;  if (Bias < 0)  Bias = 0;
  693.   for (I = 0; I < (XL.N - Bias); I++)  Z.V[I] = 0; // Clear lower of Z
  694.   KeyHit = False;
  695.   for (J = 0; J < YL.N; J++) {
  696.     II = Bias - J;  if (II < 0)  II = 0;
  697.     IJ = II + J - Bias;  K = 0;
  698.     for (I = II; I < XL.N; I++) {  // Mult, add and carry next Digit
  699.       T = XL.V[I] * YL.V[J] + Z.V[ IJ] + K;
  700.       modf( (double)(T / Base), &D);  K = D;  // integer part
  701.       Z.V[ IJ] = T - K * Base;
  702.       IJ++;
  703.     }
  704.     Z.V[ IJ] = K;
  705.     MuInterrupt;
  706.     if (MuAbort) {
  707.       Clear();  XL.V = YL.V = NULL;  return;
  708.     }
  709.     if (KeyHit) {
  710.       if (Echo)  LogF << "Floating Multiply: " 
  711.                       << (int)((100.0 * (J+1)) / YL.N) << " Percent done" << nL;
  712.       cout << "Floating Multiply: " 
  713.            << (int)((100.0 * (J+1)) / YL.N) << " Percent done" << nL;
  714.       KeyHit = False;
  715.     }
  716.   }
  717.   Z.N = IJ;
  718.   if (K)  Z.N++;
  719.   D = XL.C;  D = D + YL.C;  D = D + Bias;
  720.   if (D > LONG_MAX) {
  721.     Z.Clear();
  722.     MuWriteErr("Floating multiply overflow, continuing...");
  723.   }
  724.   else
  725.     if (-D > LONG_MAX) {
  726.       Z.Clear();
  727.     }
  728.     else
  729.       Z.C = D;
  730.   SetTo(Z);
  731.   Norm();
  732.   XL.V = YL.V = NULL;
  733. } // --end-- MultiFSlow::Mul
  734.  
  735. // ------- this = X * Y (Fast)
  736. void MultiF::MulFast( MultiF& X, MultiF& Y)
  737. // X or Y or both may have the same address in memory as this
  738. // X, Y are not changed unless they share memory with this
  739. // Sets MuErr True if overflow
  740. {
  741. //  MulSlow( X, Y); // Stub for now
  742.   double CD = (double)X.C + Y.C;
  743.   MultiF Z(X.N + Y.N); // Temp on heap
  744.   Z.MultiSI::Mul(X, Y);  
  745.   Z.Norm();
  746.   SetTo(Z);
  747.   CD += C;
  748.   if (CD > LONG_MAX) {
  749.     Clear();  MuWriteErr("Floating multiply overflow, continuing...");  return;
  750.   }
  751.   if (-CD > LONG_MAX) {
  752.     Clear();  return;  
  753.   }
  754.   C = CD;
  755. } // --end-- MultiF::MulFast
  756.  
  757. // ------- this = X * Y (Fast or Slow)
  758. void MultiF::Mul( MultiF& X, MultiF& Y)
  759. // X or Y or both may have the same address in memory as this
  760. // X, Y are not changed unless they share memory with this
  761. // Sets MuErr True if overflow
  762. {
  763.   float prod = (float) X.N * Y.N;
  764.   if (prod < 20.0)
  765.     MulSlow(X, Y);
  766.   else
  767.     MulFast(X, Y);
  768. } // --end-- MultiF::Mul
  769.  
  770. // ------- this = this * X
  771. void MultiF::RMul( MultiF& X)
  772. // X may have the same address in memory as this
  773. // X is not changed unless it shares memory with this
  774. // Sets MuErr True if overflow
  775. {
  776.   Mul( *this, X);
  777. } // --end-- MultiF::RMul
  778.  
  779. // ------- this = X * X
  780. void MultiF::Sq( MultiF& X)
  781. // X may have the same address in memory as this
  782. // X is not changed unless it shares memory with this
  783. // Sets MuErr True if overflow
  784. {
  785.   Mul(X, X);
  786. } // --end-- MultiF::Sq
  787.  
  788. // ------- this = this * this
  789. void MultiF::RSq()
  790. // Sets MuErr True if overflow
  791. {
  792.   Mul( *this, *this);
  793. } // --end-- MultiF::RSq
  794.  
  795. // ------- this = B ** P
  796. void MultiF::Pow( MultiF& B, MultiF& P)
  797. // B or P or both may have the same address in memory as this
  798. // B, P are not changed unless they share memory with this
  799. // Sets MuErr True if overflow or error
  800. {
  801.   double D, DI;
  802.  
  803.   MultiF PL( P.N + P.C);  PL.SetTo(P);
  804.   PL.ShiftL( PL.C);
  805.   MultiF Z(M);
  806.   Z.PowMB(B, PL);
  807.   Z.Norm();
  808.   PL.GetD(D);
  809.   modf( D, &DI); // Integer part
  810.   D = DI * B.C + Z.C;
  811.   if (D > LONG_MAX) {
  812.     Z.Clear();  MuWriteErr("B ** P overflow, continuing...");
  813.   }
  814.   else {
  815.     if (-D > LONG_MAX) {
  816.       Z.Clear();
  817.     }
  818.     else  Z.C = D;
  819.   }
  820.   SetTo(Z);
  821. } // --end-- MultiF::Pow
  822.  
  823. // ------- this = this ** P
  824. void MultiF::RPow( MultiF& P)
  825. // P may have the same address in memory as this
  826. // P is not changed unless it shares memory with this
  827. // Sets MuErr True if overflow or error
  828. {
  829.   Pow( *this, P);
  830. } // --end-- MultiF::RPow
  831.  
  832. // ------- this = this / D1
  833. void MultiF::RDiv1( Mu1Digit D1)
  834. // Multi-precision one super digit divide, D1 is not changed
  835. // Allows literal D1, MuErr set True if D1 = 0
  836. {
  837.   MultiF X(1);  X.SetTo1( D1);  // Temp on heap
  838.   RDiv(X);
  839. } // --end-- MultiF::RDiv1
  840.  
  841. // ------- this = U / D1
  842. void MultiF::Div1( MultiF& U, Mu1Digit D1)
  843. // Multi-precision one super digit divide, D1 and U not changed
  844. // U and this may have same address in memory, then U is changed
  845. // Allows literal D1, MuErr set True if D1 = 0
  846. {
  847.   SetTo( U);
  848.   RDiv1( D1);
  849. } // --end-- MultiF::Div1
  850.  
  851. // ------- Adjust .C, part of MultiF::Divi
  852. void AdjustC( MultiF& Z, MultiF& D)
  853. {
  854.   if ((Z.C > 0) && (D.C < 0))
  855.     if (Z.C > (LONG_MAX + D.C)) {
  856.       Z.Clear();  MuWriteErr("Floating divide overflow, continuing...");
  857.     }
  858.     else
  859.       Z.C -= D.C;
  860.   else
  861.     if ((Z.C < 0) && (D.C > 0)) 
  862.       if (-Z.C > (LONG_MAX - D.C)) {
  863.     Z.Clear();
  864.       }
  865.       else
  866.     Z.C -= D.C;
  867.     else
  868.       if (!Z.ZTest())  Z.C -= D.C;
  869. } // --end-- AdjustC
  870.  
  871. // ------- this = U / D
  872. void MultiF::Divi( MultiF& U, MultiF& D)
  873. // U or D or both may have the same address in memory as this
  874. // U, D are not changed unless they share memory with this
  875. // Sets MuErr True if D = 0
  876. {
  877.   long   K, ML;
  878.   MultiF Z;  // Temp on heap
  879.  
  880.   K = MinL(M, FMC);
  881.   if (!Expanded)  K += 4; // Four guard "digits", no sticky bit
  882.   if (K > MuNMax)  K = MuNMax;
  883.  
  884.   if ((K < 13) || (D.N < (K/3))) {
  885.     Z.NMax(K + D.N);  Z.SetTo(U);
  886.     Z.ShiftL(K - U.N + D.N );
  887.     Z.MultiSI::RDiv(D);
  888.     Z.Norm();
  889.     AdjustC(Z, D);
  890.     SetTo(Z);
  891.   } else {
  892.     long KL = 1 + K / 2;
  893.     long DN = MinL( KL, D.N);
  894.     if ((V != U.V) && (V != D.V)) { ML = M;  NMax(0); }
  895.     Z.NMax(KL + DN);  Z.SetTo1(1);
  896.     Z.ShiftL(KL - 1 + DN );
  897.     long DV = D.N - DN;
  898.     D.V += DV;  D.N -= DV;
  899.     Z.MultiSI::RDiv(D);
  900.     D.V -= DV;  D.N += DV;
  901.     Z.C = -(KL - 1 + DN + DV);
  902.     Z.Norm();
  903.     AdjustC(Z, D);  // Z ~= 1 / D
  904.     Z.NMax( Z.N);
  905.  
  906.     MultiF ZZ(K);  // ZZ = Z + Z - (Z*Z)*D
  907.     ZZ.Sq(Z);
  908.     ZZ.RMul(D);
  909.     ZZ.ChS();
  910.     ZZ.RAdd(Z);
  911.     ZZ.RAdd(Z);
  912.     Z.NMax(0);
  913.     ZZ.RMul(U);
  914.     if ((V != U.V) && (V != D.V))  NMax( ML);
  915.     SetTo(ZZ);
  916.   }
  917. } // --end-- MultiF::Divi
  918.  
  919. // ------- this = this / D
  920. void MultiF::RDiv( MultiF& D)
  921. // D may have the same address in memory as this
  922. // D is not changed unless it shares memory with this
  923. // Sets MuErr True if D = 0
  924. {
  925.   Divi( *this, D);
  926. } // --end-- MultiF::RDiv
  927.  
  928. // ------- this = X ! = 1*2*3*...*X
  929. void MultiF::Fac( MultiF& X)
  930. // X may have the same address in memory as this
  931. // X is not changed unless it shares memory with this
  932. // Sets MuErr True if overflow or error
  933. {
  934.   double NN;
  935.  
  936.   if (X.S) {
  937.     MuWriteErr("Cannot take factorial of number < zero");
  938.     Clear();  return;
  939.   }
  940.   if ((X.N + X.C) <= 2)  X.GetD( NN);
  941.   else  NN = 1.0 + LONG_MAX;
  942.   if (NN > LONG_MAX) {
  943.     MuWriteErr("Number too large for factorial function");
  944.     Clear();  return;
  945.   }
  946.   MultiF XL(2);
  947.   SetTo1(1);
  948.   TracN = 0;
  949.   while (NN > 1) {
  950.     XL.SetToD( NN);
  951.     Mul( XL, *this);
  952.     if (MuAbort) {
  953.       if (Echo)  LogF << "MultiF::Fac abort, N = " << NN << nL;
  954.       cout << "MultiF::Fac abort, N = " << NN << nL;
  955.       Clear();  NN = 0;
  956.     }
  957.     NN = NN - 1;
  958.     Trace = NN;
  959.     TracN++;
  960.   }
  961.   Trace = 0;
  962. } // --end-- MultiF::Fac
  963.  
  964. // ------- this = this ! = 1*2*3*...*X
  965. void MultiF::RFac()
  966. // Sets MuErr True if overflow or error
  967. {
  968.   Fac( *this);
  969. } // --end-- MultiF::RFac
  970.  
  971. // ------- this = SqRt(X)
  972. void MultiF::SqRoot( MultiF& X)
  973. // X is not changed, it is ok for X and this to have the same location in memory
  974. // but then X will be changed.  Sets MuErr True if error
  975. {
  976.   long Sh, ZC, ML, SaveFMC;
  977.  
  978.   if (X.ZTest()) { Clear();  return; }
  979.   Sh = MinL(M, FMC) + 1;
  980.   if (!Expanded)  Sh += 4; // Four guard "digits", no sticky bit
  981.   if (Sh > MuNMax)  Sh = MuNMax;
  982.   if ((Sh - X.N - X.C) & 1)  Sh--; // if odd
  983.   if (V != X.V) { ML = M;  NMax(0); }
  984.   MultiF Z( Sh);  Z.SetTo(X);
  985.   Z.ShiftL( Sh - Z.N);
  986.   Z.C = Z.C / 2;  ZC = Z.C;
  987.   Z.MultiSI::RSqRt();  Z.C = ZC;
  988.   Z.Norm();
  989.   if (Z.ZTest()) {
  990.     if (V != X.V) { NMax( ML); }
  991.     Clear();
  992.   }
  993.   else {
  994.     Z.NMax( Z.N); 
  995.     MultiF Y( Sh);
  996.     SaveFMC = FMC;
  997.     FMC = Sh;
  998.     Y.Divi(X, Z);
  999.     Y.RSub(Z);
  1000.     Y.RDiv1(2);
  1001.     Y.RAdd(Z); // (X/Z - Z) / 2  + Z
  1002.     Z.NMax( 0);
  1003.     if (V != X.V)  NMax( ML);
  1004.     SetTo(Y);
  1005.     FMC = SaveFMC;
  1006.     Norm();
  1007.   }
  1008. } // --end-- MultiF::SqRoot
  1009.  
  1010. // ------- this = SqRt( this)
  1011. void MultiF::RSqRt()
  1012. // Sets MuErr True if error
  1013. {
  1014.   SqRoot( *this);
  1015. } // --end-- MultiF::RSqRt
  1016.  
  1017. // ------- this = Integer part of X
  1018. void MultiF::Inte( MultiF& X)
  1019. {
  1020.   long I, J;
  1021.   Mu1Digit SvFRound; // Save Floating point rounding constant, Base / 2 or 0
  1022.  
  1023.   SvFRound = FRound;
  1024.   FRound = 0;
  1025.   SetTo(X);
  1026.   if (-C >= N) {
  1027.     FRound = SvFRound;
  1028.     Clear();  return;
  1029.   }
  1030.   J = -C - 1;
  1031.   for (I = 0; I <= J; I++)  V[I] = 0;
  1032.   Norm();
  1033.   FRound = SvFRound;
  1034. } // --end-- MultiF::Inte
  1035.  
  1036. // ------- this = Integer part of this
  1037. void MultiF::RInt()
  1038. {
  1039.   Inte( *this);
  1040. } // --end-- MultiF::RInt
  1041.  
  1042. // ------- this = Fractional part of X
  1043. void MultiF::Fract( MultiF& X)
  1044. {
  1045.   long I, J;
  1046.   Mu1Digit SvFRound; // Save Floating point rounding constant, Base / 2 or 0
  1047.  
  1048.   SvFRound = FRound;
  1049.   FRound = 0;
  1050.   SetTo(X);
  1051.   J = -C;
  1052.   if (J < 0)  J = 0;
  1053.   for (I = J; I < N; I++)  V[I] = 0;
  1054.   Norm();
  1055.   FRound = SvFRound;
  1056. } // --end-- MultiF::Fract
  1057.  
  1058. // ------- this = Fractional part of this
  1059. void MultiF::RFrac()
  1060. {
  1061.   Fract( *this);
  1062. } // --end-- MultiF::RFrac
  1063.  
  1064. // ------- this = X with one less super digit
  1065. void MultiF::Lop( MultiF& X)
  1066. {
  1067.   long L;
  1068.   Bool RoundIt;
  1069.  
  1070.   SetTo(X);
  1071.   if ((FRound > 0) && (V[0] >= FRound)) {
  1072.     RoundIt = True;
  1073.     if (V[0] == FRound) {
  1074.       L = V[1];
  1075.       if (!(L & 1))  RoundIt = False; // if not odd
  1076.     }
  1077.     if (RoundIt) {
  1078.       MultiF T(2);
  1079.       T.V[0] = 0;
  1080.       T.V[1] = 1;
  1081.       T.N = 2;
  1082.       T.S = S;
  1083.       T.C = C;
  1084.       V[0] = 0;
  1085.       RAdd(T);
  1086.       return;
  1087.     }
  1088.   }
  1089.   V[0] = 0;
  1090.   Norm();
  1091. } // --end-- MultiF::Lop
  1092.  
  1093. // ------- this = this with one less least significant "digit"
  1094. void MultiF::RLop()
  1095. {
  1096.   Lop( *this);
  1097. } // --end-- MultiF::RLop
  1098.  
  1099. // --end-- MultiF's methods implementation
  1100.  
  1101. // --end-- Method implementation
  1102.  
  1103. // ------- Other services
  1104.  
  1105. // ------- Returns max of X, Y
  1106. int Max(int X, int Y)
  1107. {
  1108.   if (X > Y)  return X;
  1109.   else  return Y;
  1110. } // --end-- Max
  1111.  
  1112. // ------- Returns min of X, Y
  1113. int Min(int X, int Y)
  1114. {
  1115.   if (X < Y)  return X;
  1116.   else  return Y;
  1117. } // --end-- Min
  1118.  
  1119. // ------- Returns max of X, Y
  1120. long MaxL(long X, long Y)
  1121. {
  1122.   if (X > Y)  return X;
  1123.   else  return Y;
  1124. } // --end-- MaxL
  1125.  
  1126. // ------- Returns min of X, Y
  1127. long MinL(long X, long Y)
  1128. {
  1129.   if (X < Y)  return X;
  1130.   else  return Y;
  1131. } // --end-- MinL
  1132.  
  1133. // ------- Returns max of X, Y
  1134. double MaxD( double X, double Y) 
  1135. {
  1136.   if (X > Y)  return X;
  1137.   else  return Y;
  1138. } // --end-- MaxD
  1139.  
  1140. // ------- Returns min of X, Y
  1141. double MinD( double X, double Y)
  1142. {
  1143.   if (X < Y)  return X;
  1144.   else  return Y;
  1145. } // --end-- MinD
  1146.  
  1147. // ------- Compute X = Pi by algorithm a
  1148. void PiAA( MultiF& X, Bool WReset, Bool ReStart)
  1149. // Uses 3 large "Registers" X, Y, and Z.  L is small.
  1150. // if WReset then Write reset data after each iteration
  1151. // if ReStart then Read reset data before starting
  1152. {
  1153.   int J;
  1154.   int It; // Number of iterations(...
  1155.   char* ResetSt = "PiReset.Txt";
  1156.   char* OldReSt = "PiReset.Old";
  1157.   ofstream WritF;  // Reset text file, output
  1158.   FILE*    ReadF; // Reset text file, input
  1159.   char Ch;
  1160.   char Name[21];
  1161.   Bool OK;
  1162.  
  1163.   MultiF Z( FMC); // Temp
  1164.   MultiF Y( FMC); // y0, y1, ...
  1165.   MultiF L( 2);   // -(2**N)
  1166.   It = log( 1024.0 * (FMC * (long)( MuDMax) + 3) / 2788.0 ) / log( 2.0) + 2;
  1167.   cout << "FMC * MuDMax = " << FMC * MuDMax << "  It = " << It << nL;
  1168.   if (!ReStart) {
  1169.     Diag("Starting Pi");
  1170.     if (Echo)  LogF << "  Starting to compute Pi by algorithm a" << nL;
  1171.     cout << "  Starting to compute Pi by algorithm a" << nL;
  1172.     int i;
  1173.     X.Value("0.5", i); // X0 = 1/2 = 0.5
  1174.     Diag("Start Y.SqRoot(0.5);");
  1175.     Y.SqRoot(X);       // Y0 = 1/SqRt(2) = SqRt(0.5)
  1176.     Diag("End Y.SqRoot(0.5);");
  1177.     L.SetTo1(-1);      // L  = -1
  1178.     J = 1;
  1179.   } else {
  1180.     Diag("Restarting Pi");
  1181.     do {
  1182.       if ((ReadF = fopen( ResetSt, "rt")) == NULL) { // Open Reset file for read
  1183.         cout << "Cannot open " << ResetSt
  1184.              << ", Type a character, Enter to retry" << nL;
  1185.         cin >> Ch;
  1186.       }
  1187.     } while (ReadF == NULL);
  1188.     fscanf( ReadF, "%d", &J);
  1189.     L.MultiSI::ReadSI( ReadF, Name, OK);
  1190.     Y.MultiSI::ReadSI( ReadF, Name, OK);
  1191.     fscanf( ReadF, " (%*ld) %ld", &Y.C);
  1192.     X.MultiSI::ReadSI( ReadF, Name, OK);
  1193.     fscanf( ReadF, " (%*ld) %ld", &X.C);
  1194.     fclose( ReadF);
  1195.     J++;
  1196.   }
  1197.   for ( ; J <= It; J++) {
  1198.     Diag("Starting Iteration");
  1199.     if (Echo) LogF << "  Start of iteration number " << J << " of " << It << nL;
  1200.     cout << "  Start of iteration number " << J << " of " << It << nL;
  1201.     // Y = (1 - SqRoot(1 - Y**2)) / (1 + SqRoot(1 - Y**2))
  1202.     Diag("Start Y.RSq();");
  1203.     Y.RSq();
  1204.     Diag("End Y.RSq();");
  1205.     Y.ChS();
  1206.     Y.RAdd1(1);
  1207.     Diag("Start Y.RSqRt();");
  1208.     Y.RSqRt();
  1209.     Diag("End Y.RSqRt();");
  1210.     Z.Add1(Y, 1);          // Z = 1 + SqRoot(1 - Y**2)
  1211.     Y.ChS();
  1212.     Y.RAdd1(1);            // = 1 - SqRoot(1 - Y**2)
  1213.     Diag("Start Y.RDiv();");
  1214.     Y.RDiv(Z);
  1215.     Diag("End Y.RDiv();");
  1216.  
  1217.     // X = (1 + Y)**2 * X - 2**N * Y
  1218.     L.RAdd(L);            // L = -2**1, -2**2, ...
  1219.     Z.Add1(Y, 1);         // Z = (1 + Y) 
  1220.     Diag("Start Z.RSq();");
  1221.     Z.RSq();
  1222.     Diag("Start X.RMul(Z);");
  1223.     X.RMul(Z);            // X = (1 + Y)**2 * X
  1224.     Diag("Start Z.Mul(Y, L);");
  1225.     Z.Mul(Y, L);          // Z = - 2**N * Y
  1226.     Diag("Start X.RAdd(Z);");
  1227.     X.RAdd(Z);
  1228.     Diag("End X.RAdd(Z);");
  1229.  
  1230.     if (MuAbort) {
  1231.       Diag("Operator abort");  break;
  1232.     }
  1233.     if (WReset && !MuAbort) {
  1234.       unlink( OldReSt);         // Delete backup file
  1235.       rename(ResetSt, OldReSt); // Rename the Reset file
  1236.       do {
  1237.     WritF.open( ResetSt, ios::trunc); // Open Reset file for rewrite
  1238.     if (WritF.fail()) {
  1239.       cout << "Cannot open " << ResetSt
  1240.                << ", Type a character, Enter to retry" << nL;
  1241.            cin >> Ch;
  1242.         }
  1243.       } while (WritF.fail());
  1244.       WritF << nL << J << " : J" << dL << "L =" << nL;
  1245.       MuTot = 0; L.MultiSI::Writ( WritF); WritF << dL << "Y =" << nL;
  1246.       MuTot = 0; Y.MultiSI::Writ( WritF); WritF << dL << Y.C << dL
  1247.     << "X =" << nL;
  1248.       MuTot = 0; X.MultiSI::Writ( WritF); WritF << dL << X.C << dL;
  1249.       WritF.close();
  1250.     }
  1251.   }
  1252.   if (!MuAbort) {
  1253.     Diag("All done but...");
  1254.     if (Echo)  LogF << "  All done but inverting alpha" << nL;
  1255.     cout << "  All done but inverting alpha" << nL;
  1256.     Z.SetTo1(1);         // PI ~= 1 / X 
  1257.     X.Divi(Z, X);
  1258.     Diag("Pi is done");
  1259.     if (Echo)  LogF << "  All done computing Pi" << dL;
  1260.     cout << "  All done computing Pi" << dL;
  1261.   }
  1262. } // --end-- PiAA
  1263.  
  1264. // ------- Compute X = Pi by algorithm b
  1265. void PiAB( MultiF& X, Bool WReset, Bool ReStart)
  1266. // Uses 4 large "Registers" X, Y, Z, and T.  L is small.
  1267. // if WReset then Write reset data after each iteration
  1268. // if ReStart then Read reset data before starting
  1269. {
  1270.   int J;
  1271.   int It; // Number of iterations(...
  1272.   char* ResetSt = "PiReset.Txt";
  1273.   char* OldReSt = "PiReset.Old";
  1274.   ofstream WritF;  // Reset text file, output
  1275.   FILE*    ReadF; // Reset text file, input
  1276.   char Ch;
  1277.   char Name[21];
  1278.   Bool OK;
  1279.  
  1280.   MultiF T( FMC); // Temp
  1281.   MultiF Z( FMC); // Temp
  1282.   MultiF Y( FMC); // y0, y1, ...
  1283.   MultiF L( 2);   // -(2**(2 * N + 1))
  1284.   It = log( 1024.0 * (FMC * (long)( MuDMax) + 3) / 2788.0 ) / log( 4.0) + 1;
  1285.   cout << "FMC * MuDMax = " << FMC * MuDMax << "  It = " << It << nL;
  1286.   if (!ReStart) {
  1287.     Diag("Starting Pi");
  1288.     if (Echo)  LogF << "  Starting to compute Pi by algorithm b" << nL;
  1289.     cout << "  Starting to compute Pi by algorithm b" << nL;
  1290.     Y.SetTo1(2); // Y0 = SqRt(2) - 1
  1291.     Y.RSqRt();
  1292.     X.SetTo(Y);
  1293.     Y.RAdd1(-1);
  1294.     X.RMul1(-4); // X0 = 6 - 4 * SqRt(2);
  1295.     X.RAdd1(6);
  1296.     L.SetTo1(-2);
  1297.     J = 1;
  1298.   } else {
  1299.     Diag("Restarting Pi");
  1300.     do {
  1301.       if ((ReadF = fopen( ResetSt, "rt")) == NULL) { // Open Reset file for read
  1302.         cout << "Cannot open " << ResetSt
  1303.              << ", Type a character, Enter to retry" << nL;
  1304.         cin >> Ch;
  1305.       }
  1306.     } while (ReadF == NULL);
  1307.     fscanf( ReadF, "%d", &J);
  1308.     L.MultiSI::ReadSI( ReadF, Name, OK);
  1309.     Y.MultiSI::ReadSI( ReadF, Name, OK);
  1310.     fscanf( ReadF, " (%*ld) %ld", &Y.C);
  1311.     X.MultiSI::ReadSI( ReadF, Name, OK);
  1312.     fscanf( ReadF, " (%*ld) %ld", &X.C);
  1313.     fclose( ReadF);
  1314.     J++;
  1315.   }
  1316.   for ( ; J <= It; J++) {
  1317.     Diag("Starting Iteration");
  1318.     if (Echo) LogF << "  Start of iteration number " << J << " of " << It << nL;
  1319.     cout << "  Start of iteration number " << J << " of " << It << nL;
  1320.     // Y = (1 - 4thRoot(1 - Y**4)) / (1 + 4thRoot(1 - Y**4))
  1321.     Y.RSq();  Y.RSq();
  1322.     Y.ChS();
  1323.     Y.RAdd1(1);
  1324.     Y.RSqRt();  Y.RSqRt();
  1325.     Z.Add1(Y, 1);          // Z = 1 + 4thRoot(1 - Y**4 
  1326.     Y.ChS();
  1327.     Y.RAdd1(1);            // = 1 - 4thRoot(1 - Y**4
  1328.     Y.RDiv(Z);
  1329.  
  1330.     // X = (1 + Y)**4 * X - 2**(2*N + 1) * Y * (1 + (1 + Y)*Y) 
  1331.     L.RMul1(4);            // L = -2**3, -2**5, ...
  1332.     Z.Add1(Y, 1);          // Z = (1 + Y) 
  1333.     T.Mul(Z, Y);
  1334.     T.RAdd1(1);
  1335.     T.RMul(Y);
  1336.     T.RMul(L);             // T = - 2**(2*N + 1) * Y * (1 + (1 + Y)*Y 
  1337.     Z.RSq();  Z.RSq();
  1338.     X.RMul(Z);             // = (1 + Y)**4 * X 
  1339.     X.RAdd(T);
  1340.  
  1341.     if (MuAbort) {
  1342.       Diag("Operator abort");  break;
  1343.     }
  1344.     if (WReset && !MuAbort) {
  1345.       unlink( OldReSt);         // Delete backup file
  1346.       rename(ResetSt, OldReSt); // Rename the Reset file
  1347.       do {
  1348.     WritF.open( ResetSt, ios::trunc); // Open Reset file for rewrite
  1349.     if (WritF.fail()) {
  1350.       cout << "Cannot open " << ResetSt
  1351.                << ", Type a character, Enter to retry" << nL;
  1352.            cin >> Ch;
  1353.         }
  1354.       } while (WritF.fail());
  1355.       WritF << nL << J << " : J" << dL << "L =" << nL;
  1356.       MuTot = 0; L.MultiSI::Writ( WritF); WritF << dL << "Y =" << nL;
  1357.       MuTot = 0; Y.MultiSI::Writ( WritF); WritF << dL << Y.C << dL
  1358.     << "X =" << nL;
  1359.       MuTot = 0; X.MultiSI::Writ( WritF); WritF << dL << X.C << dL;
  1360.       WritF.close();
  1361.     }
  1362.   }
  1363.   if (!MuAbort) {
  1364.     Diag("All done but...");
  1365.     if (Echo)  LogF << "  All done but inverting alpha" << nL;
  1366.     cout << "  All done but inverting alpha" << nL;
  1367.     T.SetTo1(1);         // PI ~= 1 / X 
  1368.     X.Divi(T, X);
  1369.     Diag("Pi is done");
  1370.     if (Echo)  LogF << "  All done computing Pi" << dL;
  1371.     cout << "  All done computing Pi" << dL;
  1372.   }
  1373. } // --end-- PiAB
  1374.  
  1375. // ------- Compute X = Pi by algorithm b, set Pi to its computed value
  1376. void PiTo( MultiF& X)
  1377. {
  1378.   Bool WReset = False;
  1379.   Bool ReStart = False;
  1380.   PiAB( X, WReset, ReStart); // Compute X = Pi by algorithm b
  1381. } // --end-- PiTo
  1382.  
  1383. // ------- End of other services
  1384.  
  1385. // ------- Initialization section
  1386.  
  1387. // ------- Initialize Multi-precision floating decimal package
  1388. void InitMultiF()
  1389. {
  1390.   InitMulti(); // Initialize Multi-precision package
  1391.   FRound = Base / 2; // Set floating point rounding constant to round}
  1392.   FMC = 2 + 23 / MuDMax; // Set for at least 25 decimal digits}
  1393.   FTn = 1; // Set to truncate 1 super digits in display}
  1394.   FDn = (MuNMax - 2) * MuDMax; // Set max digits to display to max}
  1395.   ScieN = False;
  1396.   Expanded = False;
  1397.   FMCB = FMC;
  1398. } // --end-- InitMultiF
  1399.  
  1400. /*
  1401. Revisions made -
  1402. --------
  1403. Converted from Turbo Pascal 6.0 to Borland C++ 3.1 for Windows  HJS 1992-12-03
  1404. --------
  1405. */
  1406. // --end-- MultiFDW.Cpp  Multiple-precision floating decimal algorithms
  1407.