#undef  SPEEDUP				// 0 // 0 or 1 (default is 1)
#define SPEEDUP				// 1 // 0 or 1 (default is 1)

using System;
using System.IO;			// for StreamWriter I/O to a file
using System.Windows.Forms; // for RichTextBox

namespace MyNamespace
{
	/// <summary>
	///		Multiple precision integer algorithms
	///		Developed in TURBO Pascal 5.5 and 6.0, converted to Visual Basic 6
	///		and now in Visual C# .Net
	///
	/// ---------------- *** LEGAL NOTICE: *** ------------------
	///
	/// This program is free software; you can redistribute it and/or modify
	/// it under the terms of the GNU General Public License (GPL) as published
	/// by the Free Software Foundation. cf. the file COPYING.txt.
	///
	/// ------------ *** end of legal notice *** ------------------
	/// </summary>

	// ------- A multiple-precision unsigned integer
	public class MultiUI
	{
		public int N;		// Current length in super digits, 1..MuNMax
		public int M;		// Max number of super digits ever in this number
		public int[] V;		// string of super digits, Value, LSD at index o

		public MultiUI()	// initialize with M = 1
		{
			N = 1;			// constructor logic
			M = 1;
			V = new int[2];
		}

		public MultiUI(int nMaxI)	// initialize with M = nMaxI
		{
			N = 1;					// allocate m + 1 Super digits
			M = nMaxI;				// indexed from 0 thru m-1
			if (M < 1) M = 1;		// with index m for overflow calculations
			if (M > MultiID.mMax) M = MultiID.mMax;
			try
			{
				V = new int[M+1];
			}
			catch (Exception ex)
			{
				long totalMemory = GC.GetTotalMemory(false);
				string st = ex.Message + ", M = " + M + ", z set to zero, Memory used = " + totalMemory;
				Common.ReportError("MultiUI", st);
				M = 1;
				V = new int[M+1];
			}
		}

//		public static void operator = (MultiUI left, MultiUI right)		// not allowed to overload =
//		{
//			MultiID.SetToUI(right, left);
//		}

		public static MultiUI operator + (MultiUI left, MultiUI right)	// can be used for a += b
		{
			int m = Math.Max(left.N, right.N) + 1;
			MultiUI z = new MultiUI(m);
			MultiID.AddUI(left, right, z);
			return z;
		}

		public static bool operator == (int d1, MultiUI z)
		{
			return MultiID.EQ1UI(d1, z);
		}

		public static bool operator == (MultiUI z, int d1)
		{
			return MultiID.EQ1UI(d1, z);
		}

		public static bool operator != (int d1, MultiUI z)
		{
			return !MultiID.EQ1UI(d1, z);
		}

		public static bool operator != (MultiUI z, int d1)
		{
			return !MultiID.EQ1UI(d1, z);
		}

		public override bool Equals(object o)	// dummy
		{
			return true;
		}

		public override int GetHashCode()		// dummy
		{
			return 0;
		}
	}

	// ------- A multiple-precision signed integer
	public class MultiSI : MultiUI
	{
		public bool S;		// sign, true = negative

		public MultiSI()	// initialize with M = 0
			: base()		// calls MultiUI()
		{
			S = false;	// constructor logic
		}

		public MultiSI(int nMaxI)	// initialize with M = nMaxI
			: base(nMaxI)			// calls MultiUI(nMaxI)
		{
			S = false;
		}

//		public static void operator = (MultiSI left, MultiSI right)		// not allowed to overload =
//		{
//			MultiID.SetToSI(right, left);
//		}

		public static MultiSI operator & (MultiSI left, MultiSI right)	// use &= for = assignment
		{
			MultiID.SetToSI(right, left);
			return left;
		}

		public static MultiSI operator & (MultiSI left, int right)	// use &= for = assignment
		{
			MultiID.SetTo1SI(right, left);
			return left;
		}

		public static MultiSI operator & (MultiSI left, long right)	// use &= for = assignment
		{
			MultiID.SetToExSI(right, left);
			return left;
		}

		public static MultiSI operator & (MultiSI left, double right)	// use &= for = assignment
		{
			MultiID.SetToDbSI(right, left);
			return left;
		}

		public static MultiSI operator & (MultiSI left, decimal right)	// use &= for = assignment
		{
			MultiID.SetToDeSI(right, left);
			return left;
		}

		public static MultiSI operator + (MultiSI left, MultiSI right) // can be used for a += b
		{
			int m = Math.Max(left.N, right.N) + 1;
			MultiSI z = new MultiSI(m);
			MultiID.AddSI(left, right, z);
			return z;
		}

		public static MultiSI operator + (MultiSI left, int right) // can be used for a += b
		{
			MultiSI z = new MultiSI(left.N + 1);
			MultiID.Add1SI(left, right, z);
			return z;
		}

		public static MultiSI operator + (int left, MultiSI right) // can be used for a += b
		{
			MultiSI z = new MultiSI(right.N + 1);
			MultiID.Add1SI(right, left, z);
			return z;
		}

		public static MultiSI operator - (MultiSI left, MultiSI right) // can be used for a -= b
		{
			int m = Math.Max(left.N, right.N) + 1;
			MultiSI z = new MultiSI(m);
			MultiID.SubSI(left, right, z);
			return z;
		}

		public static MultiSI operator - (MultiSI left, int right) // can be used for a -= b
		{
			MultiSI z = new MultiSI(left.N + 1);
			MultiID.Add1SI(left, -right, z);
			return z;
		}

		public static MultiSI operator - (int left, MultiSI right) // can be used for a -= b
		{
			MultiSI z = new MultiSI(right.N + 1);
			MultiID.Add1SI(right, -left, z);
			MultiID.ChSSI(z);
			return z;
		}

		public static MultiSI operator - (MultiSI x)
		{
			MultiSI z = new MultiSI(x.N);
			MultiID.SetToSI(x, z);
			MultiID.ChSSI(z);
			return z;
		}

		public static MultiSI operator * (MultiSI left, MultiSI right) // can be used for a *= b
		{
			MultiSI z = new MultiSI(left.N + right.N);
			MultiID.MulSI(left, right, z);
			return z;
		}

		public static MultiSI operator * (MultiSI left, int right) // can be used for a *= b
		{
			MultiSI z = new MultiSI(left.N + MultiID.sdInInt);
			MultiID.Mul1SI(left, right, z);
			return z;
		}

		public static MultiSI operator * (int left, MultiSI right)
		{
			MultiSI z = new MultiSI(right.N + MultiID.sdInInt);
			MultiID.Mul1SI(right, left, z);
			return z;
		}

		public static MultiSI operator / (MultiSI left, MultiSI right) // can be used for a /= b
		{
			MultiSI q = new MultiSI(left.N);
			MultiID.remSI =	new MultiSI(right.N);
			MultiID.DivQRSI(left, right, MultiID.remSI, q);
			return q;
		}

		public static MultiSI operator / (MultiSI left, int right) // can be used for a /= b
		{
			MultiSI q = new MultiSI(left.N);
			MultiID.Div1QRSI(left, right, out MultiID.rem1, q);
			return q;
		}

		public static MultiSI operator % (MultiSI left, MultiSI right) // can be used for a %= b
		{
			MultiSI r =	new MultiSI(right.N);
			MultiID.ModSI(left, right, r);
			return r;
		}

		public static int operator % (MultiSI left, int right) // can be used for a %= b
		{
			MultiSI q = new MultiSI(left.N);
			int r;
			MultiID.Div1QRSI(left, right, out r, q);
			return r;
		}

		public static MultiSI operator ^ (MultiSI left, MultiSI right) // can be used for a ^= b
		{
			int m = MultiID.mMax;
			if (MultiID.useFMB && (MultiID.fMB != 0)) m = 2 * Math.Max(left.N, MultiID.fMB.N);
			MultiSI z = new MultiSI(m);
			if (MultiID.useFMB) MultiID.PowMBSI(left, right, MultiID.fMB, z);
			else				MultiID.PowMBSI(left, right, MultiID.sZero, z);
			return z;
		}

		public static MultiSI operator ^ (MultiSI left, int right) // can be used for a ^= b
		{
			int m = MultiID.mMax;
			if (right <= 2) m = 2 * left.N;
			if (MultiID.useFMB && (MultiID.fMB != 0)) m = 2 * Math.Max(left.N, MultiID.fMB.N);
			MultiSI z = new MultiSI(m);
			if (right == 2)
			{
				if (MultiID.useFMB) MultiID.PowMBSI(left, MultiID.sTwo, MultiID.fMB, z);
				else				MultiID.MulSI(left, left, z);
				return z;
			}
			MultiID.SetTo1SI(right, z);
			if (MultiID.useFMB) MultiID.PowMBSI(left, z, MultiID.fMB, z);
			else				MultiID.PowMBSI(left, z, MultiID.sZero, z);
			return z;
		}

		public static bool operator == (int d1, MultiSI z)
		{
			return MultiID.EQ1SI(d1, z);
		}

		public static bool operator != (int d1, MultiSI z)
		{
			return !MultiID.EQ1SI(d1, z);
		}

		public static bool operator == (MultiSI z, int d1)
		{
			return MultiID.EQ1SI(d1, z);
		}

		public static bool operator != (MultiSI z, int d1)
		{
			return !MultiID.EQ1SI(d1, z);
		}

		public static bool operator == (long ex, MultiSI z)
		{
			return MultiID.EQExSI(ex, z);
		}

		public static bool operator != (long ex, MultiSI z)
		{
			return !MultiID.EQExSI(ex, z);
		}

		public static bool operator == (MultiSI z, long ex)
		{
			return MultiID.EQExSI(ex, z);
		}

		public static bool operator != (MultiSI z, long ex)
		{
			return !MultiID.EQExSI(ex, z);
		}

		public static bool operator == (MultiSI x, MultiSI y)
		{
			return (MultiID.CompSI(y, x) == 0);
		}

		public static bool operator != (MultiSI x, MultiSI y)
		{
			return (MultiID.CompSI(y, x) != 0);
		}

		public static bool operator < (MultiSI x, MultiSI y)
		{
			return (MultiID.CompSI(y, x) < 0);
		}

		public static bool operator > (MultiSI x, MultiSI y)
		{
			return (MultiID.CompSI(y, x) > 0);
		}

		public static bool operator <= (MultiSI x, MultiSI y)
		{
			return (MultiID.CompSI(y, x) <= 0);
		}

		public static bool operator >= (MultiSI x, MultiSI y)
		{
			return (MultiID.CompSI(y, x) >= 0);
		}

		public static MultiSI operator ++ (MultiSI x)
		{
			MultiID.RAdd1SI(1, x);
			return x;
		}

		public static MultiSI operator -- (MultiSI x)
		{
			MultiID.RAdd1SI(-1, x);
			return x;
		}

		public override bool Equals(object o) // dummy
		{
			return true;
		}

		public override int GetHashCode() // dummy
		{
			return 0;
		}
	}

	public class MultiID
	{
		public MultiID() // constructor logic
		// initialize Multi-precision integer package
		// input: muDMax, muNMax, baseM
		{
			baseM = (int)Math.Pow(10, muDMax);
			baseMr = (int)Math.Pow(10, muDMax/2);
			baseML = baseM;
			maxDigit = baseM - 1;
			baseSq = Math.BigMul(baseM, baseM);
			bSMB = baseSq - baseM;
			muDpG = 5;
			muLenD = 10;
			muErr = false;
			muUnderflow = false;
			muErrRep = true;
			muAbort = false;
			softAbort = false;
			statReq = false;
			statReq2 = false;
			statReq3 = false;
			mMax = Common.nMax;
			maxDis = 200;
			maxDis = Math.Min(maxDis, mMax * muDMax);
			nDone = 0;
			nToGo = 0;
			time0 = timeP = Environment.TickCount/1000.0;
			tV1 = tV2 = 0;
			baseI = baseO = 10;
			nL = Environment.NewLine;
			nL2 = nL + nL;
			nLc = nL;					// line continuation mark for WriteT nL or "-" + nL
			muICount = 0;				// belongs to MuInterrupt()
			muITimeSaved = 0;			// belongs to MuInterrupt()
			fMB = new MultiSI(1);		// Multi-precision modulo base, 0 => normal arithmetic.
			remSI = new MultiSI(1);		// holds the remainder of infix operator SI / SI
			sMOne = new MultiSI(1);
			SetTo1SI(-1, sMOne);		// sMOne   = -1, minus works in SetTo1SI
			sZero = new MultiSI(1);		// sZero   =  0, zero is default value
			sOne = new MultiSI(1);
			sOne.V[0] = 1;				// sOne    =  1, direct method
			sTwo = new MultiSI(1);
			SetTo1SI(2, sTwo);			// sTwo    =  2, function method
			sThree = new MultiSI(1);
			SetTo1SI(3, sThree);		// sThree  =  3
			sFour = new MultiSI(1);
			SetTo1SI(4, sFour);			// sFour   =  4
			sFive = new MultiSI(1);
			SetTo1SI(5, sFive);			// sFive   =  5
			s100 = new MultiSI(1);
			SetTo1SI(100, s100);		// s100    =  100
			dblMaxValue = Math.Pow(2, 53) - 1; // 2^53-1 = 9,00719,92547,40991
			decMaxValue = Decimal.MaxValue; // 2^96-1 = 7922,81625,14264,33759,35439,50335
			sdInInt = (int)Math.Ceiling(Math.Log(Int32.MaxValue, baseM));
			sdInLong = (int)Math.Ceiling(Math.Log(Int64.MaxValue, baseM));
			sdInDbl = (int)Math.Ceiling(Math.Log(dblMaxValue, baseM));
			sdInDec =(int)Math.Ceiling(Math.Log((double)decMaxValue, baseM));
			d1LineOn = false;
			dBoolOn = true;
			dDigOn = true;
			diagOn = false;
			debug = false;
			useFMB = false;
			max64 = new MultiSI(sdInLong);
			SetToExSI(Int64.MaxValue, max64);
			min64 = new MultiSI(sdInLong);
			SetToExSI(-Int64.MaxValue, min64);
			writeBool = true;
		}

		private static void IpAdd(MultiUI x, MultiUI y, out int k) // |y| = |y| + |x|
		// Multi-precision integer add absolute values in place y is changed
		// y.N Must be >= x.N, x.S and y.S not used or changed, x not changed
		// overflow super digit k returned not zero if and only if y overflows
		// muErr flag not set. Address of x and y can be the same, then x is changed
		{
			if (muErr || muAbort)
			{
				ClearUI(y);
				k = 0;
				return;
			}
			int t;
			int i = 0;
			k = 0;
			while (i < x.N)
			{
				t = x.V[i] + y.V[i] + k;	// add next super digit with carry
				if (t >= baseM)
				{
					k = 1;
					y.V[i] = t - baseM;
				}
				else
				{
					k = 0;
					y.V[i] = t;
				}
				i++;
			}
			if (k == 0) return;
			while (i < y.N) // x = 0 from here on
			{
				if (y.V[i] == maxDigit)		// k = 1, carry to next super digit of y
					y.V[i] = 0;
				else
				{
					y.V[i]++;
					k = 0;
					return;					// return k = 0
					}
				i++;
			}
			y.N++;
			if (y.N > y.M)					// if overflow return k = 1
			{
				y.N = y.M;
				NormUI(y);
			}
			else
			{
				y.V[y.N-1] = 1;
				k = 0; // return k = 0
			}
		} // IpAdd

		private static void IpSub(MultiUI x, MultiUI y) // |y| = |y| - |x|
		// Multi-precision integer subtract absolute values in place, y is changed
		// |y| Must be >= |x|, x.S and y.S not used or changed, x not changed
		// muErr flag not set. Address of x and y can be the same, then x is changed
		{
			if (muErr || muAbort)
			{
				ClearUI(y);
				return;
			}
			int t;
			int i = 0;
			int k = 0;
			while (i < x.N)
			{
				t = y.V[i] - x.V[i] + k;	// subtract next super digit with borrow
				if (t < 0)
				{
					k = -1;
					y.V[i] = t + baseM;
				}
				else
				{
					k = 0;
					y.V[i] = t;
				}
				i++;
			}
			if (k == 0 || x.N == y.N)
			{
				NormUI(y);
				return;
			}
			while (i < y.N)					// y = 0 from here on
			{
				if (y.V[i] == 0)			// k = -1, borrow from next super digit of y
					y.V[i] = maxDigit;
				else
				{
					y.V[i]--;
					NormUI(y);
					return;
				}
				i++;
			}
		} // IpSub

		private static void IpMul1(int d1, MultiUI x, out int k) // x = |d1| * x
		// Multi-precision one super digit multiply in place, x is changed
		// k = overflow super digit returned > zero if and only if x overflows
		// signs not used
		{
			if (muErr || muAbort)
			{
				ClearUI(x);
				k = 0;
				return;
			}
			int i;
			long t, rem;
			int aD1 = Math.Abs(d1);
			if (aD1 > maxDigit)
			{
				MuWriteErr("IpMul1: Multiplier too large");
				ClearUI(x);
				k = 0;
				return;
			}
			long kL = 0;
			for (i = 0; i < x.N; i++)
			{
				t = Math.BigMul(aD1, x.V[i]) + kL;		// multiply next super digit and add carry
				kL = Math.DivRem(t, baseML, out rem);	// k = t/base
				x.V[i] = (int)rem;
			}
			k = (int)kL;
			if (k == 0) return;
			x.N++;
			if (x.N > x.M)					// if overflow return k > 0
				x.N = x.M;
			else
			{
				x.V[x.N-1] = k;
				k = 0;
			}
		} // IpMul1

		private static void IpMul(MultiUI x, MultiUI y, MultiUI z) // z = x * y
		// x or y may NOT have the same address in memory as z
		// x and y may have the same address, they are not changed
		// for efficiency have x.N >= y.N
		// sets muErr true if overflow
		{
			if (muErr || muAbort)
			{
				ClearUI(z);
				return;
			}
			int i;
			int j;
			long t, rem;
			int iJ = 0;
			long k = 0;
			for (i = 0; i < x.N; i++)
				z.V[i] = 0;					// clear lower of z
			statReq = false;
			for (j = 0; j < y.N; j++)
			{
				k = 0;
				i = 0;
				while (i < x.N)
				{
					iJ = i + j;
					if (iJ < z.M)			// multiply, add and carry next Digit
					{
						t = Math.BigMul(x.V[i], y.V[j]) + z.V[iJ] + k;
						k = Math.DivRem(t, baseML, out rem);	// k = t/base
						z.V[iJ] = (int)rem;
					}
					else
					{
						i = x.N;			// break Out of inner loop
					}
					i++;
					MuInterrupt();
					if (muAbort)
					{
						ClearUI(z);
						return;
					}
					if (statReq)
						WriteStatus("Integer Multiply: " + (Math.Floor(100.0 * (j + 1) / y.N)).ToString() + " Percent done");
				}
				iJ++;
				if (iJ < z.M) z.V[iJ] = (int)k;
			}
			if (k != 0) iJ++;
			if (iJ > z.M)
			{
				MuWriteErr("Multiplication overflow, continuing... (N=" +
							iJ + " > M=" + z.M + ")");
				iJ = z.M;
			}
			z.N = iJ;
			NormUI(z);
		} // IpMul

		private static void IpDiv1(int d1, out int r1, MultiSI u) // u = u / d1
		// Multi-precision one super digit divide in place, u is changed // r1 = Remainder
		// allows literal d1, muErr set true if d1 = 0, sign(r1) = sign(d1)
		{
			if (muErr || muAbort)
			{
				ClearSI(u);
				r1 = 0;
				return;
			}
			long t, rem;
			int i;
			long aD1;
			r1 = 0;
			if (d1 == 0)
			{
				MuWriteErr("Cannot divide by zero, continuing...");
				return;
			}
			MuInterrupt();
			if (muAbort)
			{
				ClearSI(u);
				return;
			}
			aD1 = Math.Abs(d1);
			i = u.N - 1;
			while (i >= 0)
			{
				t = Math.BigMul(r1, baseM) + u.V[i];		// divide next super digit and save remainder
				u.V[i] = (int)Math.DivRem(t, aD1, out rem);	// V[i] = t/aD1
				r1 = (int)rem;
				i--;
			}
			if (u.S) r1 = -r1;
			u.S = u.S ^ (d1 < 0);			// ^ is Xor, exclusive Or
			if (u.S && (r1 != 0))
			{
				r1 += d1;
				RAdd1SI(-1, u);
			}
			NormSI(u);
		} // IpDiv1

		private static void IpDivSlow(MultiSI u, MultiSI d, MultiSI r, MultiSI q)
		// q = u / d, r = Remainder, Slow Method (but fast for small numbers)
		// u and d are not changed, u is moved to r and then r is divided
		// in place to give q and a remainder in r, it is ok for u and r
		// to have the same location in memory, but then u will be changed
		// sets muErr true if d = 0, sign(r) = sign(d)
		// See Knuth Vol. 2 of "The Art..." Algorithm d, pg. 237.
		{
			if (muErr || muAbort)
			{
				ClearSI(q);
				ClearSI(r);
				return;
			}
			long rem;
			int v1;
			int v2;
			int d1;
			int k;
			long kL;
			int qH;
			long t;
			int i;
			int j;
			int iJ;
			if (d.N == 1)			// use IpDiv1 for a one super digit divide
			{
				SetToSI(u, q);
				d1 = d.V[0];
				if (d.S && (d1 != 0)) d1 = -d1;
				IpDiv1(d1, out k, q);
				r.S = (k < 0);
				r.V[0] = Math.Abs(k);
				r.N = 1;
			}
			else
			{
				SetToSI(u, r);		// SetToSI is noOp if same address
				if (r.N < d.N)
				{
					ClearSI(q);
					return;
				}
				q.N = r.N - d.N + 1;
				j = r.N;
				r.V[r.N] = 0;
				d.V[d.N] = 0;
				v1 = d.V[d.N-1];
				v2 = d.V[d.N-2];
				d1 = baseM / (v1 + 1);	// d1. normalize
				if (d1 > 1)
				{
					IpMul1(d1, d, out k);
					v1 = d.V[d.N-1];
					v2 = d.V[d.N-2];
					IpMul1(d1, r, out k);
					if (k != 0) r.V[j] = k;
				}											// d2. init j done above
				statReq = false;
				while (j >= d.N)							// d3. calculate q hat
				{
					t = Math.BigMul(r.V[j], baseM) + r.V[j-1];
					if (v1 == r.V[j])
						qH = maxDigit;
					else
						qH = (int)(t / v1);
					while ( Math.BigMul(v2, qH) > (Math.BigMul((int)(t - Math.BigMul(qH, v1)), baseM) + r.V[j-2]) )
					{	qH--; }
					// qH =qH+1 // diag for add back test
					kL = baseM;
					for (i = 0; i <= d.N; i++)				// d4. multiply and subtract
					{
						iJ = i + j - d.N;
						t = r.V[iJ] - Math.BigMul(qH, d.V[i]) + kL + bSMB; // + baseSq - baseM
						kL = Math.DivRem(t, baseM, out rem);
						r.V[iJ] = (int)rem;
					}
					q.V[j - d.N] = qH;						// d5. test remainder
					if (kL == maxDigit)
					{
						q.V[j - d.N] = qH - 1;				// d6. add back
						k = 0;
						iJ = j - d.N;
						for (i = 0; i <= d.N; i++)
						{
							t = d.V[i] + r.V[iJ] + k;
							k = (int)Math.DivRem(t, baseM, out rem);
							r.V[iJ] = (int)rem;
							iJ++;
						}
					}
					j--;									// d7. loop on j
					MuInterrupt();
					if (muAbort)
					{
						ClearSI(q);
						ClearSI(r);
						return;
					}
					if (statReq)
						WriteStatus("Integer Divide: " + (Math.Floor(100.0 * (r.N - j) / (r.N - d.N + 1))).ToString() + " Percent done");
				}
				q.S = u.S ^ d.S; // exclusive or
				r.S = u.S;
				NormSI(q);
				NormSI(r);
				if (d1 > 1)									// d8. un-normalize
				{
					IpDiv1(d1, out k, r);
					IpDiv1(d1, out k, d);
				}
				if (r.S != d.S && r != 0)
				{
					RAddSI(d, r);	
					RAdd1SI(-1, q);	
				}
			
			}
			NormSI(q);
			NormSI(r);
		} // IpDivSlow

		public static void ReducePrecision(int ps, MultiSI d, MultiSI dL)
		// remove lower ps super digits, d and dL may be the same object
		{

			int j = d.N - ps;
			if (j <= 0)
				ClearSI(dL);
			else
			{
				if (dL.M < j) ReSizeSI(j, dL);
				for (int i = ps; i < d.N; i++)
					dL.V[i-ps] = d.V[i];			// dL is d at a reduced precision
				dL.N = j;
			}
		}

#if SPEEDUP
		private static void IpDivFast(MultiSI u, MultiSI d, MultiSI r, MultiSI q)
		// q = u / d, r = Rem, Fast method for large numbers, sign(r) = sign(d)
		// u and d are not changed, u is divided by d
		// to give q and a remainder r, it is OK for u and r
		// to have the same location in memory, but then u will be changed
		// Sets MuErr true if d = 0
		{
			if (muErr || muAbort)
			{
				ClearSI(q);
				ClearSI(r);
				return;
			}
			int i, j, pl, po, ps;
			if ( ((d.N <= 2 || u.N - d.N <= 2) && !forceFHT) || !allowFHT )
			{
				IpDivSlow(u, d, r, q);
				return;
			}
			MultiSI mm = new MultiSI();			// m * m
			MultiSI mmd = new MultiSI();		// (m*m) * d
			pl = 5;
			if (pl > (u.N - d.N))  pl = u.N - d.N + 1;
			j = (pl+1 > d.N) ? d.N : pl+1;
			MultiSI dL = new MultiSI(j);
			ps = d.N-j;
			ReducePrecision(ps, d, dL);			// dL is d at a reduced precision
			MultiSI e = new MultiSI(pl+1+dL.N);	// e is 1.0 for 1 / d and Error in q
			e.N = e.M;
			e.V[ e.N-1] = 1;
			MultiSI m = new MultiSI(pl+12);		// m is 1 / d and 2*m
			IpDivSlow(e, dL, e, m);
			if (e.M < 5)
				ReSizeSI(5, e);
			DoneSI(dL);
			while (pl <= (u.N - d.N))			// m = 2*m - d*m^2 = m + m - (m*m) * d
			{
				po = pl;  pl += pl;				// Should use m = m + m * (1 - m*d)
				if (pl > (u.N - d.N))  pl = (u.N - d.N) + 1;
				int p;
				for(p = u.N - d.N + 1;  p > pl;  p=(p+1)/2);
				pl = p;
				ReSizeSI(m.N + m.N, mm);
				MulSI(m, m, mm);				// mm = m * m
				ps = po+po-pl;
				ReducePrecision(ps, mm, mm);
				j = (pl+1 > d.N) ? d.N : pl+1;
				ReSizeSI(j, dL);
				ps = d.N-j;
				ReducePrecision(ps, d, dL);		// dL is d at a reduced precision
				ReSizeSI(mm.N + dL.N, mmd);
				MulSI(mm, dL, mmd);				// mmd = (m*m) * d
				ps = dL.N;
				ReducePrecision(ps, mmd, mmd);
				DoneSI(mm);
				DoneSI(dL);
				RAddSI(m, m);					// m = m + m
				j = m.N - 1;
				i = j + pl - po;
				Common.ReDimPreserve(ref m.V, i+2); // preserve and re-dimention arrays
				m.N = m.M = i+1;
				for ( ;  i >= 0;  i--, j--)
				{
					if (j >= 0)
						m.V[i] = m.V[j];
					else
						m.V[i] = 0;
				}
				RSubSI(mmd, m);					// m = (m + m) - ((m*m) * d)
				DoneSI(mmd);
			}
			ps = d.N - 1;
			MultiSI uL = new MultiSI(u.N - ps);		// uL = Upper part of u (Points to u.V)
			ReducePrecision(ps, u, uL);
			MultiSI qL = new MultiSI(uL.N + m.N);	// qL = u * m
			MulSI(uL, m, qL);
			DoneSI(uL);
			DoneSI(m);
			ps = pl+1;
			ReducePrecision(ps, qL, qL);
			int ii = Math.Max(d.N + qL.N, u.N);
			MultiSI rL = new MultiSI(ii);			// rL = Local remainder
			MulSI(d, qL, rL);
			SubSI(u, rL, rL);
			if (muAbort || muErr)
			{
				ClearSI(q);
				ClearSI(r);
				return;
			}
			if (!ZTestSI(rL))					// if Rem != 0
			{
				IpDivSlow(rL, d, rL, e);		// now 0 < rL < |d|
				RAddSI(e, qL);
			}
			SetToSI(qL, q);
			SetToSI(rL, r);
			DoneSI(rL);
			DoneSI(qL);
			DoneSI(e);
			NormSI(q);
			NormSI(r);
		} // IpDivFast
#endif // SPEEDUP

		private static void IpDiv(MultiSI u, MultiSI d, MultiSI r, MultiSI q)
		// q = u / d, r = Rem, sign(r) = sign(d)
		// u and d are not changed, u is divided by d
		// to give q and a remainder r, it is OK for u and r
		// to have the same location in memory, but then u will be changed
		// Sets MuErr true if d = 0
		{
			if (muErr || muAbort)
			{
				ClearSI(q);
				ClearSI(r);
				return;
			}
#if SPEEDUP
			float prod = (float)(u.N - d.N) * d.N;
			if ( (prod / u.N < 344.0 && !forceFHT) || !allowFHT )
				IpDivSlow(u, d, r, q);
			else
				IpDivFast(u, d, r, q);
#else
			IpDivSlow(u, d, r, q);
#endif // SPEEDUP
		}

		public static void DoneUI(MultiUI z)	// free up space, set to zero
		{
			z.V = new int[2];					// 1 + spare Super digits
			z.M = 1;
			z.N = 1;
		}

		public static void ClearUI(MultiUI z)	// set Value = 0
		{
			z.V[0] = 0;
			z.N = 1;
		}

		public static void ReSizeUI(int m, MultiUI z)
		// re-dimension UI, m is new length, value = 0
		{
			if (m == z.M && z.V != null && m+1 == z.V.GetLength(0))
				return;
			z.M = m;
			z.N = 1;
			try
			{
				z.V = new int[m+1];		// re-dimention V
			}
			catch (Exception ex)
			{
				long totalMemory = GC.GetTotalMemory(false);
				string st = ex.Message + ", M = " + m + ", M set to one, Memory used = " + totalMemory;
				Common.ReportError("ReSizeUI", st);
				z.V = new int[2];		// z.M + 1 Super digits
			}
		}

		public static void SetToUI(MultiUI x, MultiUI z) // z = x (final z.M >= x.M)
		{
			int i;
			// how do you tell if two variables are the same? answer: Object.Equals(x, z)
			// or Object.ReferenceEquals(x, z) or x == z, but ReferenceEquals preferred here
			if (Object.ReferenceEquals(x, z)) return;
			z.N = x.N;
			if (z.N > z.M)
			{
				z.V = new int[z.N+1];	// resize V (ReDim)
				z.M = z.N;
			}
			for (i = 0; i < z.N; i++)
			{
				z.V[i] = x.V[i];
			}
		}

		public static void NormUI(MultiUI z) // normalize
		{
			int i;
			if (z.V.Length < 2) DoneUI(z);
			if (z.N < 1) ClearUI(z);
			if (z.N >= z.V.Length) z.N = z.V.Length - 1;
			i = z.N - 1;						// delete leading zero super digits
			while (z.V[i] == 0 && i > 0)
				i--;
			z.N = i + 1;
		}

		public static int AbCompUI(MultiUI x, MultiUI z)
		// AbCompUI = sign(|z| - |x|)
		// sign = -1 if |z| < |x|
		// sign = 0 if |z| = |x|
		// sign = +1 if |z| > |x|
		// x.S and y.S not used or changed, z and x not changed
		{
			int i;
			if (z.N < x.N) return -1;
			if (z.N > x.N) return 1;
			i = z.N - 1;
			while (z.V[i] == x.V[i] && i > 0)
				i--;
			if (z.V[i] < x.V[i]) return -1;
			if (z.V[i] > x.V[i]) return 1;
			return 0;
		}

		public static void RAddUI(MultiUI x, MultiUI z) // z = z + x
		// sets muErr true if overflow
		{
			int k;
			int i;
			if (z.N < x.N)					// if digits in z < digits in x
			{
				if (z.M >= x.N)
				{
					for (i = z.N; i < x.N; i++)
						z.V[i] = 0;
					z.N = x.N;
				}
				else
				{
					MuWriteErr("X too big to add, continuing... (N=" +
								x.N + " > M=" + z.M + ")");
					return;
				}
			}
			IpAdd(x, z, out k);
			if (k > 0) MuWriteErr("Addition overflow, continuing... (M=" + z.M + ")");
		}

		public static void RSubUI(MultiUI x, MultiUI z) // z = z - x
		// sets muErr true if n > z
		{
			int i = AbCompUI(x, z);
			if (i == 0)
				ClearUI(z);
			else if (i < 0)						// if z < x, Error
				MuWriteErr("Unsigned subtraction error, continuing...");
			else
				IpSub(x, z);
		}

		public static void AddUI(MultiUI x, MultiUI y, MultiUI z) // z = x + y
		// sets muErr true if overflow
		{
			if (Object.ReferenceEquals(z, x))
				RAddUI(y, z);
			else if (Object.ReferenceEquals(z, y))
				RAddUI(x, z);
			else if (x.N >= y.N)
			{
				SetToUI(x, z);
				RAddUI(y, z);
			}
			else
			{
				SetToUI(y, z);
				RAddUI(x, z);
			}
		}

		public static void SubUI(MultiUI x, MultiUI y, MultiUI z) // z = x - y
		// sets muErr true if overflow
		{
			if (Object.ReferenceEquals(x, y))
				ClearUI(z);
			else if (Object.ReferenceEquals(z, x))
				RSubUI(y, z);
			else if (Object.ReferenceEquals(z, y))
			{
				MultiUI yL = new MultiUI(y.N);
				SetToUI(y, yL);
				SetToUI(x, z);
				RSubUI(yL, z);
			}
			else
			{
				SetToUI(x, z);
				RSubUI(y, z);
			}
		}

		public static bool ZTestUI(MultiUI z) // true if z is Zero
		{
			return (z.N == 1) && (z.V[0] == 0);
		}

		public static bool ZTestSI(MultiSI z) // true if z is Zero
		{
			if ((z.N != 1) || (z.V[0] != 0)) return false;
			z.S = false;	// side effect, prevents -0
			return true;
		}

		public static bool EQ1UI(int d1, MultiUI z) // true if |z| = d1
		// z, d1 are not changed, allows literal d1
		{
			return (z.N == 1) && (z.V[0] == d1);
		}

		public static bool EQ1SI(int d1, MultiSI z) // true if z = d1
		// z, d1 are not changed, allows literal d1
		{
			return (z.N == 1) && (z.V[0] == Math.Abs(d1)) && (z.S == (d1 < 0));
		}

		public static bool EQExSI(long ex, MultiSI z) // true if z = ex
		// z, ex are not changed, allows literal ex
		{
			int exN = sdInLong;
			if (z.N > exN) return false;
			MultiSI x = new MultiSI(exN);
			SetToExSI(ex, x);
			bool b = (CompSI(x, z) == 0);
			DoneSI(x);
			return b;
		}

		public static bool GT1UI(int d1, MultiUI z) // true if |z| > d1
		// z, d1 are not changed, allows literal d1
		{
			return (z.N > 1) || (z.V[0] > d1);
		}

		public static bool GTDbSI(double db, MultiSI z) // true if z > db
		// z, db are not changed, allows literal db
		{
			MultiSI DbSI = new MultiSI(sdInDbl); // up to 16 digits
			SetToDbSI(db, DbSI);
			bool b = (z > DbSI);
			DoneSI(DbSI);
			return b;
		}

		public static void ReSizeSI(int m, MultiSI z)
		// re-dimension SI, m is new length, value = 0
		{
			ReSizeUI(m, z);
			z.S = false;
		}

		public static void DoneSI(MultiSI z) // free up space
		{
			DoneUI(z);
			z.S = false;
		}

		public static void ClearSI(MultiSI z) // set Value = 0
		{
			ClearUI(z);
			z.S = false;
		}

		public static void ChSSI(MultiSI z) // change sign, z = - z
		{
			if (!ZTestSI(z)) z.S = !z.S;
		}

		public static void SetToSI(MultiSI x, MultiSI z) // z = x
		{
			SetToUI(x, z);
			z.S = x.S;
		}

		public static void SwapSI(ref MultiSI x, ref MultiSI y) // temp = x, x = y, y = temp
		{
			MultiSI temp; // only pointers are changed
			temp = x;
			x = y;
			y = temp;
		}

		public static void SetToSIS(MultiSI x, MultiSI z) // z = x
		// copy everything but .Nm name, but don't make it smaller
		{
			try {int zM = z.M;}
			catch (Exception ex)
			{
				long totalMemory = GC.GetTotalMemory(false);
				string st = ex.Message + ", z set to zero, Memory used = " + totalMemory;
				Common.ReportError("SetToSIS", st);
				DoneSI(z);
				return;
			}
			int m = x.M;				// z = x; but you can't do it that way
			if (m > z.M)
			{
				try
				{
					z.V = new int[m+1];	// resize V (ReDim)
					z.M = m;
				}
				catch (Exception ex)
				{
					long totalMemory = GC.GetTotalMemory(false);
					string st = ex.Message + ", M = " + m + ", z set to zero, Memory used = " + totalMemory;
					Common.ReportError("SetToSIS", st);
					DoneSI(z);
					return;
				}
			}
			SetToSI(x, z);
		}

		public static void SetToSIN(MultiSI x, MultiSI z) // z = x, final z.M = x.N
		// copy everything but .Nm name
		{
			int m = x.N;
			if (m != z.M)
			{
				z.V = new int[m+1];	// resize V (ReDim)
				z.M = m;
			}
			SetToSI(x, z);
		}

		public static void NormSI(MultiSI z) // normalize
		{
			NormUI(z);
			if (z.N == 1 && z.V[0] == 0)
				z.S = false;	// prevent -0
		}

		public static void RAddSI(MultiSI x, MultiSI z) // z = z + x
		// sets muErr true if overflow
		{
			int i;
			if (z.S == x.S)							// same signs
				RAddUI(x, z);
			else									// opposite signs
			{
				i = AbCompUI(x, z);
				if (i < 0)							// if |z| < |x|, |z| = |x| - |z|
				{
					MultiUI temp = new MultiUI(z.N);
					SetToUI(z, temp);				// temp = z
					SetToSI(x, z);					// z = Large x
					IpSub(temp, z);					// z = z - temp (z = old x - old z)
				}
				else if (i > 0)						// if |z| > |x|, |z| = |z| - |x|
				{
					IpSub(x, z);
				}
				else
				{
					ClearSI(z);						// if |z| = |x|, z = 0
				}
			}
		}

		public static void RSubSI(MultiSI x, MultiSI z) // z = z - x
		// sets muErr true if overflow
		{
			if (Object.ReferenceEquals(z, x))
				ClearSI(z);
			else
			{
				ChSSI(x);
				RAddSI(x, z);
				ChSSI(x);
			}
		}

		public static void AddSI(MultiSI x, MultiSI y, MultiSI z) // z = x + y
			// sets muErr true if overflow
		{
			int i;
			if (Object.ReferenceEquals(z, x))
				RAddSI(y, z);
			else if (Object.ReferenceEquals(z, y))
				RAddSI(x, z);
			else if (x.S == y.S)
			{
				AddUI(x, y, z);
				z.S = x.S;
			}
			else
			{
				i = AbCompUI(y, x);		// signs are different
				if (i == 0)
					ClearSI(z);
				else if (i > 0)
				{
					SetToSI(x, z);
					IpSub(y, z);		// |x| > |y|, so |x| - |y|
				}
				else
				{
					SetToSI(y, z);
					IpSub(x, z);
				}
			}
		}

		public static void SubSI(MultiSI x, MultiSI y, MultiSI z) // z = x - y
		// sets muErr true if overflow
		{
			if (Object.ReferenceEquals(x, y))
				ClearSI(z);
			else if (Object.ReferenceEquals(z, x))
				RSubSI(y, z);
			else if (Object.ReferenceEquals(z, y))
			{
				RSubSI(x, z);
				ChSSI(z);				// prevent -0
			}
			else
			{
				ChSSI(y);
				AddSI(x, y, z);
				ChSSI(y);
			}
		}

		public static void ValueSI10(string st, MultiSI z)
		// convert String to z, returns used = number of character used
		// allows st to be a literal, muErr set true if overflow
		// works for base 10 only and st must be all digits, no commas or spaces
		{
			ClearSI(z);
			int leng = st.Length;
			if (leng == 0) return;
			int sD;
			if (st[0] == '-')
			{
				st = st.Substring(1);
				leng--;
				z.S = true;
			}
			int r;
			int q = Math.DivRem(leng, muDMax, out r);
			int n = q;
			if (r > 0) n++;
			if (n > z.M)
			{
				MuWriteErr("Input number overflow, continuing... (M=" + z.M + ")");
				return;
			}
			z.N = n;
			if (r > 0)
			{
				sD = Convert.ToInt32(st.Substring(0, r));
				z.V[q] = sD;
			}
			for (int i = q - 1; i >= 0; i--)
			{
				sD = Convert.ToInt32(st.Substring(r, muDMax));
				z.V[i] = sD;
				r += muDMax;
			}
			NormSI(z);
		} // ValueSI10

		public static void ValueSI(string st, out int used, MultiSI z)
		// convert String to z, returns used = number of character used
		// allows st to be a literal, muErr set true if overflow
		// works for all bases [2, 36]
		{
			char ch;
			int k;
			int j;
			int C;
			int d;
			bool fini;
			bool overflow;
			string maxD = MaxDig(baseI);	// max Input digit "1" to "Z"
			int leng = st.Length;
			ClearSI(z);
			used = 0;
			if (leng == 0) return;
			fini = false;
			overflow = false;
			j = 0;
			while (st[j] == ' ')
				j++;						// HJS 91/04/24, 91/09/24
			if (st[j] == '-')
			{
				j++;
				z.S = true;
			}
			else if (st[j] == '+')
				j++;
			C = 'A' - '0' - 10;
			while (!fini && (j < leng))
			{
				ch = Char.ToUpper(st[j]);	// lower case allowed for screen input
				d = ch - '0';
				if (d > 9)
				{
					d -= C;
					if (d < 10) d = -1;
				}
				if (0 <= d && d < baseI)
				{
					if (z.N == z.M && z.V[z.N-1] >= baseM / baseI)
					{
						overflow = true;
						fini = true;
					}
					else
					{
						IpMul1(baseI, z, out k);
						if (k > 0) overflow = true;
						z.V[0] += d;
					}
				}
				else
				{
					if (st[j] != ',')			// only extra ,s allowed for screen input
					{
						j--;
						fini = true;
					}
				}
				j++;
			}
			while (j > 0 && st[j-1] == ',')
				j--;							// don't use trailing commas
			used = j;							// return length used
			if (overflow)
				MuWriteErr("Input number overflow, continuing... (M=" + z.M + ")");
			NormSI(z);
		} // ValueSI

		public static void SetTo1SI(int d1, MultiSI z) // z = d1
		// allows a literal d1, d1 an integer, d1 not changed
		{
			if (d1 < baseM && d1 > -baseM) // protect against d1 = Int32.MinValue
			{
				int aD1 = Math.Abs(d1);
				z.V[0] = aD1;
				z.N = 1;
				z.S = (d1 < 0) && (aD1 != 0);	// prevent -0
			}
			else
			{
				SetToExSI((Int64)d1, z);
			}
		}

		public static void SetToDeSI(decimal de, MultiSI z)	// z = ext
		// allows a literal ext, ext not changed, max decimal integer for a decimal is 2^96-1 = 7922,81625,14264,33759,35439,50335
		// was Extended then Variant, then object, then Int64 now decimal
		{
			string st = de.ToString();
			ValueSI10(st, z);
		}

		public static void SetToExSI(Int64 ext, MultiSI z)	// z = ext
		// allows a literal ext, ext not changed, max Int64 = 9223,37203,68547,75807 = 2^63-1
		// was Extended then Variant, then object now Int64
		{
			string st = ext.ToString();
			ValueSI10(st, z);
		}

		public static void SetToDbSI(double db, MultiSI z) // z = db
		// allows a literal db, db not changed
		// db > 2^53-1 = 9,00719,92547,40991 may not be totally accurate
		{
			if (Math.Abs(db) <= Int64.MaxValue) // max long
			{
				Int64 ext = (Int64)db;
				SetToExSI(ext, z);
				return;
			}
			z.S = (db < 0);
			db = Math.Abs(db);
			for (z.N = 0; db > 0; z.N++)
			{
				if (z.N >= z.M)
				{
					Common.ReportError("SetToDbSI", "double too large");
					ClearSI(z);
					return;
				}
				double q = Math.Floor(db / baseM);
				z.V[z.N] = (int)(db - q * baseM);
				db = q;
			}
		}

		public static void Get1SI(out int d1, MultiSI z) // d1 = z % baseM
		// z is not changed
		{
			d1 = z.V[0];
			if (z.S) d1 = -d1;
		}

		public static void GetIntSI(out int i, MultiSI z) // i = z
		// z is not changed
		{
			if (z.N == 1)
			{
				i = z.V[0];
				if (z.S) i = -i;
				return;
			}
			long ext = z.V[1] * baseM + z.V[0];
			if (z.S) ext = -ext;
			if (ext < Int32.MinValue || ext > Int32.MaxValue || z.N > 2)
			{
				Common.ReportError("GetIntSI", "Value out of bounds");
				i = 0;
				return;
			}
			i = (int)ext;
		}

		public static void GetDbSI(out double db, MultiSI z) // db = z
		// z is not changed
		// db > 2^53-1 = 9,00719,92547,40991 may not be totally accurate
		{
			try
			{
				int i = 1;
				double t = 1.0;
				db = z.V[0];
				while (i < z.N)
				{
					t *= baseM;
					db += t * z.V[i];		// add next super digit
					i++;
				}
				if (z.S) db = -db;
			}
			catch (Exception ex)
			{
				Common.ReportError("GetDbSI", ex.Message);
				db = 0;
			}
		}

		public static void GetExSI(out Int64 ext, MultiSI z) // ext = z % baseM^3
		// z is not changed, Added procedure GetEx HJS 91/05/02
		// max Int64 = 9223,37203,68547,75807 = 2^63-1
		{
			try
			{
				if (z < min64 || z > max64)
					throw(new ArgumentOutOfRangeException());
				int i = 1;
				Int64 t = 1;
				ext = (Int64)(z.V[0]);
				while (i < z.N)
				{
					t *= baseML;
					ext += t * (Int64)z.V[i];	// add next super digit
					i++;
				}
				if (z.S) ext = -ext;
			}
			catch (Exception ex)
			{
				Common.ReportError("GetExSI", ex.Message);
				ext = 0;
			}
		}

		public static void GetDeSI(out decimal de, MultiSI z) // ext = z
		// z is not changed, Added procedure GetEx HJS 91/05/02
		// max integer decimal is 2^96-1 = 7922,81625,14264,33759,35439,50335
		{
			try
			{
				int i = 1;
				decimal t = 1;
				de = (decimal)(z.V[0]);
				int ii = 1 + 28 / muDMax;
				while (i < z.N && i < ii)
				{
					t *= (decimal)baseM;
					de += t * (decimal)z.V[i];	// add next super digit
					i++;
				}
				if (z.S) de = -de;
			}
			catch (Exception ex)
			{
				Common.ReportError("GetExSI", ex.Message);
				de = 0;
			}
		}

		public static void RAdd1SI(int d1, MultiSI z) // z = z + d1
		// allows literal d1, d1 is not changed, sets muErr true if overflow
		{
			MultiSI dL = new MultiSI(sdInInt);
			SetTo1SI(d1, dL);
			RAddSI(dL, z);
			DoneSI(dL);
		}

		public static void Add1SI(MultiSI x, int d1, MultiSI z) // z = x + d1
		// allows literal d1, d1 is not changed, sets muErr true if overflow
		{
			SetToSI(x, z);
			RAdd1SI(d1, z);
		}

		public static int SignSI(MultiSI z) // sign(z) = -1 if z < 0, 0 if z = 0, +1 if z > 0
		// z.S and x.S used but not changed, z not changed
		{
			if (ZTestSI(z))
				return 0;
			else if (z.S)
				return -1;
			else
				return 1;
		}
		
		public static int CompSI(MultiSI x, MultiSI z) // CompSI = sign(z - x)
		// sign = -1 if z < x sign = 0 if z = x sign = +1 if z > x
		// z.S and x.S used but not changed, z and x not changed
		{
			if (!z.S && !x.S)
				return AbCompUI(x, z);	// both positive
			else if (z.S && x.S)
				return AbCompUI(z, x);	// both negative
			else if (x.S)
				return 1;				// z pos, x neg
			else
			{
				return -1;				// z neg, x pos
			}
		}

		public static void WriteSIFLn(MultiSI z, ref int tot) // write z and base to outS at standard width with writeBool = false
		{
			bool save = writeBool;
			writeBool = false;
			WriteSILn(outS, null, z, ref tot, 0, true);
			writeBool = save;
		}

		public static void WriteSIF(MultiSI z, ref int tot) // write z and base to outS at standard width with writeBool = false
		{
			bool save = writeBool;
			writeBool = false;
			WriteSI(outS, null, z, ref tot, 0, true);
			writeBool = save;
		}

		public static void WriteSI(MultiSI z, ref int tot) // write z and base to outS at standard width
		{
			WriteSI(outS, null, z, ref tot, 0, true);
		}

		public static void WriteSI(MultiSI z, ref int tot, int width, bool WriteBase) // write z to outS
		{
			WriteSI(outS, null, z, ref tot, width, WriteBase);
		}

		public static void WriteSI(RichTextBox outB, StreamWriter outF, MultiSI z, ref int tot) // write z and base at standard width
		{
			WriteSI(outB, outF, z, ref tot, 0, true);
		}

		public static void WriteSILn(MultiSI z, ref int tot) // write z and base to outS at standard width and linefeed
		{
			WriteSILn(outS, null, z, ref tot, 0, true);
		}

		public static void WriteSILn(MultiSI z, ref int tot, int width, bool WriteBase) // write z to outS and linefeed
		{
			WriteSILn(outS, null, z, ref tot, width, WriteBase);
		}

		public static void WriteSILn(RichTextBox outB, StreamWriter outF, MultiSI z, ref int tot) // write z, base and linefeed
		{
			WriteSILn(outB, outF, z, ref tot, 0, true);
		}

		public static void WriteSILn(RichTextBox outB, StreamWriter outF, MultiSI z, ref int tot, int width, bool WriteBase)
		// output z and a new line
		{
			WriteSI(outB, outF, z, ref tot, width, WriteBase);
			WriteLn(outB, outF);
			tot = 0;
		}

		public static void WriteSI(RichTextBox outB, StreamWriter outF, MultiSI z, ref int tot, int width, bool WriteBase)
		// output z as a line of text, z not modified
		// good for any base = baseO [2, 36], base 10 much faster
		{
			int d;
			int g;
			int h;
			int i;
			int j;
			int k;
			int nOut;
			int toGo;
			int lenD;
			string st = "";
			string stL;
			int r;
			char ch;
			int maxDisL;
			if (d1LineOn)
				width = Int32.MaxValue;
			if (width == 0)
			{
				if (outB == null)
					width = 60;			// 60 for files only
				else
				{
					width = 80;			// 80 for screen
				}
			}
			if (z == 0 && dBoolOn)		// EQ1UI(0, z)
			{
				if (writeBool && outF == null) WriteT(outB, outF, "0 (False)", ref tot, width);
				else WriteT(outB, outF, "0", ref tot, width);
				if (outF != null && WriteBase && baseO != 10)
					WriteT(outB, outF, nL + "(base " + baseO + ")", ref tot, width);
				return;
			}
			else if (1 == z && dBoolOn)	// EQ1SI(1, z)
			{
				if (writeBool && outF == null) WriteT(outB, outF, "1 (True)", ref tot, width);
				else WriteT(outB, outF, "1", ref tot, width);
				if (outF != null && WriteBase && baseO != 10)
					WriteT(outB, outF, nL + "(base " + baseO + ")", ref tot, width);
				return;
			}
			maxDisL = maxDis;			// if to a file only, output all digits
			if (outB == null) maxDisL = Int32.MaxValue;
			stL = "";
			if (z.S) stL = "-";
			if (baseO != 10)
			{
				k = 'A' - 10;
				MultiSI t1 = new MultiSI(1);	// t1 is temp for z
				SetToSI(z, t1);
				t1.S = false;
				int chASize = (int)((z.N + 1) * muDMax / Math.Log10(baseO)) + 100;
				char[] chA = new char[chASize]; 
				i = 0;
				do
				{
					IpDiv1(baseO, out r, t1);	// t1 = t1 / baseO, r = Remainder
					if (r < 10)
						ch = r.ToString()[0];
					else
						ch = (char)(r + k);
					chA[i] = ch;
					i++;
					MuInterrupt();
					if (muAbort) break;			// exit do
					if (statReq)
						WriteStatus("Converting number to output base: " + (Math.Floor(100.0 - (100 * t1.N) / z.N)).ToString() + " Percent done");
				}
				while (t1.N != 1 || t1.V[0] != 0);
				DoneSI(t1);
				lenD = i;
				ch = chA[lenD - 1];			// first digit
				if (ch > '9') stL += "0";
				stL += ch;
				toGo = lenD - 1;
				nOut = 1;
				j = lenD;
				if (lenD > maxDisL) j = maxDisL / 2;
				i = lenD - 1;
				while (i != 0)
				{
					if (nOut >= j)
					{
						if (muDpG > 0)
						{
							if (toGo % muDpG == 0) stL += ",";
						}
						stL += "...";
						i = maxDis - j;
						toGo = i;
						j = Int32.MaxValue;
					}
					if (muDpG > 0)
						if (toGo % muDpG == 0) stL += ",";
					stL += chA[i - 1];
					if (stL.Length >= 120)
					{
						WriteT(outB, outF, stL, ref tot, width);	// don't let stL get too big
						stL = "";
						MuInterrupt();
						if (statReq && outB == null)
							WriteStatus("WriteSI: " + toGo + " to go of " + lenD + " digits");
						if (muAbort) break;		// exit do
					}
					nOut += 1;
					toGo -= 1;
					i -= 1;
				}
			}
			else // baseO = 10
			{
				toGo = z.N * muDMax;
				st = (baseM + z.V[z.N-1]).ToString();
				j = 1;
				do
					j++;
				while (st[j-1] == '0' && j != muDMax + 1);
				stL += st.Substring(j-1, 1);
				nOut = 1;
				toGo -= j - 1;
				lenD = toGo + 1;
				while (j <= muDMax && nOut < maxDisL)		// first few digits
				{
					j += 1;
					if (muDpG > 0)
					{
						if (toGo % muDpG == 0) stL += ",";
					}
					stL += st.Substring(j-1, 1);
					nOut += 1;
					toGo -= 1;
				}
				g = lenD;
				h = lenD;
				if (lenD > maxDisL)
				{
					g = maxDisL / 2;
					h = lenD - maxDisL + g;
				}
				for (i = 2; i <= z.N; i++)						// for each super digit
				{
					if (z.V[z.N - i] < 0)
					{
						st = "Negative number in a super digit";
						stL += "... " + st;
						if (outB != null)	// if output to screen
							st = "";
						Common.ReportError("WriteSI", st);
						break;
					}
					st = (baseM + z.V[z.N-i]).ToString();
					for (j = 2; j <= muDMax + 1; j++)			// for each decimal digit
					{
						if (nOut == g)
						{
							if (muDpG > 0)
							{
								if (toGo % muDpG == 0) stL += ",";
							}
							stL += "...";
							d = (h - g) / muDMax - 1;			// skip d super digits
							if (d > 0)
							{
								i += d;
								nOut += muDMax * d;
								toGo -= muDMax * d;
							}
						}
						if (nOut < g || nOut >= h)
						{
							if (muDpG > 0)
							{
								if (toGo % muDpG == 0) stL += ",";
							}
							stL += st.Substring(j-1, 1);
							if (stL.Length >= 120)
							{
								WriteT(outB, outF, stL, ref tot, width);	// don't let stL get too big
								stL = "";
								if (statReq && outB == null)
									WriteStatus("WriteSI: " + toGo + " to go of " + lenD + " digits");
							}
						}
						nOut += 1;
						toGo -= 1;
					}
					MuInterrupt();
					if (muAbort)
					{
						if (i < z.N) stL += "...";
						break;
					}
				}
			}
			WriteT(outB, outF, stL, ref tot, width);			// write the rest of it
			if ( (lenD >= muLenD || nOut < lenD || baseO != 10) && WriteBase )
			{
				// if output to a file only, put "(xxx Digits) on a new line
				if (outB == null)
				{
					Write(null, outF, nL);
					tot = 0;
					st = "";
				}
				else
				{
					st = " ";
				}
				if (baseO == 10 && dDigOn)
					WriteT(outB, outF, st + "(" + lenD + " Digits)", ref tot, width);
				if (baseO != 10)
				{
					if (lenD >= muLenD || lenD > maxDisL)
						WriteT(outB, outF, st + "(" + lenD + " Digits base " + baseO + ")", ref tot, width);
					else
					{
						WriteT(outB, outF, st + "(base " + baseO + ")", ref tot, width);
					}
				}
			}
			nOutSI = nOut;
		} // WriteSI

		public static void ShortWrSI(int shortD, MultiSI z, ref int tot)
		// output z in short form to outS
		{
			ShortWrSI(outS, null, shortD, z, ref tot);
		}

		public static void ShortWrSILn(RichTextBox outB, StreamWriter outF, int shortD, MultiSI z, ref int tot)
		// output z in short form and a new line
		{
			ShortWrSI(outB, outF, shortD, z, ref tot);
			WriteLn(outB, outF);
			tot = 0;
		}

		public static void ShortWrSI(RichTextBox outB, StreamWriter outF, int shortD, MultiSI z, ref int tot)
		// output z in short form ... if more than shortD digits
		// z is not modified}
		{
			int did;
			int j;
			int toGo;
			int lenD;
			string st;
			int nn;
			nn = z.N;
			if (nn < 3 || nn <= shortD / muDMax)
			{
				WriteSI(outB, outF, z, ref tot);
				return;
			}
			toGo = nn * muDMax;
			if (z.S) WriteT(outB, outF, "-", ref tot);
			st = (baseM + z.V[nn-1]).ToString();
			j = 1;
			while (st[j] == '0' && j != muDMax)
				j++;
			WriteT(outB, outF, st.Substring(j, 1), ref tot);
			did = 1;
			toGo -= j;
			lenD = toGo + 1;
			while (j < muDMax && did < 5)
			{
				j++;
				WriteT(outB, outF, st.Substring(j, 1), ref tot);
				did++;
			}
			st = (baseM + z.V[nn-2]).ToString();
			for (j = 1; j < 6 - did; j++)
				WriteT(outB, outF, st.Substring(j, 1), ref tot);
			WriteT(outB, outF, ",...,", ref tot);
			st = (baseM + z.V[0]).ToString();
			for (j = muDMax - 3; j <= muDMax + 1; j++)
				WriteT(outB, outF, st.Substring(j-1, 1), ref tot);
			WriteT(outB, outF, " (" + lenD + ")", ref tot);
		} // ShortWrSI

		public static void RMul1SI(int d1, MultiSI z) // z = z * d1
		// Multi-precision one super digit multiply d1 is not changed
		// muErr set true if z overflows
		{
			if (d1 < baseM)
			{
				int k;
				z.S ^= (d1 < 0);		// exclusive or
				IpMul1(d1, z, out k);
				if (k > 0)
					MuWriteErr("One digit multiply overflow, continuing... (M=" + z.M + ")");
			}
			else
			{
				MultiSI d = new MultiSI(sdInInt);
				d &= d1;
				MulSI(z, d, z); // z = z * d
			}
		}

		public static void Mul1SI(MultiSI x, int d1, MultiSI z) // z = x * d1
		// Multi-precision one super digit multiply x and d1 are not changed
		// muErr set true if z overflows
		{
			SetToSI(x, z);
			RMul1SI(d1, z);
		}

		public static void MulSI(MultiSI x, MultiSI y, MultiSI z) // z = x * y
		// x or y or both may have the same address in memory as z
		// x, y are not changed unless they share memory with z
		// sets muErr true if overflow
		{
#if SPEEDUP
			bool slow;
			double t1, t2;
			int n, ldn;
			if (!allowFHT)
				slow = true;
			else if (forceFHT)
				slow = false;
			else
			{
				int xN = x.N;
				int yN = y.N;
				t1 = xN * (double)yN;						// t1 proportional to runtime of normal multiply
				if (t1 < 3387.25 || xN <= 2 || yN <= 2)		// 3387.25 = 58.2^2
					slow = true;
				else
				{
					XIMult.FHTnLdn(xN, yN, out n, out ldn); // compute n for FHT
					n /= 2;
					t2 = 9.0 * n * (ldn - 1);				// t2 proportional to runtime of FHT multiply
					// the 9.0 above is a heuristic that causes a crossover at about xN = 58.2 if xN = yN
					// or about 232.5 decimal digits
					slow = (t1 < t2);						// t1 smaller so do normal multiply
				}
			}
			if (slow)
				MulSlowSI(x, y, z);
			else
			{
				z.S = x.S ^ y.S;
				XIMult.XIMul(x, y, z); // do a FHT multiply
			}
#else
			MulSlowSI(x, y, z);
#endif // SPEEDUP
		} // MulSI

		public static void MulSlowSI(MultiSI x, MultiSI y, MultiSI z) // z = x * y
		// x or y or both may have the same address in memory as z
		// x, y are not changed unless they share memory with z
		// sets muErr true if overflow
		{
			double timeAtStart = TimerF();
			z.S = x.S ^ y.S;								// exclusive or
			int xN = x.N;
			int yN = y.N;
			if (!Object.ReferenceEquals(z, x) && !Object.ReferenceEquals(z, y))
			{
				if (x.N >= y.N)
				{
					IpMul(x, y, z);			// 1 of 7 paths to be tested
					goto EndMulSI;
				}
				else
				{
					IpMul(y, x, z);			// 2
					goto EndMulSI;
				}
			}
			MultiUI zL = new MultiUI(z.N);
			SetToUI(z, zL);
			if (Object.ReferenceEquals(z, x) && Object.ReferenceEquals(z, y))
				IpMul(zL, zL, z);			// 3
			else if (Object.ReferenceEquals(x, z))
			{
				if (zL.N >= y.N)
					IpMul(zL, y, z);		// 4
				else
				{
					IpMul(y, zL, z);		// 5
				}
			}
			else // Object.ReferenceEquals(y, z)
			{
				if (zL.N >= x.N)
					IpMul(zL, x, z);		// 6
				else
				{
					IpMul(x, zL, z);		// 7
				}
			}
		EndMulSI:
			if (z.N == 1 && z.V[0] == 0) z.S = false; // prevent -0
		} // MulSlowSI

		public static void RMulSI(MultiSI x, MultiSI z) // z = z * x
		// x may have the same address in memory as z
		// x is not changed unless it shares memory with z
		// sets muErr true if overflow
		{
			MulSI(z, x, z);
		}

		public static void SqSI(MultiSI x, MultiSI z) // z = x * x
		// x may have the same address in memory as z
		// x is not changed unless it shares memory with z
		// sets muErr true if overflow
		{
			MulSI(x, x, z);
		}

		public static void RSqSI(MultiSI z) // z = z * z
		// sets muErr true if overflow
		{
			MulSI(z, z, z);
		}

		public static void ModMBSI(MultiSI z) // z = z % fMB
		// z is changed
		{
			ModSI(z, fMB, z);
		}

		public static void ModMBSI(MultiSI x, MultiSI z) // z = x % fMB
		// z is changed, sign(z) = sign(fMB)
		{
			ModSI(x, fMB, z);
		}

		public static void ModSI(MultiSI z, MultiSI mb) // z = z % mb
		// z is changed, sign(z) = sign(mb)
		{
			ModSI(z, mb, z);
		}

		public static void ModSI(MultiSI x, MultiSI mb, MultiSI z) // z = x % mb
		// Rem(x / mb), sign(z) = sign(fMB)
		// x or mb or both may have the same address in memory as z
		// x, mb are not changed unless they share memory with z
		/// z must be large enough to hold Rem, does not need to hold x
		// sets z = x if mb = 0
		{
			if (ZTestSI(mb))
			{
				SetToSI(x, z);
				return;
			}
			int i = AbCompUI(mb, x);
			if (i > 0)				// if |x| > |mb|
			{
				MultiSI qL = new MultiSI(x.N);
				if (!Object.ReferenceEquals(mb, z))
					IpDiv(x, mb, z, qL);
				else
				{
					MultiSI rL = new MultiSI(x.N);
					IpDiv(x, mb, rL, qL);
					SetToSI(rL, z);
					DoneSI(rL);
				}
				DoneSI(qL);
			}
			else if (i == 0)		// if |x| == |mb|
			{
				ClearSI(z);
			}
			else					// if |x| < |mb|
			{
				SetToSI(x, z);
				if (z.S != mb.S && z != 0)
					RAddSI(mb, z);
			}
		} // ModSI

		public static void RModSI(MultiSI mb, MultiSI z) // z = Rem(z / mb)
		// mb may have the same address in memory as z
		// mb is not changed unless it shares memory with z
		// no-op if mb = 0
		{
			ModSI(z, mb, z);
		}

		public static void RDiv1SI(int d1, MultiSI z) // z = z / d1
		// Multi-precision one super digit divide, d1 is not changed
		// allows literal d1, muErr set true if d1 = 0
		{
			int rL;
			IpDiv1(d1, out rL, z);
		}

		public static void Div1SI(MultiSI u, int d1, MultiSI z) // z = u / d1
		// Multi-precision one super digit divide, d1 and u not changed
		// u and z may have same address in memory, then u is changed
		// allows literal d1, muErr set true if d1 = 0
		{
			int rL;
			SetToSI(u, z);
			IpDiv1(d1, out rL, z);
		}

		public static void Mod1SI(int d1, out int r1, MultiSI z) // r1 = Rem(z / d1)
		// Multi-precision one super digit Modulo, d1 and z not changed
		// d1 and r1 may have same address in memory, then d1 is changed
		// allows literal d1, muErr set true if d1 = 0
		{
			MultiSI qL = new MultiSI(z.N);
			SetToSI(z, qL);
			IpDiv1(d1, out r1, qL);
			DoneSI(qL);
		}

		public static void RDiv1QRSI(int d1, out int r1, MultiSI z)
		// z = z / d1, r1 = Rem
		// Multi-precision one super digit Divide q + r. d1 and u not changed
		// d1 and r1 may have same address in memory, then d1 is changed
		// allows literal d1, muErr set true if d1 = 0
		{
			IpDiv1(d1, out r1, z);
		}

		public static void Div1QRSI(MultiSI u, int d1, out int r1, MultiSI z)
		// z = u / d1, r1 = Rem
		// Multi-precision one super digit Divide q + r. d1 and u not changed
		// d1 and r1 may have same address in memory, then d1 is changed
		// allows literal d1, muErr set true if d1 = 0
		{
			SetToSI(u, z);
			IpDiv1(d1, out r1, z);
		}

		public static void DivSI(MultiSI u, MultiSI d, MultiSI z) // z = u / d
		// u or d or both may have the same address in memory as z
		// u, d are not changed unless they share memory with z
		// sets muErr true if d = 0
		{
			MultiSI rL = new MultiSI(u.N);
			SetToSI(u, rL);
			if (Object.ReferenceEquals(d, z))
			{
				MultiSI dL = new MultiSI(d.N);
				SetToSI(d, dL);
				IpDiv(rL, dL, rL, z);
				DoneSI(dL);
			}
			else
			{
				IpDiv(rL, d, rL, z);
			}
			DoneSI(rL);
		}

		public static void RDivSI(MultiSI d, MultiSI z) // z = z / d
		// d may have the same address in memory as z
		// d is not changed unless it shares memory with z
		// sets muErr true if d = 0
		{
			DivSI(z, d, z);
		}

		public static void DivQRSI(MultiSI u, MultiSI d, MultiSI r, MultiSI z)
		// z = u / d, r = Rem
		// u or d or both may have the same address in memory as z or r
		// u, d are not changed unless they share memory with z or r
		// z and r may not have the same location in memory. Sets muErr true if d = 0
		{
			if (Object.ReferenceEquals(z, r))
			{
				MuWriteErr("Cannot set same location to both quotient and remainder, continuing...");
				return;
			}
			MultiSI rL = new MultiSI(u.N);
			if (!Object.ReferenceEquals(d, z))
				IpDiv(u, d, rL, z);		// 1 of 2 paths to test
			else
			{
				MultiSI qL = new MultiSI(u.N);
				IpDiv(u, d, rL, qL);	// 2
				SetToSI(qL, z);
				DoneSI(qL);
			}
			SetToSI(rL, r);
			DoneSI(rL);
			if (MultiID.muErr)
			{
				SetToSI(MultiID.sZero, r);
				SetToSI(MultiID.sZero, z);
			}

		}

		public static void RDivQRSI(MultiSI d, MultiSI r, MultiSI z)
		// z = z / d, r = Rem
		// z may have the same address in memory as d or r
		// d is not changed unless it shares memory with z or r
		// sets muErr true if d = 0
		{
			DivQRSI(z, d, r, z);
		}

		public static MultiSI PowSI(MultiSI b, MultiSI p)
		// returns b^p
		// b, p are not changed
		// sets muErr true if overflow or error
		{
			MultiSI z = new MultiSI(mMax);
			PowMBSI(b, p, sZero, z);
			return z;
		}

		public static void PowSI(MultiSI b, MultiSI p, MultiSI z)
		// z = (b^p)
		// b or p or both may have the same address in memory as z
		// b, p are not changed unless they share memory with z
		// sets muErr true if overflow or error
		{
			PowMBSI(b, p, sZero, z);
		}

		public static MultiSI ModPow(MultiSI b, MultiSI p, MultiSI mb)
		// returns b ^ p % mb
		{
			int m = mMax;
			if (mb != 0) m = 2 * Math.Max(b.N, mb.N);
			MultiSI z = new MultiSI(m);
			PowMBSI(b, p, mb, z);
			return z;
		}

		public static void PowMBSI(MultiSI b, MultiSI p, MultiSI z)
		// z = (b^p) % fMB
		// b or p or both may have the same address in memory as z
		// b, p are not changed unless they share memory with z
		// sets muErr true if overflow or error
		{
			PowMBSI(b, p, fMB, z);
		}

		public static void PowMBSI(int b, MultiSI p, MultiSI mb, MultiSI z)
		// z = (b^p) % mb
		// p may have the same address in memory as z
		// p are not changed unless it share memory with z
		// sets muErr true if overflow or error
		{
			MultiSI bL = new MultiSI(sdInInt);	// local value of b
			SetTo1SI(b, bL);
			PowMBSI(bL, p, mb, z);
			DoneSI(bL);
		}

		public static void Pow10SI(MultiSI p, MultiSI z)	// z = (10^p)
		// sets muErr true if overflow or error
		{
			long pL;
			long rem;
			if (p > MultiID.max64)
				MultiID.GetExSI(out pL, MultiID.max64);
			else
				MultiID.GetExSI(out pL, p);
			long k = Math.DivRem(pL, MultiID.muDMax, out rem);   // k = p/decimal digits in a super digit
			z &= (int)Math.Pow(10.0, (double)rem);
			int m = Math.Min(MultiID.mMax, (int)k + z.N);
			if (z.M < m)
			{
				z.M = m;
				Common.ReDimPreserve(ref z.V, m + 1);
			}
			XICommon.ShiftL(k, z);
		}

		public static void PowMBSI(MultiSI b, MultiSI p, MultiSI mb, MultiSI z)
		// z = (b^p) % mb
		// b or p or both may have the same address in memory as z
		// b, p are not changed unless they share memory with z
		// sets muErr true if overflow or error
		{

			if ((b == 10 || b == -10) && !p.S && ZTestSI(mb))
			{
				int lsbP = p.V[0] & 1;
				bool sign = (b.S && lsbP == 1);
				Pow10SI(p, z);
				if (muErr)
					ClearSI(z);
				else
					z.S = sign;
				return;
			}
			int r;
			bool fini;
			if (ZTestSI(b)) // b = 0
			{
				ClearSI(z);
				if (p.S)
					MuWriteErr("Cannot raise zero to a negative power");
				else if (ZTestSI(p))	// p = 0
					z.V[0] = 1;			// can raise zero to the zero power, 0^0 = 1
				return;
			}
			if ( p.S && (b.N > 1 || b.V[0] > 1) )		// (>1)^-p
			{
				ClearSI(z);
				return;
			}
			bool mbZero = (ZTestSI(mb));
			if (b.N == 1 && b.V[0] == 1) // b = +/-1
			{
				z.S = (b.S && (p.V[0] & 1) != 0); // Odd(iL)
				z.N = 1;
				z.V[0] = 1;
				if (!mbZero) ModSI(z, mb);
				return;
			}
			fini = false;
			MultiSI bL = new MultiSI(z.M);
			ModSI(b, mb, bL);
			MultiSI pL = new MultiSI(p.N);
			SetToSI(p, pL);
			SetTo1SI(1, z);
			do
			{
				if ((pL.V[0] & 1) != 0)			// Odd(iL)
				{
					if (z == sOne)
						z &= bL;
					else
						MulSI(bL, z, z);
					if (!mbZero) ModSI(z, mb);
				}
				IpDiv1(2, out r, pL);
				if (ZTestSI(pL))
					fini = true;
				else
				{
					MulSI(bL, bL, bL);
					if (!mbZero) ModSI(bL, mb);
				}
				if (statReq2)
				{
					int tot = 0;
					Common.SingleSpace();
					WriteT("Power = ", ref tot);
					WriteSI(pL, ref tot, 0, false);
					TimeDel(ref tot);
					RunForm.FocusBothUs();
					statReq2 = false;
				}
				if (muAbort)
				{
					ClearSI(z);
					fini = true;
				}
			}
			while (!fini);
			DoneSI(pL);
			DoneSI(bL);
		} // PowMBSI

		public static void RPowMBSI(MultiSI p, MultiSI z) // z = (z^p) % fMB
		// p may have the same address in memory as z
		// p is not changed unless it shares memory with z
		// sets muErr true if overflow or error
		{
			PowMBSI(z, p, z);
		}

		public static void LCMSI(MultiSI a, MultiSI b, MultiSI z)
		// z = Least common multiple a, b
		// a, b are not modified unless they share memory with z
		// will write to standard output if diagOn is true
		{
			MultiSI gcd = new MultiSI(Math.Min(a.N, b.N));
			GCDSI(a, b, gcd);
			if (gcd == 0)
				z &= 0;
			else
				z &= (a * b) / gcd;		// lcm(a, b) = (a * b) / gcd(a, b)
			z.S = false;
			DoneSI(gcd);
		} // LCMSI

		public static void GCDSI(MultiSI a, MultiSI b, MultiSI z)
		// z = Greatest common divisor a, b
		// a, b are not modified unless they share memory with z
		// will write to standard output if diagOn is true
		{
			int nn;
			bool fini = false;
			int tot = 0;
			MultiSI aL = new MultiSI(a.N);
			SetToSI(a, aL);
			MultiSI bL = new MultiSI(b.N);
			SetToSI(b, bL);
			aL.S = false;
			bL.S = false;
			if (ZTestSI(aL))
			{
				SetToSI(bL, z);
				fini = true;
			}
			if (ZTestSI(bL))
			{
				SetToSI(aL, z);
				fini = true;
			}
			nn = 0;
			if (debug && !fini) Common.SingleSpace();
			while (!fini)
			{
				nn++;
				if ((nn & 1) != 0)			// Odd(nn)
				{
					ModSI(aL, bL, aL);
					if (debug)
					{
						WriteT(nn + " ", ref tot);
						WriteSILn(aL, ref tot, 0, false);
					}
					if (ZTestSI(aL))
					{
						SetToSI(bL, z);
						fini = true;
					}
				}
				else
				{
					ModSI(bL, aL, bL);
					if (debug)
					{
						WriteT(nn + " ", ref tot);
						WriteSILn(bL, ref tot, 0, false);
					}
					if (ZTestSI(bL))
					{
						SetToSI(aL, z);
						fini = true;
					}
				}
			}
			DoneSI(bL);
			DoneSI(aL);
		} // GCDSI

		public static void RGCDSI(MultiSI a, MultiSI z)
		// z = Greatest common divisor a, z
		// a is not modified unless it shares memory with z
		// will write to standard output if diagOn is true
		{
			GCDSI(a, z, z);
		}

		public static void RLCMSI(MultiSI a, MultiSI z)
		// z = Least common multiple a, z
		// a is not modified unless it shares memory with z
		// will write to standard output if diagOn is true
		{
			LCMSI(a, z, z);
		}

		// compute x! = 1*2*3*...*x by binary splitting
		public static long maxK, currK;
		public static double oldPct;
		public static double timeAtStart;

		public static void ReportPct() // report status of binary splitting
		{
			if (!statReq2 && !diagOn) return;
			double pct;
			pct = (100.0 * currK) / maxK;
			if (pct > 100.0) pct = 100.0;
			if (pct >= oldPct + 5.0)
			{
				double deltaTime = TimerF() - timeAtStart;
				WriteLn(string.Format("Binary Split: {0:0.000}% complete, DT = {1:0.00} sec.  ", pct, deltaTime));
				statReq2 = false;
				oldPct = pct;
			}
		} // ReportPct

		public static void BinoSI(MultiSI a, MultiSI b, MultiSI z)
		// z = Binomial coefficient (a, b) by binary splitting
		// a and/or b may have the same address in memory as z
		// a and b are not changed unless they shares memory with z
		// sets muErr true if overflow or error
		{
			timeAtStart = TimerF();
			if (b.S)							// 1: (a, b) = 0 for all a if b < 0
			{
				ClearSI(z);
				return;
			}
			MultiSI aL = new MultiSI(a.N);
			SetToSI(a, aL);
			bool sign = false;
			if (aL.S)							// 2: if a < 0
			{
				aL &= (b - a - 1);				// (a, b) = (-1)^b * (b-a-1, b)
				sign = ((b.V[0] & 1) == 1);
			}
			if (CompSI(aL, b) > 0)		// 3: b > a
			{
				ClearSI(z);				// (a, b) = 0
				DoneSI(aL);
				return;
			}
			MultiSI bL = new MultiSI(a.N);
			SetToSI(b, bL);
			MultiSI v = new MultiSI(mMax);
			v &= aL / 2;
			if (CompSI(v, bL) > 0)		// 4: b > a / 2
			{
				bL &= aL - bL;					// b = a - b
			}
			if (ZTestSI(bL))			// 5: (a, 0) = 1 for all a
			{
				SetToSI(sOne, z);
				z.S = sign;
				DoneSI(v);
				DoneSI(bL);
				DoneSI(aL);
				return;
			}
			MultiSI nnSI = new MultiSI(bL.N);
			nnSI &= bL;
			double nn;
			GetDbSI(out nn, nnSI);
			bL &= aL - bL;					// b = a - b
			double maxD = Math.Round(nn * (Math.Log(nn, 2) + 1));
			if (maxD > Int64.MaxValue) maxD = Int64.MaxValue;
			maxK = (long)maxD;
			currK = 0;
			oldPct = 0.0;
			z &= PMr(aL, bL);	// z = Fac(a) / Fac(b)
			currK = maxK;
			if(oldPct < 100.0) ReportPct();
			currK = 0;	// maxK is the save
			oldPct = 0.0;
			v &= PMr(nnSI, sZero);	// v = Fac(a - b)
			currK = maxK;
			if(oldPct < 100.0) ReportPct();
			remSI =	new MultiSI(v.N);
			MultiID.DivQRSI(z, v, remSI, z);
			z.S = sign;
			DoneSI(v);
			DoneSI(bL);
			DoneSI(aL);
		} // BinoSI

		public static void Fac2SI(MultiSI a, MultiSI b, MultiSI z)
		// z = Fac(a) / Fac(b) % fMB by binary splitting
		// a and/or b may have the same address in memory as z
		// a and b are not changed unless they shares memory with z
		// sets muErr true if overflow or error
		{
			timeAtStart = TimerF();
			if (a.S || b.S)
			{
				MuWriteErr("Cannot take factorial of number < zero");
				ClearSI(z);
				return;
			}
			if (b > a)
			{
				ClearSI(z);
				return;
			}
			MultiSI nnSI = new MultiSI(a.N);
			nnSI &= a - b;
			double nn;
			GetDbSI(out nn, nnSI);
			double maxD = Math.Round(nn * (Math.Log(nn, 2) + 1));
			if (maxD > Int64.MaxValue) maxD = Int64.MaxValue;
			maxK = (long)maxD;
			currK = 0;
			oldPct = 0.0;
			z &= PMr(a, b);		// z = Fac(a) / Fac(b)
			currK = maxK;
			if(oldPct < 100.0) ReportPct();
		} // Fac2SI

		public static void FacMBiSI(MultiSI x, MultiSI z)
		// z = x! = 1*2*3*...*x % fMB by binary splitting
		// x may have the same address in memory as z
		// x is not changed unless it shares memory with z
		// sets muErr true if overflow or error
		{
			Fac2SI(x, sZero, z);
		}

		public static MultiSI PMr(MultiSI a, MultiSI b)
		// computes a!/b! % fMB by binary splitting (a >= b)
		{
			if (a == b) return sOne;
			if (a < b || muErr || muAbort) return sZero;
			MultiSI dSI = new MultiSI(a.N);
			dSI &= a - b;
			long d;
			GetExSI(out d, dSI);
			if (muErr)			// d is too large
			{
				MuWriteErr("Number too large for factorial function");
				return sZero;
			}
			MultiSI p = new MultiSI();
			bool notZero = !ZTestSI(fMB);
			if (d > 3)
			{
				MultiSI m = new MultiSI(a.N);
				m &= (a + b) / 2;
				p &= PMr(a, m) * PMr(m, b);
				if (notZero) ModMBSI(p);
				if (currK > (maxK - d))	// currK += d could overflow
					currK = maxK;
				else
					currK += d;
			}
			else
			{
				p &= a;
				if (notZero) ModMBSI(p);
				if (currK > (maxK - 1)) currK = maxK;
				else currK += 1;
				if (d > 1)
				{
					p *= (a-1);
					if (notZero) ModMBSI(p);
					if (currK > (maxK - 3)) currK = maxK;
					else currK += 3;
				}
				if (d > 2)
				{
					p *= (a-2);
					if (notZero) ModMBSI(p);
					if (currK > (maxK - 4)) currK = maxK;
					else currK += 4;
				}
			}
			ReportPct();
			return p;
		} // PMr

		public static void FacMBSI(MultiSI x, MultiSI z)
		// z = (x!) % fMB = 1*2*3*...*x
		// x may have the same address in memory as z
		// x is not changed unless it shares memory with z
		// sets muErr true if overflow or error
		{
			double nn;
			if (x.S)
			{
				MuWriteErr("Cannot take factorial of number < zero");
				ClearSI(z);
				return;
			}
			if (x.N <= 2)
				GetDbSI(out nn, x);
			else
				nn = (double)(Int32.MaxValue) + 1;
			if (nn > Int32.MaxValue) // an arbitrary limit on n!
			{
				MuWriteErr("Number too large for factorial function");
				ClearSI(z);
				return;
			}
			MultiSI xL = new MultiSI(sdInInt);
			ReSizeSI(Calc.fMC, z);
			SetTo1SI(1, z);
			nDone = 0;
			while (nn > 1)
			{
				SetToDbSI(nn, xL);
				MulSI(xL, z, z);
				ModMBSI(z);
				if (muAbort)
				{
					WriteLn("Factorial function aborted by operator, " + nn +
						" multiplies to go");
					ClearSI(z);
					nn = 0;
				}
				nn--;
				nToGo = (long)nn;
				nDone++;
				if (statReq2) WritenDone("FacMBSI:");
			}
			DoneSI(xL);
			nn = z.N;
			Common.ReDimPreserve(ref z.V, (int)nn+1); // preserve and re-dimention arrays
			z.M = (int)nn;
			nToGo = 0;
		} // FacMBSI

		public static void RFacMBSI(MultiSI z) // z = (z!) % fMB = 1*2*3*...*x
		// sets muErr true if overflow or error
		{
			FacMBSI(z, z);
		}

		public static void SqRtRemSlowSI(MultiSI x, MultiSI r, MultiSI z)
		// z = SqRt(x), r = Rem
		// x is not changed, x is moved to r and then r is "divided"
		// in place to give the SqRt in z and a remainder in r, it is ok for x and r
		// or x and z to have the same location in memory, but then x will be changed
		// z and r may not have the same location in memory. Sets muErr true if error
		{
			int d0;
			int qH;
			double tD;
			double d4;
			long t;
			long t2;
			long rem;
			long k;
			int i;
			int j;
			int rI0;
			int dI0;
			int rI;
			int dI;
			int iJ;
			if (muErr || muAbort)
			{
				ClearSI(r);
				ClearSI(z);
				return;
			}
			if (Object.ReferenceEquals(z, r))
			{
				MuWriteErr("Cannot set same location to both square root and remainder, continuing...");
				return;
			}
			SetToSI(x, r);
			ClearSI(z);
			if (r.S)
			{
				MuWriteErr("Cannot take square root of negative number, continuing...");
				ClearSI(r);
				return;
			}
			if (r.N < 4)
			{
				for (i = r.N; i <= 3; i++)
					r.V[i] = 0;
				r.N = 4;
			}
			if ((r.N & 1) != 0)			// Odd(r.N)
			{
				r.V[r.N] = 0;
				r.N++;
			}
			z.N = r.N / 2 + 1;
			tD = 0;
			for (i = 1; i <= 4; i++)
				tD = tD * baseM + r.V[r.N-i];
			t2 = (long)Math.Floor(Math.Sqrt(tD));
			t2 += 3; // HJS 2008-08-23, this fixed a problem with
			// SqRt(91258893,03356221,57253936,60000000,00000000,00000000,00000000,00000000)
			do
			{
				d0 = (int)Math.DivRem(t2, baseML, out rem);
				qH = 1 + (int)rem;						// +1 for safety
				t = Math.BigMul(r.V[r.N-1], baseM) + r.V[r.N-2] - Math.BigMul(d0, d0);
				if (t < 0) t2--;
			}
			while (t < 0);
			k = Math.DivRem(t, baseML, out rem);
			r.V[r.N-2] = (int)rem;
			r.V[r.N-1] = (int)k;
			if (z.N >= 4)
				z.V[z.N - 4] = 0;
			z.V[z.N - 3] = qH;
			t = 2 * d0;
			k = Math.DivRem(t, baseML, out rem);
			z.V[z.N-2] = (int)(rem);
			z.V[z.N-1] = (int)k;
			statReq = false;
			for (j = 2; j <= r.N / 2; j++)
			{
				k = baseM;
				rI0 = r.N - j - j;
				dI0 = z.N - j - 1;
				rI = rI0;
				dI = dI0;
				iJ = rI0 + j;
				while (rI <= iJ)			// multiply and subtract
				{
					t = r.V[rI] - Math.BigMul(qH, z.V[dI]) + k + bSMB; // bSMB = baseSq - baseM
					k = Math.DivRem(t, baseML, out rem);
					t = rem;
					r.V[rI] = (int)t;
					rI++;
					dI++;
				}
				t = r.V[rI] + k + bSMB;
				k = Math.DivRem(t, baseML, out rem);
				r.V[rI] = (int)rem;
				while (k == maxDigit)		// test remainder
				{
					qH--; // add back
					k = qH;
					rI = rI0;
					dI = dI0;
					while (rI <= iJ)
					{
						t = r.V[rI] + z.V[dI] + k;
						k = Math.DivRem(t, baseML, out rem);
						r.V[rI] = (int)rem;
						rI++;
						dI++;
					}
					t = r.V[rI] + k;
					k = Math.DivRem(t, baseML, out rem);
					r.V[rI] = (int)rem;
					k += maxDigit;
					z.V[dI0] = qH;
				}
				i = dI0;
				iJ = i + j;
				k = qH;
				while (i <= iJ)
				{
					t = z.V[i] + k;			// carry to next super digit
					k = Math.DivRem(t, baseML, out rem);
					z.V[i] = (int)rem;
					if (k == 0)
						i = iJ;				// break out of loop
					i++;
				}
				if (j < r.N / 2)
				{
					tD = 0;
					for (i = 0; i <= 4; i++)
						tD = tD * baseM + r.V[r.N - j - i];
					d4 = 0;
					for (i = 1; i <= 4; i++)
						d4 = d4 * baseM + z.V[z.N-i];
					qH = 1 + (int)Math.Floor(tD / d4);
					z.V[z.N - j - 2] = qH;
				}
				MuInterrupt();
				if (muAbort)
				{
					ClearSI(z);
					ClearSI(r);
					return;
				}
				if (statReq)
					WriteStatus("Integer Square root: " + Math.Floor(100.0 * j / (r.N / 2)).ToString() + " Percent done");
			}
			r.N = r.N / 2 + 1;
			NormSI(r);
			NormSI(z);
			IpDiv1(2, out qH, z);
		} // SqRtRemSlowSI

		public static void SqRtRemFastSI(MultiSI x, MultiSI r, MultiSI z)
		// z = SqRt(x), r = Rem, Fast method for large numbers
		// x is not changed, x is moved to r and then r is "divided"
		// in place to give the SqRt in z and a remainder in r, it is ok for x and r
		// or x and z to have the same location in memory, but then x will be changed
		// z and r may not have the same location in memory. Sets muErr true if error
		{
			MultiSI q = new MultiSI();		// q = xL / yL
			MultiSI d = new MultiSI();		// rL / t
			MultiSI tD = new MultiSI();		// t * d
			int i, j, pl, po, ps, xs, ys, qs; // scalings, x = xL * Base^xs
			bool adjusted;	// true iff an adjustment was needed

			pl = 5;
			j = pl+pl+2;
			if ((x.N & 1) != 0)  j--;	// j = pl+pl+1
			if (x.N < j || x.S)
			{
				SqRtRemSlowSI(x, r, z);
				return;
			}
			MultiSI xL = new MultiSI(j); // xL is x at a reduced precision
			ps = x.N - j;
			ReducePrecision(ps, x, xL);
			MultiSI yL = new MultiSI(pl+2);		// yL is Local Y = SqRt(x)
			MultiSI rL = new MultiSI(j);		// rL is Local remainder
			SqRtRemSlowSI(xL, rL, yL);
			DoneSI(rL);
			DoneSI(xL);
			ys = (x.N - j) / 2;
			while (pl <= (x.N / 2)) { // y = (y + x / y) / 2
				po = pl;
				pl += pl;
				if (pl > (x.N / 2))  pl = 1 + x.N / 2;
				j = pl+po+2;
				if (j > x.N)  j = x.N;
				ps = x.N-j;
				ReducePrecision(ps, x, xL);	// xL is x at a reduced precision

				xs = x.N - xL.N;
				ReSizeSI(1 + xL.N - yL.N, q);	// q = xL / yL
				ReSizeSI(xL.N, rL);
				IpDiv(xL, yL, rL, q);
				DoneSI(xL);
				DoneSI(rL);
				qs = xs - ys;
				if (ys > qs)
				{
					j = yL.N - 1;
					i = j + ys - qs;
					Common.ReDimPreserve(ref yL.V, i+3);	// preserve and re-dimention arrays
					yL.N = yL.M = i+2;
					for ( ;  i >= 0;  i--, j--)
					{
						if (j >= 0)
							yL.V[i] = yL.V[j];
						else
							yL.V[i] = 0;
					}
					ys = qs;
				}
				if (ys < qs)
					MuWriteErr("ys < qs in SqRtRemFastSI()");
				RAddSI(q, yL);
				DoneSI(q);
				RDiv1SI(2, yL);
				if (ys < 0)
				{
					for (i = 0, j = -ys, yL.N -= j;  i < yL.N;  i++, j++)
						yL.V[i] = yL.V[j];
					ys = 0;
				}
			}
			ReSizeSI(yL.N + yL.N, xL);
			MulSI(yL, yL, xL);
			ReSizeSI(x.N+1, rL);
			SubSI(x, xL, rL);
			DoneSI(xL);
			MultiSI t = new MultiSI(yL.N+1);	// t = 2 * yL
			AddSI(yL, yL, t);
			adjusted = false;
			while (rL.S || (rL > t))			// was CompSI(t, rL) > 0)
			{
				adjusted = true;
				if (rL.S)
				{
					if (t.N > rL.N)
					{
						DoneSI(d);
					}
					else
					{
						rL.S = false;	// plus
						ReSizeSI(1 + rL.N - t.N, d);
						DivSI(rL, t, d);
						rL.S = true;	// minus;
					}
					RAdd1SI(1, d);
					RSubSI(d, t);
					ReSizeSI(t.N + d.N, tD);
					MulSI(t, d, tD);
					RAddSI(tD, rL);
					DoneSI(tD);
					RSubSI(d, t);
				}
				else if (rL > t)	// (CompSI(t, rL) > 0)
				{
					ReSizeSI(1 + rL.N - t.N, d);
					DivSI(rL, t, d);
					RAddSI(d, t);
					ReSizeSI(t.N + d.N, tD);
					MulSI(t, d, tD);
					RSubSI(tD, rL);
					DoneSI(tD);
					RAddSI(d, t);
				}
			}
			if (adjusted)  Div1SI(t, 2, yL);
			DoneSI(t);
			SetToSI(rL, r);		// set r to rL (copy)
			SetToSI(yL, z);		// set z to yL (copy)
			DoneSI(yL);
		} // SqRtRemFastSI

		public static void SqRtRemSI(MultiSI x, MultiSI r, MultiSI z)
		// z = SqRt(x), r = Rem
		// x is not changed, x is moved to r and then r is "divided"
		// in place to give the SqRt in z and a remainder in r, it is ok for x and r
		// or x and z to have the same location in memory, but then x will be changed
		// z and r may not have the same location in memory. Sets muErr true if error
		{
			if ((x.N <= 2024 && !forceFHT) || !allowFHT)
				SqRtRemSlowSI(x, r, z);
			else
				SqRtRemFastSI(x, r, z);
		}

		public static void RSqRtRemSI(MultiSI r, MultiSI z) // z = SqRt(z), r = Rem
		// z is moved to r and then r is "divided" in place to give the SqRt in z
		// and a remainder in r, z and r may not have the same location in memory
		// sets muErr true if error
		{
			SqRtRemSI(z, r, z);
		}

		public static void SqRootSI(MultiSI x, MultiSI z) // z = SqRt(x)
		// x is not changed, it is ok for x and z to have the same location in memory
		// but then x will be changed. Sets muErr true if error
		{
			MultiSI rL = new MultiSI(Math.Max(x.N, 4));
			SqRtRemSI(x, rL, z);
			DoneSI(rL);
		}

		public static void RSqRtSI(MultiSI z) // z = SqRt(z)
		// sets muErr true if error
		{
			SqRootSI(z, z);
		}

		public static void GetDbInSI(out double db, double low, double high, MultiSI z)
		// db = z % baseM^3, in a range, + or -
		{
			GetDbSI(out db, z);
			if (db < low) db = low;
			if (db > high) db = high;
		}

		public static void Write(string message) // write a short message to outS
		{
			Write(outS, null, message);
		}

		public static void WriteLn() // write a linefeed to outS
		{
			WriteLn(outS, null, "");
		}

		public static void WriteLn(string message) // write a short message and a linefeed to outS
		{
			WriteLn(outS, null, message);
		}

		public static void WriteLn(RichTextBox outB, StreamWriter outF) // Write linefeed
		{
			Write(outB, outF, nL);
		}

		public static void WriteLn(RichTextBox outB, StreamWriter outF, string st) // Write st plus linefeed
		{
			Write(outB, outF, st + nL);
		}

		public static void WriteLn2() // Write 2 linefeeds to outS
		{
			Write(outS, null, nL2);
		}

		public static void WriteLn2(string st) // Write st plus 2 linefeeds to outS
		{
			Write(outS, null, st + nL2);
		}

		public static void WriteLn2(RichTextBox outB, StreamWriter outF) // Write 2 linefeeds
		{
			Write(outB, outF, nL2);
		}

		public static void WriteLn2(RichTextBox outB, StreamWriter outF, string st) // Write st plus 2 linefeeds
		{
			Write(outB, outF, st + nL2);
		}

		public static void WriteSLn(string st) // Write SingleSpace, st, and linefeed to outS
		{
				Common.SingleSpace(outS);
				Write(outS, null, st + nL);
		}

		public static void WriteDLn(string st) // Write DoubleSpace, st, and linefeed to outS
		{
				Common.DoubleSpace(outS);
				Write(outS, null, st + nL);
		}

		public static void WriteDLn2(string st) // Write DoubleSpace, st, and 2 linefeeds to outS
		{
				Common.DoubleSpace(outS);
				Write(outS, null, st + nL2);
		}

		public static void WriteDLn(RichTextBox outB, StreamWriter outF, string st) // Write DoubleSpace, st, and linefeed
		{
				Common.DoubleSpace(outB);
				Write(outB, outF, st + nL);
		}

		public static void Write(RichTextBox outB, StreamWriter outF, string st) // print string to RichTextBox or to a file
		{
			if (st == "") return;
			if (outF != null)
			{
				try	{outF.Write(st);}
				catch (Exception ex) {Common.ReportError("Write", ex.Message);}
			}
			if (outB != null)
			{
				outB.AppendText(st);
				if (st[(st.Length - 1)] == '\n')	// if a new line, move to bottom/end of output text box
				{
					Global.runTextOut.Select(Global.runTextOut.Text.Length, 0);
					Global.runTextOut.Focus();
				}
				if (logg && (outB == outS))			// only log output to Run form textbox
				{
					try	{outL.Write(st);}
					catch {}
				}
				if (outB.Text.Length > 200000)		// clear top of screen after about 50 full pages
				{
					int i = outB.Text.IndexOf("\n", 4000);
					outB.Text = outB.Text.Substring(i + 1);	// delete partial line
				}
			}

		} // Write

		public static void WriteT(string st)
		// WriteT with 0 tot and 0 width to outS
		{
			int tot = 0;
			Common.SingleSpace();
			WriteT(outS, null, st, ref tot, 0);
		}

		public static void WriteT(string st, ref int tot)
		// WriteT with 0 width to outS
		{
			WriteT(outS, null, st, ref tot, 0);
		}

		public static void WriteT(RichTextBox outB, StreamWriter outF, string st)
		// WriteT with 0 tot and 0 width
		{
			int tot = 0;
			Common.SingleSpace(outB);
			WriteT(outB, outF, st, ref tot, 0);
		}

		public static void WriteT(RichTextBox outB, StreamWriter outF, string st, ref int tot)
		// WriteT with 0 width
		{
			WriteT(outB, outF, st, ref tot, 0);
		}

		public static void WriteTLn()
		// WriteTLn with "", 0 tot and 0 width to outS
		{
			Write(outS, null, nL);
		}

		public static void WriteTLn(string st)
		// WriteTLn with 0 tot and 0 width to outS
		{
			int tot = 0;
			Common.SingleSpace();
			if (st != "") WriteT(outS, null, st, ref tot, 0);
			Write(outS, null, nL);
		}

		public static void WriteTLn(string st, ref int tot)
		// WriteT and linefeed to outS
		{
			if (st != "") WriteT(outS, null, st, ref tot, 0);
			Write(outS, null, nL);
			tot = 0;
		}

		public static void WriteTLn(RichTextBox outB, StreamWriter outF, string st)
		// WriteT and linefeed with 0 tot and 0 width
		{
			int tot = 0;
			Common.SingleSpace(outB);
			if (st != "") WriteT(outB, outF, st, ref tot, 0);
			Write(outB, outF, nL);
		}

		public static void WriteTLn(RichTextBox outB, StreamWriter outF, string st, ref int tot)
		// WriteT and linefeed with 0 width
		{
			if (st != "") WriteT(outB, outF, st, ref tot, 0);
			Write(outB, outF, nL);
			tot = 0;
		}

		public static void WriteTLn(RichTextBox outB, StreamWriter outF, string st, ref int tot, int width)
		// WriteT and linefeed
		{
			if (st != "") WriteT(outB, outF, st, ref tot, width);
			Write(outB, outF, nL);
			tot = 0;
		}

		public static void WriteTLn2(string st)
		// WriteTLn2 with 0 tot and 0 width to outS
		{
			int tot = 0;
			Common.SingleSpace();
			if (st != "") WriteT(outS, null, st, ref tot, 0);
			Write(outS, null, nL2);
		}

		public static void WriteTLn2(RichTextBox outB, StreamWriter outF, string st)
		// WriteTLn2 with 0 tot and 0 width
		{
			int tot = 0;
			Common.SingleSpace(outB);
			if (st != "") WriteT(outB, outF, st, ref tot, 0);
			Write(outB, outF, nL2);
		}

		public static void WriteTLn2(RichTextBox outB, StreamWriter outF, string st, ref int tot, int width)
		// WriteT and a double linefeed after
		{
			if (st != "") WriteT(outB, outF, st, ref tot, width);
			Write(outB, outF, nL2);
			tot = 0;
		}

		public static void WriteTD(RichTextBox outB, StreamWriter outF, string st)
		// force double space and WriteT with 0 tot and 0 width
		{
			int tot = 0;
			Common.SingleSpace(outB);
			WriteT(outB, outF, st, ref tot, 0);
		}

		public static void WriteTDLn(string st)
		// WriteT with a guaranteed double linefeed before and single linefeed after to outS
		{
			Common.DoubleSpace(outS);
			WriteTLn(outS, null, st);
		}

		public static void WriteTDLn(RichTextBox outB, StreamWriter outF, string st)
		// WriteT with a guaranteed double linefeed before and single linefeed after
		{
			Common.DoubleSpace(outB);
			WriteTLn(outB, outF, st);
		}

		public static void WriteT(RichTextBox outB, StreamWriter outF, string st, ref int tot, int width)
		// output String and linefeed every 60, 80 or width characters
		// tot = characters on line so far
		{
			int i;
			int j;
			int iL;
			int n;
			string stE;
			if (width > 0)
				n = width;
			else
			{
				if (outB == null) n = 60;	// 80 for screen, 60 for files
				else n = 80;
			}
			i = 1;
			j = n - tot;
			iL = st.Length;
			int ix = 0;
			while (iL > 0)
			{
				stE = "";
				if (tot >= n)
				{
					stE = nLc;				// format n characters per line
					j = n;
					tot = 0;
				}
				if (iL <= j)
				{
					Write(outB, outF, stE + st.Substring(i - 1));	// built for speed, don't call Write too often
					ix = st.Substring(i - 1).LastIndexOf("\n");
					tot += iL;
					if (ix >= 0)
						tot = st.Substring(i - 1).Length - ix - 1;
					iL = 0;
				}
				else
				{
					Write(outB, outF, stE + st.Substring(i-1, j));
					ix = st.Substring(i - 1, j).LastIndexOf("\n");
					tot = n;
					if (ix >= 0)
						tot = st.Substring(i - 1, j).Length - ix - 1;
					i += j;
					iL -= j;
					j = n;
				}
			}
		} // WriteT

		public static void MuWriteErr(string st) // write error statement to outS and set error flag
		{
			muErr = true;
			if (muErrRep) WriteDLn(st);
		}

		public static void DiagT(string message, ref int tot) // output a diagnostic message to outS
		{
			if (diagOn) WriteT(message, ref tot);
		}

		public static void ADiag(string message) // output a diagnostic message to outS
		{
			if (diagOn) WriteLn(message);
		}

		public static void ADebug(string message) // output a debug message to outS
		{
			if (debug) WriteLn(message);
		}

		public static void ADDiag(string message) // output DoubleSpace and diagnostic message to outS
		{
			if (diagOn) WriteDLn(message);
		}

		public static void WritenDone(string st) // write number nDone ...
		{
			int tot = 0;
			if (st != "") WriteT(st + " ", ref tot);
			WriteT("(Done " + nDone + " of " + (nDone + nToGo) + ", " + nToGo + " to go)", ref tot);
			TimeDel(ref tot);
			RunForm.FocusBothUs();
			statReq2 = false;
		}

		public static void WriteStatus(string st) // write status message
		{
			int tot = 0;
			Common.SingleSpace();
			WriteT(st, ref tot);
			TimeDel(ref tot);
			RunForm.FocusBothUs();
			statReq = false;
		}

		public static void MuInterrupt() // test for ESC key pressed to abort operation
		// muAbort set true if Abort selected, unchanged if not
		// softAbort set true if Retry selected
		{
			muICount++;
			if (muICount > 100 || (!intCalc && Calc.fMC > 100000))
			{
				muICount = 0;
				int Time = Environment.TickCount;	// time in milliseconds
				if (Time - muITimeSaved > 500 || Time < muITimeSaved)
				{
					// don't hog the processor, but don't slowdown too much
					Application.DoEvents();		// process user requested actions;
					muITimeSaved = Time;
				}
			}
			if (!abortCalc) return;
			if (nToGo > 0) WritenDone("");
			if (muAbort) return;
			double TimeSaved = TimerF();	// allow for time paused
			DialogResult MsgButton = MessageBox.Show("*** INTERRUPT: To Continue select Ignore\n" +
									"To Abort Computation select Abort\n" +
									"To Set softAbort flag select Retry\n",
				Global.title,
				MessageBoxButtons.AbortRetryIgnore,
				MessageBoxIcon.Information);
			// ----------
			double TimeOut = TimerF() - TimeSaved;	// allow for time-out
			time0 += TimeOut;
			abortCalc = false;
			if (MsgButton == DialogResult.Ignore)
				WriteDLn("Continuing...");
			if (MsgButton == DialogResult.Abort)
			{
				muAbort = true;
				WriteDLn("Computation aborted by operator!");
			}
			if (MsgButton == DialogResult.Retry)
			{
				softAbort = true;
				WriteDLn("softAbort flag set by operator!");
			}
			RunForm.FocusBothUs();
		} // MuInterrupt

		public static double TimerF() // get the value of the system clock in total seconds from start of run
		// good to one clock tick ~= 0.0115 seconds (0.055 on some systems)
		// must call once each 24 days to maintain time0
		{
			double time1 = Environment.TickCount/1000.0; // TickCount = number of milliseconds since Windows started.
			if (time1 < timeP)							 // TickCount cycles 2^31 mills = 24.855... days
				time0 -= Math.Pow(2, 31)/1000;			 // time0 = 2^53 milliseconds = 285,420.9... years (long enough for me)
			timeP = time1;
			return time1 - time0;
		}

		public static void TimeDel() // output delta time only
		{
			int tot = 0;
			TimeDel(ref tot);
		}

		public static void TimeDel(ref int tot) // output delta time only
		{
			GC.Collect();
			long totalMemory = GC.GetTotalMemory(false);
			if (timeIt)
			{
				string st = ", Memory in use = " + totalMemory + ", DT = ";
				double tV2L = TimerF();
				double delta = tV2L - tV2;	// delta time from last TimeIt call
				WriteTLn(st + delta.ToString("0.00") + " sec.", ref tot);
			}
			else
			{
				WriteTLn("", ref tot);
			}
		}

		public static void TimeIt(string message) // output a diagnostic message with delta times
		{
			if (timeIt)
			{
				tV1 = tV2;
				tV2 = TimerF();
				double delta = tV2 - tV1;	// delta time from last TimeIt call
				WriteSLn("T = " + tV2.ToString("0.00") +
					"  DT = " + delta.ToString("0.00") + " sec.  " + message);
			}
		}

		public static void Diag(string message) // output a diagnostic message with delta times
		{
			if (diagOn) TimeIt(message);
		}

		public static string MaxDig(int iBase)
		{
			if (iBase <= 10)
			{
				return Convert.ToString(iBase - 1);
			}
			else
			{
				return "" + (char)(iBase - 11 + (Int16)('A'));
			}
		}

		public static void Delete(ref string st, int index, int count)
		// delete count characters starting at index, count >= 0
		{
			st = st.Substring(0, index) + st.Substring(index + count);
		}

		// in Intel micro processors the max odd integer to store
		// in an Int16 is 2^15-1 = 32767 (= short)
		// in a float is 2^24-1 = 167,77215
		// in a Int32 is 2^31-1 = 21474,83647 (= int)
		// for a double is 2^53-1 = 9,00719,92547,40991
		// for Int64 is 2^63-1 = 9223,37203,68547,75807 (= long)
		// for an Extended or Long Double is 2^64-1 = 18446,74407,37095,51615
		// for a decimal is 2^96-1 = 7922,81625,14264,33759,35439,50335
		// c has Long Double, Pascal calls it Extended, but not available in Visual Basic or C#

		public static int muDMax = 8;		// max number of decimal digits in a super digit
		public static int muNMax = 8388650;	// max number of super digits in a multi-precision // can be changed in StartupForm
											// number. Allows for 67,109,200 = 2^26 + 336 decimal digits
											// each super digit is from 0 to 99,999,999, i.e. base = 10^8
		private static int muICount;		// belongs to MuInterrupt()
		private static int muITimeSaved;	// belongs to MuInterrupt()

		public static bool abortCalc;	// true if abort calc button clicked
		public static bool allowFHT;	// the FHT multiply is supported in this project and currently used
		public static bool d1LineOn;		// display all on one line to screen
		public static bool dBoolOn;		// output bool turned on
		public static bool dDigOn;		// output digits turned on
		public static bool debug;		// debug mode turned on
		public static bool diagOn;		// diagnostics turned on
		public static bool forceFHT;	// FHT is to be used for all multiplies
		public static bool logg;		// true if Logging screen output to FILE
		public static bool muAbort;		// true if abort calc confirmed by operator
		public static bool muErr;		// true if multi-precision error occurred
		public static bool muErrRep;	// true if multi-precision error reports are desired
		public static bool muUnderflow;	// true if multi-precision underflow occurred
		public static bool softAbort;	// true if softAbort Flag set by operator
		public static bool statReq;		// true if Status of calculation is requested, Level 1
		public static bool statReq2;	// true if Status of calculation is requested, Level 2
		public static bool statReq3;	// true if Status of calculation is requested, Level 3
		public static bool timeIt;		// timing mode turned on
		public static bool useFMB;		// true if fMB is to be used by operator ^
		public static bool writeBool;	// true if "0 (False)" and "1 (True)" to be output for 0 and 1
		public static bool intCalc;		// true iff Integer calculator
		public static decimal decMaxValue;	// max decimal holding an integer // 2^96-1 = 7922,81625,14264,33759,35439,50335
		public static double dblMaxValue;	// max double holding an integer // 2^53-1 = 9,00719,92547,40991
		public static double time0;		// time in seconds of start of run
		public static double timeP;		// previous value of system clock for time recycle test
		public static double tV1;		// time value for delta time, previous value of tV2
		public static double tV2;		// time value for total time, delta = tV2 - tV1
		public static int baseI;		// display base used for input [2, 35]
		public static int baseM;		// radix of multi-precision numbers, 10^muDMax
		public static int baseMr;		// square soot of radix baseM
		public static int baseO;		// display base used for output [2, 35]
		public static int maxDigit;		// baseM - 1 (=99999999 when muDMax=8, =9999 when muDMax=4)
		public static int maxDis;		// current max digits to display an integer, >= 2
		public static int mMax;			// current Max number of super digits <= muNMax
		public static int muDpG;		// digits per group in number output display >= 0
		public static int muLenD;		// min length in digits to display length in digits
		public static int nOutSI;		// nOutSI = nOut from Sub WriteSI, digits output
		public static int rem1;			// holds the remainder of infix operator SI / int
		public static int sdInDbl;		// Super digits in a double holding an integer
		public static int sdInDec;		// Super digits in a decimal holding an integer
		public static int sdInInt;		// Super digits in an Int32
		public static int sdInLong;		// Super digits in an Int64
		public static long baseML;		// radix baseM as type long
		public static long baseSq;		// base squared
		public static long bSMB;		// (baseSq - baseM) for MuDiv and MuSqRt
		public static long nDone;		// number if iterations of an innermost loop
		public static long nToGo;		// number of "Digits" left to compute
		public static MultiSI fMB;		// Multi-precision modulo base, 0 => normal arithmetic.
		public static MultiSI max64;	// an SI =  Int64.MaxValue
		public static MultiSI min64;	// an SI = -Int64.MaxValue
		public static MultiSI remSI;	// holds the remainder of infix operator SI / SI
		public static MultiSI sMOne;	// sMOne   = -1
		public static MultiSI sZero;	// sZero   =  0
		public static MultiSI sOne;		// sOne    =  1
		public static MultiSI sTwo;		// sTwo    =  2
		public static MultiSI sThree;	// sThree  =  3
		public static MultiSI sFour;	// sFour   =  4
		public static MultiSI sFive;	// sFive   =  5
		public static MultiSI s100;		// s100    =  100
		public static RichTextBox outS; // screen Output destination
		public static StreamWriter outL;// log file output #
		public static string nL;		// one Environment.NewLine \n or \r or \r\n
		public static string nL2;		// nL + nL
		public static string nLc;		// line continuation mark for WriteT nL or "-" + nL
	}
}
