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

namespace MyNamespace
{
	/// <summary>
	///		Multiple precision floating-point algorithms for computing Pi = 3.14159...
	///		Developed in Turbo Pascal 5.5.  Converted to Borland C++ 3.1, 4.0 and
	///		5.02 for Windows. Now in 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>

	public class MultiPi
	{
		public MultiPi()
		{
			fA = new MultiF();
			fB = new MultiF();
			fJ = new MultiF();
		}

		public static void QuRtIvF(MultiF x, MultiF z) // z = x ^ (-1/4), inverse 4th root
		// x is not changed, it is ok for x and this to have the same location in memory
		// but then x will be changed. Sets muErr True if error
		// Sets muErr True if x = 0
		{
			int sh, saveFMC;
			int mL = 0;
			sh = Math.Min(z.M, Calc.fMC) + 1;
			if (!MultiFD.expanded)  sh += MultiFD.fGuard;	// n guard "digits", no sticky bit
			if (sh > MultiID.muNMax)  sh = MultiID.muNMax;
			if (!Object.ReferenceEquals(x, z))	{ mL = z.M;  MultiFD.ReSizeF(1, z); }
			MultiF t = new MultiF(sh);
			if (sh > 14) {
				MultiF xL = new MultiF(1);
				MultiF y = new MultiF(sh);
				int ps, pl = 7;
				saveFMC = Calc.fMC;
				Calc.fMC = pl+2;
				QuRtIvF(x, t);			// t ~= 1/SqRt(SqRt(x))
				while (pl <= sh)
				{
					pl += pl;
					if (pl > sh)  pl = sh + 1;
					int p;
					for(p = sh + 1; p > pl; p=(p+1)/2)  ;
					pl = p;
					Calc.fMC = pl+2;
					MultiFD.SqF(t, y);
					MultiFD.SqF(y, y);
					if (x.N > pl+2)
					{
						ps = x.N - (pl+2);
						MultiFD.ReducePrecisionF(ps, x, xL);
						MultiFD.MulF(xL, y, y);
					}
					else
					{
						MultiFD.MulF(x, y, y);
					}
					MultiFD.SubF(MultiFD.fOne, y, y);
					MultiFD.MulF(y, t, y);
					MultiFD.RDiv1F(4, y);
					MultiFD.RAddF(y, t); // (1 - x * t^4) * t / 4 + t // computing x ^ (-1/4)
				}
				xL.V = null;
				Calc.fMC = saveFMC;
				MultiFD.ReSizeF(1, y);
			}
			else
			{
				MultiFD.SqRootF(x, t);
				MultiFD.SqRootF(t, t);
				MultiFD.DivF(MultiFD.fOne, t, t);
			}
			if (!Object.ReferenceEquals(x, z))  MultiFD.ReSizeF(mL, z);
			MultiFD.SetToF(t, z);
		} // QuRtIv

		public static void SqRtIvF(MultiF x, MultiF z) // z = 1 / SqRt(x)
		// x is not changed, it is ok for x and this to have the same location in memory
		// but then x will be changed.  Sets muErr True if error
		// Sets muErr True if x = 0
		{
			int sh, saveFMC;
			int mL = 0;
			sh = Math.Min(z.M, Calc.fMC) + 1;
			if (!MultiFD.expanded)  sh += MultiFD.fGuard; // Four guard "digits", no sticky bit
			if (sh > MultiID.muNMax)  sh = MultiID.muNMax;
			if (!Object.ReferenceEquals(x, z))  { mL = z.M;  MultiFD.ReSizeF(1, z); }
			MultiF t = new MultiF(sh);
			if (sh > 14)
			{
				MultiF xL = new MultiF(1);
				MultiF y = new MultiF(sh);
				int ps, pl = 7;
				saveFMC = Calc.fMC;
				Calc.fMC = pl+2;
				SqRtIvF(x, t); // t ~= 1/SqRt(x)
				while (pl <= sh)
				{
					pl += pl;
					if (pl > sh)  pl = sh + 1;
					int p;
					for(p = sh + 1; p > pl; p=(p+1)/2)  ;
					pl = p;
					Calc.fMC = pl+2;
					if (x.N > pl+2)
					{
						ps = x.N - (pl+2);
						MultiFD.ReducePrecisionF(ps, x, xL);
						MultiFD.MulF(xL, t, y);
					}
					else
					{
						MultiFD.MulF(x, t, y);
					}
					MultiFD.MulF(y, t, y);
					MultiFD.SubF(MultiFD.fOne, y, y);
					MultiFD.MulF(y, t, y);
					MultiFD.RDiv1F(2, y);
					MultiFD.RAddF(y, t); // (1 - x * t * t) * t / 2 + t // computing 1/SqRt(x)
				}
				xL.V = null;
				Calc.fMC = saveFMC;
				MultiFD.ReSizeF(1, y);		// release memory
			}
			else
			{
				MultiFD.SqRootF(x, t);
				MultiFD.DivF(MultiFD.fOne, t, t);
			}
			if (!Object.ReferenceEquals(x, z)) MultiFD.ReSizeF(mL, z);		// NMax(mL);
			MultiFD.SetToF(t, z);
		} // SqRtIv

		public static void PiAa(MultiF x) // Compute x = Pi by Borwein algorithm a
		// Uses 3 large "Registers" x, y, and z. l is small.
		{
			int i, j;
			int it; // Number of iterations(...
			string st;
			double timeAtStart;
			double deltaTime;
			timeAtStart = MultiID.TimerF();
			MultiFD.ReSizeF(Calc.fMC, x);		// 1: working alpha and final output
			MultiF z = new MultiF(Calc.fMC);	// 2: Temp
			MultiF y = new MultiF(Calc.fMC);	// 3: y0, y1, ...
			MultiF l = new MultiF(2);			// -(2**N)
			it = (int)(Math.Log(1024.0 * (Calc.fMC * MultiID.muDMax + 3) / 2788.0) / Math.Log(2.0)) + 2;
			Common.DoubleSpace();
			int digits = (Calc.fMC - 1) * MultiID.muDMax + 1;
			st = "  Starting to compute Pi by Borwein algorithm a, " + digits + " decimal digits";
			MultiID.WriteLn(st);
			if (MultiID.muAbort || MultiID.muErr) // if not enough memory
			{
				MultiFD.ReSizeF(1, x);			// clear X
				return;
			}
			MultiFD.ValueF10("0.5", out i, x);	// X0 = 1/2 = 0.5
			MultiFD.SetTo1F(2, y);
			SqRtIvF(y, y);						// Y0 = 1/SqRt(2) = SqRt(0.5)
			MultiFD.SetTo1F(-1, l);				// l  = -1
			for (j = 1; j <= it; j++)
			{
				MultiID.MuInterrupt();
				if (MultiID.muErr) break;
				if (MultiID.muAbort)
				{
					MultiID.Diag("Operator abort");
					break;
				}
				st = "  Start of iteration " + j + " of " + it;
				MultiID.Write(st);
				MultiID.TimeDel();
				// y = (1 - SqRoot(1 - y**2)) / (1 + SqRoot(1 - y**2))
				// y = (SqRtIv(1 - y**2) - 1) / (SqRtIv(1 - y**2) + 1)
				MultiFD.RSqF(y);
				MultiFD.ChSF(y);
				MultiFD.RAdd1F(1, y);
				SqRtIvF(y, y);
				MultiFD.Add1F(y, 1, z);			// z = SqRtIv(1 - y**2) + 1
				MultiFD.RAdd1F(-1, y);			// = SqRtIv(1 - y**2) - 1
				MultiFD.RDivF(z, y);
				// x = (1 + y)**2 * x - 2**N * y
				MultiFD.RAddF(l, l);			// l = -2**1, -2**2, ...
				MultiFD.Add1F(y, 1, z);			// z = (1 + y)
				MultiFD.RSqF(z);
				MultiFD.RMulF(z, x);			// x = (1 + y)**2 * x
				MultiFD.MulF(y, l, z);			// z = - 2**N * y
				MultiFD.RAddF(z, x);
			}
			if (!MultiID.muAbort && !MultiID.muErr)
			{
				MultiID.Write("  All done but inverting alpha");
				MultiID.TimeDel();
				MultiFD.SetTo1F(1, z);			// PI ~= 1 / x
				MultiFD.DivF(z, x, x);
				MultiID.WriteLn2("  All done computing Pi");
				deltaTime = MultiID.TimerF() - timeAtStart;
				st = string.Format("PiAa, DT = {0:0.00} sec.  ", deltaTime);
				MultiID.WriteTDLn(st);
			}
			else
			{
				MultiFD.ReSizeF(1, x);	// clear X
			}
		} // PiAa

		public static void PiAb(MultiF x) // Compute x = Pi by Borwein algorithm b
		// Uses 3 large "Registers" x, y, and z. l is small.
		{
			int i, j;
			int it; // Number of iterations(...
			string st;
			double timeAtStart;
			double deltaTime;
			timeAtStart = MultiID.TimerF();
			MultiFD.ReSizeF(Calc.fMC, x);		// 1: working alpha and final output
			MultiF z = new MultiF(Calc.fMC);	// 2: Temp
			MultiF y = new MultiF(Calc.fMC);	// 3: y0, y1, ...
			MultiF l = new MultiF(2);			// -(2**(2 * N + 1))
			it = (int)(Math.Log(1024.0 * (Calc.fMC * MultiID.muDMax + 3) / 2788.0) / Math.Log(4.0)) + 1;
			Common.DoubleSpace();
			int digits = (Calc.fMC - 1) * MultiID.muDMax + 1;
			st = "  Starting to compute Pi by Borwein algorithm b, " + digits + " decimal digits";
			MultiID.WriteLn(st);
			if (MultiID.muAbort || MultiID.muErr) // if not enough memory
			{
				MultiFD.ReSizeF(1, x);			// clear X
				return;
			}
			MultiFD.ValueF10("0.5", out i, y);	// Y0 = SqRt(2) - 1 = 1/SqRt(0.5) - 1
			SqRtIvF(y, y);
			MultiFD.SetToF(y, x);				// Temp x = SqRt(2)
			MultiFD.RAdd1F(-1, y);
			MultiFD.RMul1F(-4, x);				// X0 = 6 - 4 * SqRt(2);
			MultiFD.RAdd1F(6, x);
			MultiFD.SetTo1F(-2, l);
			for (j = 1; j <= it; j++)
			{
				MultiID.MuInterrupt();
				if (MultiID.muErr) break;
				if (MultiID.muAbort)
				{
					MultiID.Diag("Operator abort");
					break;
				}
				st = "  Start of iteration " + j + " of " + it;
				MultiID.Write(st);
				MultiID.TimeDel();
				// y = (1 - 4thRoot(1 - y**4)) / (1 + 4thRoot(1 - y**4))
				// y = (Inverse4thRoot(1 - y**4) - 1) / (inverse4thRoot(1 - y**4) + 1)
				MultiFD.RSqF(y);
				MultiFD.RSqF(y);
				MultiFD.ChSF(y);
				MultiFD.RAdd1F(1, y);
				QuRtIvF(y, y);
				MultiFD.Add1F(y, 1, z);			// z = Inverse4thRoot(1 - y**4) + 1
				MultiFD.RAdd1F(-1, y);			// = Inverse4thRoot(1 - y**4) - 1
				MultiFD.RDivF(z, y);
				// x = (1 + y)**4 * x - 2**(2*N + 1) * y * (1 + (1 + y)*y)
				MultiFD.RMul1F(4, l);			// l = -2**3, -2**5, ...
				MultiFD.Add1F(y, 1, z);			// z = (1 + y)
				MultiFD.RSqF(z);
				MultiFD.RSqF(z);
				MultiFD.RMulF(z, x);			// = (1 + y)**4 * x
				MultiFD.Add1F(y, 1, z);			// z = (1 + y)
				MultiFD.RMulF(y, z);
				MultiFD.RAdd1F(1, z);
				MultiFD.RMulF(y, z);
				MultiFD.RMulF(l, z);			// z = - 2**(2*N + 1) * y * (1 + (1 + y)*y
				MultiFD.RAddF(z, x);
			}
			if (!MultiID.muAbort && !MultiID.muErr)
			{
				MultiID.Write("  All done but inverting alpha");
				MultiID.TimeDel();
				MultiFD.SetTo1F(1, z);			// PI ~= 1 / x
				MultiFD.DivF(z, x, x);
				MultiID.WriteLn2("  All done computing Pi");
				deltaTime = MultiID.TimerF() - timeAtStart;
				st = string.Format("PiAb, DT = {0:0.00} sec.  ", deltaTime);
				MultiID.WriteTDLn(st);
			}
			else
			{
				MultiFD.ReSizeF(1, x);			// clear x
			}
		} // PiAb

		public static void PiGL(MultiF x, bool SpeedUp)
		// Compute x = Pi by Gauss-Legendre algorithm.
		// See Yasumasa KANADA Computer Centre, University of Tokyo.
		// Uses 4 large "Registers" x, y, z, and t. l is small.
		// if SpeedUp then speed up SqRootF() with guess, optional, default = True
		{
			int i, j;
			int it;		// Number of iterations(...
			float plG;	// Places the guess to SqRoot is good, in "Super" digits
			float plG0 = (float)((1396.0 / 1024.0) / MultiID.muDMax);
			string st;
			double timeAtStart;
			double deltaTime;
			timeAtStart = MultiID.TimerF();
			MultiFD.ReSizeF(Calc.fMC, x);		// 1: working number and final output
			MultiF t = new MultiF(Calc.fMC);	// 2: Temp
			MultiF z = new MultiF(Calc.fMC);	// 3: Temp
			MultiF y = new MultiF(Calc.fMC);	// 4: y0, y1, ...
			MultiF l = new MultiF(2);			// -(2**N)
			it = (int)(Math.Log(1024.0 * (Calc.fMC * MultiID.muDMax + 3) / 2788.0) / Math.Log(2.0)) + 1;
			plG = 0;
			if (SpeedUp)  plG = plG0;
			Common.DoubleSpace();
			int digits = (Calc.fMC - 1) * MultiID.muDMax + 1;
			st = "  Starting to compute Pi by Gauss-Legendre algorithm, " + digits + " decimal digits";
			int tot = 0;
			MultiID.WriteT(st, ref tot);
			if (SpeedUp)
			{
				MultiID.WriteTLn("", ref tot);
			}
			else
			{
				MultiID.WriteTLn(" without speed up", ref tot);
			}
			MultiFD.SetTo1F(1, x);				// X0 = 1
			MultiFD.SetTo1F(2, y);				// Y0 = 1/SqRt(2) = 0.707...
			SqRtIvF(y, y);
			MultiFD.ValueF10("0.25", out i, t);	// T0 = 1/4
			MultiFD.SetTo1F(-1, l);				// L0 = -1
			for (j = 1; j <= it; j++)
			{
				MultiID.MuInterrupt();
				if (MultiID.muErr) break;
				if (MultiID.muAbort)
				{
					MultiID.Diag("Operator abort");
					break;
				}
				st = "  Start of iteration " + j + " of " + it;
				MultiID.Write(st);
				MultiID.TimeDel();
				MultiFD.SetToF(x, z);		// z = x
				MultiFD.RAddF(y, x);		// x = (x+y) / 2
				MultiFD.RDiv1F(2, x);
				MultiFD.RMulF(z, y);		// y = SqRt(y * z)
				plG += plG;
				MultiFD.SqRootF(y, x, (int)plG, y);
				MultiFD.RSubF(x, z);		// z - x
				MultiFD.RSqF(z);			// (z - x) ** 2
				MultiFD.RMulF(l, z);		// l * (z - x) ** 2
				MultiFD.RAddF(z, t);		// t = t + l * (z - x) ** 2
				MultiFD.RAddF(l, l);		// l = 2 * l
			}
			if (!MultiID.muAbort && !MultiID.muErr)
			{
				MultiID.Write("  All done but ((x+y)^2) / (4*t)");
				MultiID.TimeDel();
				MultiFD.RAddF(y, x);		// PI ~= ((x+y)^2) / (4*t)
				MultiFD.RSqF(x);
				MultiFD.RDivF(t, x);
				MultiFD.ValueF10("0.25", out i, l);
				MultiFD.RMulF(l, x);
				MultiID.WriteLn2("  All done computing Pi");
				deltaTime = MultiID.TimerF() - timeAtStart;
				st = string.Format("PiGL, DT = {0:0.00} sec.  ", deltaTime);
				MultiID.WriteTDLn(st);
			}
			else
			{
				MultiFD.ReSizeF(1, x);		// clear x
			}
		} // PiGL

		public static void PiAGM(MultiF x, bool SpeedUp)
		// Compute x = Pi by Schoenhage's AGM algorithm.
		// See C++ author = Joerg Arndt		email: arndt@jjj.de
		// the C++ software is online at	http://www.jjj.de/
		// Uses 5 large "Registers" x, a, b, t, and s.
		// if SpeedUp then initialize with first iteration
		{
			int i, j, k;
			int it;			// Number of iterations(...
			string st;
			int pow2;	// 2^(k+1)
			double timeAtStart;
			double deltaTime;
			int n = Calc.fMC;
			timeAtStart = MultiID.TimerF();
			MultiFD.ReSizeF(Calc.fMC, x);		// 1: working number and final output
			MultiF a = new MultiF(Calc.fMC);	// 2: arithmetic mean
			MultiF b = new MultiF(Calc.fMC);	// 3: geometric mean
			MultiF t = new MultiF(Calc.fMC);	// 4: t
			MultiF s = new MultiF(Calc.fMC);	// 5: s
			it = (int)(Math.Log(1024.0 * (Calc.fMC * MultiID.muDMax + 3) / 2788.0) / Math.Log(2.0)) + 1;
			Common.DoubleSpace();
			int digits = (Calc.fMC - 1) * MultiID.muDMax + 1;
			st = "  Starting to compute Pi by Schoenhage's AGM algorithm, " + digits + " decimal digits";
			MultiID.WriteLn(st);
			if (MultiID.muAbort || MultiID.muErr)	// if not enough memory
			{
				MultiFD.ReSizeF(1, x);				// clear x
				return;
			}
			if (SpeedUp)
			{
				MultiFD.SetTo1F(2, b);				// b = 2
				SqRtIvF(b, b);						// b = 1/SqRt(2) = SqRt(1/2)	(SqRt inverse function)
				MultiFD.Div1F(b, 2, x);				// x = SqRt(1/2)/2
				MultiFD.SetToF(x, t);				// t = x
				MultiFD.ValueF10("0.5", out i, a);	// a = 1/2 = 0.5
				MultiFD.AddF(a, x, a);				// a = sqrt(1/2)/2 + 1/2
				MultiFD.ValueF10("0.125", out i, s);	// s = 1/8 = 0.125
				MultiFD.SubF(t, s, t);				// t = sqrt(1/2)/2 - 1/8
				MultiFD.Mul1F(s, 3, s);				// s = 3/8
				MultiFD.AddF(x, s, x);				// x = sqrt(1/2)/2 + 3/8
				MultiFD.Mul1F(t, 2, t);				// t = sqrt(1/2) - 1/4
				k = 1;
			}
			else
			{
				MultiFD.SetTo1F(1, a);				// a = 1
				MultiFD.SetTo1F(1, x);				// x = a^2 = 1
				MultiFD.ValueF10("0.5", out i, b);	// b = 1/2 = 0.5
				MultiFD.SetToF(b, t);				// t = 1/2 = 0.5
				k = 0;
				it++;
			}
			long lastPrec = 0;
			for (j = 1;  k < 30;  j++, k++)
			{
				MultiID.MuInterrupt();
				if (MultiID.muErr) break;
				if (MultiID.muAbort)
				{
					MultiID.Diag("Operator abort");
					break;
				}
				st = "  Start of iteration " + j + " of " + it;
				MultiID.Write(st);
				MultiID.TimeDel();
				MultiFD.AddF(x, b, s);			// s = x + b
				MultiFD.Div1F(s, 4, s);			// s = (x+b) / 4
				MultiFD.SqRootF(b, b);			// b = sqrt(b)		(Full square root)
				MultiFD.AddF(a, b, a);			// a' = a + b
				MultiFD.Div1F(a, 2, a);			// a' = (a+b) / 2
				MultiFD.SqF(a, x);				// x = a'^2			(Full square multiply)
				MultiFD.SubF(x, s, b);			// b = x - s
				MultiFD.Mul1F(b, 2, b);			// b = 2*(x-s)
				MultiFD.SubF(x, b, s);			// s = x - b
				pow2 = 1 << (k+1);				// pow2 = 2^(k+1)
				MultiFD.Mul1F(s, pow2, s);		// s = 2^(k+1) * (x-b)
				MultiFD.SubF(t, s, t);			// t = t - 2^(k+1) * (x-b)

				if (lastPrec >= n / 2) break;
				lastPrec = -s.N - s.C - 2;
			}
			if (!MultiID.muAbort && !MultiID.muErr)
			{
				MultiID.Write("  All done but calculation pi = 2*a*a/t");
				MultiID.TimeDel();
				MultiFD.SqF(a, a);				// a' = a*a
				MultiFD.Mul1F(a, 2, a);			// a' = 2*a*a
				MultiFD.DivF(a, t, x);			// pi = 2*a*a/t
				MultiID.WriteLn2("  All done computing Pi");
				deltaTime = MultiID.TimerF() - timeAtStart;
				st = string.Format("PiAGM, DT = {0:0.00} sec.  ", deltaTime);
				MultiID.WriteTDLn(st);
			}
			else
			{
				MultiFD.ReSizeF(1, x);			// clear x
			}
		} // PiAGM

		// Calculates pi using the binary splitting algorithm of the Chudnovsky brothers
		// Code is essentially taken from CLN by Bruno Haible
		// Edited by Harry J. Smith, 2006-03-27

		private static MultiF fA, fB, fJ;
//		private static long maxK, currK;
//		private static double oldPct;
//		private static double timeAtStart;

		private static MultiF Ba(int n)
		{
			MultiF v = new MultiF(Calc.fMC);
			v &= fA + fB * n;
			v.S = ((n & 1) == 1);
			return v;
		}

		private static MultiF Bp(int n)
		{
			if (n == 0) return MultiFD.fOne;
			MultiF f = new MultiF(2);
			f &= n;
			MultiF sixf = new MultiF(Calc.fMC);
			sixf &= MultiFD.fSix * f;
			return (sixf - MultiFD.fOne) * (MultiFD.fTwo * f - MultiFD.fOne) * (sixf - MultiFD.fFive);
		}

		private static MultiF Bq(int n)
		{
			if (n == 0) return MultiFD.fOne;
			MultiF f = new MultiF(2);
			f &= n;
			return fJ * f * f * f;
		}

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

		private static void PiR(int n1, int n2, MultiF fT, MultiF fQ, MultiF fP)
		{
			if (MultiID.muErr) return;
			if (MultiID.muAbort)
			{
				MultiID.Diag("Operator abort");
				return;
			}
			MultiF p0;
			MultiF p01;
			MultiF p012;
			switch (n2 - n1)
			{
				case 0:
					n2 = n1;
					MultiID.Write("PiCh case 0: n1 = " + n1);
					break;
				case 1:
					p0 = new MultiF(Calc.fMC);
					p0 &= Bp(n1);
					fT &= Ba(n1) * p0;
					fQ &= Bq(n1);
					if (fP != MultiFD.fZero) fP &= p0;
					MultiID.currK += 1;
					MultiID.Write("PiCh case 1: n1 = " + n1);
					break;
				case 2:
					p0 = new MultiF(Calc.fMC);
					p0 &= Bp(n1);
					p01 = new MultiF(Calc.fMC);
					p01 &= p0 * Bp(n1 + 1);
					MultiF q1 = new MultiF(Calc.fMC);
					q1 &= Bq(n1 + 1);
					fT &= q1 * Ba(n1) * p0 + Ba(n1 + 1) * p01;
					fQ &= Bq(n1) * q1;
					if (fP != MultiFD.fZero) fP &= p01;
					MultiID.currK += 4;
					break;
				case 3:
					p0 = new MultiF(Calc.fMC);
					p0 &= Bp(n1);
					p01 = new MultiF(Calc.fMC);
					p01 &= p0 * Bp(n1 + 1);
					p012 = new MultiF(Calc.fMC);
					p012 &= p01 * Bp(n1 + 2);
					MultiF q2 = new MultiF(Calc.fMC);
					q2 &= Bq(n1 + 2);
					MultiF q12 = new MultiF(Calc.fMC);
					q12 &= Bq(n1 + 1) * q2;
					fT &= q12 * Ba(n1) * p0 + q2 * Ba(n1 + 1) * p01 + Ba(n1 + 2) * p012;
					fQ &= Bq(n1) * q12;
					if (fP != MultiFD.fZero) fP &= p012;
					MultiID.currK += 8;
					break;
				case 4:
					p0 = new MultiF(Calc.fMC);
					p0 &= Bp(n1);
					p01 = new MultiF(Calc.fMC);
					p01 &= p0 * Bp(n1 + 1);
					p012 = new MultiF(Calc.fMC);
					p012 &= p01 * Bp(n1 + 2);
					MultiF p0123 = new MultiF(Calc.fMC);
					p0123 &= p012 * Bp(n1 + 3);
					MultiF q3 = new MultiF(Calc.fMC);
					q3 &= Bq(n1 + 3);
					MultiF q23 = new MultiF(Calc.fMC);
					q23 &= Bq(n1 + 2) * q3;
					MultiF q123 = new MultiF(Calc.fMC);
					q123 &= Bq(n1 + 1) * q23;
					fT &= q123 * Ba(n1) * p0 +
						q23 * Ba(n1 + 1) * p01 +
						q3 * Ba(n1 + 2) * p012 +
						Ba(n1 + 3) * p0123;
					fQ &= Bq(n1) * q123;
					if (fP != MultiFD.fZero) fP &= p0123;
					MultiID.currK += 12;
					break;
				default:
					int nM = (n1 + n2) / 2;
					MultiF fLT = new MultiF(Calc.fMC);
					MultiF fLQ = new MultiF(Calc.fMC);
					MultiF fRT = new MultiF(Calc.fMC);
					MultiF fRQ = new MultiF(Calc.fMC);
					MultiF fLP = new MultiF(Calc.fMC);
					MultiF fRP = new MultiF(Calc.fMC);
					fLP &= MultiFD.fOne;
					fRP &= fP;
					if (fP == MultiFD.fZero) fRP &= MultiFD.fZero;
					else fRP &= MultiFD.fOne;
					PiR(n1, nM, fLT, fLQ, fLP);
					PiR(nM, n2, fRT, fRQ, fRP);
					fT &= fRQ * fLT + fLP * fRT;
					fQ &= fLQ * fRQ;
					if (fP != MultiFD.fZero) fP &= fLP * fRP;
					MultiID.currK += n2 - n1;
					MultiID.MuInterrupt();	// allow for abort
					break;
			}
			ReportPiCh();
		} // PiR

		public static void PiCh(MultiF x) // Compute x = Pi by Chudnovsky brothers' binary splitting algorithm
		{
			MultiID.timeAtStart = MultiID.TimerF();
			int n;
			int Basedigit = 10;
			MultiID.WriteLn("Using the Chudnovsky brothers' binary splitting algorithm to compute Pi");
			MultiFD.ReSizeF(Calc.fMC, fA);
			MultiFD.ReSizeF(Calc.fMC, fB);
			MultiFD.ReSizeF(Calc.fMC, fJ);
			MultiF fT = new MultiF(Calc.fMC);
			MultiF fQ = new MultiF(Calc.fMC);
			fA &= 13591409;
			fB &= 545140134;
			fJ &= 107701824;
			fJ *= 101568;
			fJ *= 1000;		// fJ = 10939058860032000
			int destprec = Calc.fMC * MultiID.muDMax;
			n = (int)(destprec * Math.Log((double)Basedigit) / 32.65445004177);
			MultiID.maxK = (int)(n * (Math.Log((double)n) / Math.Log(2.0) + 1));
			MultiID.currK = 0;
			MultiID.oldPct = 0.0;
			MultiF fZeroL = new MultiF(); 
			PiR(0, n + 1, fT, fQ, fZeroL);
			MultiID.currK = MultiID.maxK;
			if(MultiID.oldPct < 100.0) ReportPiCh();
			fA &= 640320;
			SqRtIvF(fA, fB);
			fA &= fB * fT;
			MultiFD.ReSizeF(Calc.fMC, x);
			x &= 53360 * fQ / fA;
			MultiFD.ReSizeF(0, fA);
			MultiFD.ReSizeF(0, fB);
			MultiFD.ReSizeF(0, fJ);
		} // PiCh
	}
}
