using System;
using System.Windows.Forms; // for Application


namespace MyNamespace
{
	/// <summary>
	///		PrimeFA - Multi-precision prime factors algorithm
	///		Factoring into primes. See Knuth, Vol. 2, Page 348
	///		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>

	public class PrimeFA
	{
		public PrimeFA()				// constructor logic
		{
			lenPTab = 33000;			// 200 low, 3512 at max for Int16 primes, was 32766
			rS = 8;						// number of sieve tables (r) for SieveMethod
			maxPs = 2000; 				// constant, Max number of prime factors, 24 for 36 digit numbers
										// 54 for 101 digit numbers, about 1400 for 5000 digit #s
			cutoff = new MultiSI();		// limit for using MethodX
			delX = new MultiSI();		// for SieveMethod
			remS = new MultiSI();		// for SieveMethod
			root2 = new MultiSI();		// for SieveMethod
			sieveN = new MultiSI();		// for SieveMethod
			start = new MultiSI();		// minimum factor to start prime test, normally 27
			z = new MultiSI();			// factor of zIn not yet factored
			zIn = new MultiSI();		// original input number to be factored
			log2OfBase = Math.Log(MultiID.baseM , 2);
			sdBits = Len(MultiID.maxDigit) - 1;		// bits in maxDigit - 1
			twoToSDB = (int)Math.Pow(2, sdBits);	// 2 ^ sdBits = largest power of 2 < maxDigit
			complete = false;
			exitEarly = false;
			isSquare = true;			// assume true until proven otherwise
			squareFree = true;
			testSqFree = false;
			testSquare = false;
			testPrimR = false;
			byECM = false;
			stopDiv = 131072;			// 2^17 for PrimeDiv function, value can change
			PhiLSIDepth = 0;
			PiLDepth = 0;
			pTabFilled = false;
			pTab = new int[0];			// pTab[1 To lenPTab], pTab[1] = 2
		}

		public static void FillPTab()	// fill the prime table pTab[1 TO lenPTab]
		{
			MultiID.ADiag("Wait, filling the " + lenPTab + " entry prime table");
			int i;						// index to old entry in prime table
			int mC = lenPTab;			// max value allowed for mC (mC = mCMax)
			pTab = new int[lenPTab+1];	// pTab[1 To lenPTab], pTab[1] = 2
			pTab[0] = 1;				// [0] = 1 to be compatible with UBASIC prm(n) function
			new NPrimes(mC, 115);		// init prime number generator, 115 is the nr returned last time
			int nr = NPrimes.NPrime(mC);// nr = the current length of the internal tables returned
			i = nr;						// needed so nr can be monitored
			for (i = 1; i <= lenPTab; i++)
				pTab[i] = (int)NPrimes.iP[i];
			NPrimes.NPrimeDone();
			squareMaxP = pTab[lenPTab];
			squareMaxP *= squareMaxP;
			pTabFilled = true;
			if (MultiID.diagOn)
			{
				for (i = lenPTab - 1; i <= lenPTab; i++)
					MultiID.WriteLn(i + " " + pTab[i]);
				MultiID.WriteLn2("Square Max P = " + squareMaxP);
			}
		}

		private static void InitP()		// init NPrime prime number generator
		{
			int mC = lenPTab;			// max value allowed for mC (mC = mCMax)
			new NPrimes(mC, 200000);	// init prime number generator
			int nr = NPrimes.NPrime(mC);// nr = the current length of the internal tables returned
			xI = 1;
		}

		private static void NextP(ref Int64 p) // p = first prime > p
		{
			int nr;						// the current length of the internal tables returned
			int mC = lenPTab;			// max value allowed for mC (mC = mCMax)
			if (NPrimes.jLimM != mC)
				InitP();
			if (p < 2)
			{
				p = 2;
				xI = 1;
				return;
			}
			if (p == NPrimes.iP[xI])	// normal case
			{
				xI++;
				if (xI > mC)
				{
					nr = NPrimes.NPrime(mC);
					xI = 1;
				}
				p = (Int64)NPrimes.iP[xI];
				return;
			}
			if (p < NPrimes.iP[xI])
			{
				p = (Int64)NPrimes.iP[xI];
				return;
			}
			if (p < NPrimes.iP[1])
			{
				new NPrimes(mC, 200000);		// init prime number generator
				nr = NPrimes.NPrime(mC);
			}
			while (p > NPrimes.iP[mC])
				nr = NPrimes.NPrime(mC);
			xI = 1;
			while (p > NPrimes.iP[xI])
				xI++;
			if (p == NPrimes.iP[xI]) xI++;
			if (xI > mC)
			{
				nr = NPrimes.NPrime(mC);
				xI = 1;
			}
			p = (Int64)NPrimes.iP[xI];
		} // NextP

		private static void Init1() // initialize program basics
		{
			nF = 0;								// no factors yet
			primes.a = new MultiSI[maxPs+1];	// prime factor table
			powers.a = new int[maxPs+1];		// factor power table
			MultiID.ReSizeSI(size, zIn);		// init object
		}

		private static void Init() // initialize program
		{
			Init1();
			int mRM = 22;						// mRM= pTab[rS + 1] - 1, (m sub r - 1) for SieveMethod
			sieve = new bool[rS+1, mRM+1];		// for SieveMethod // rS tables
			squareMod = new int[mRM+1];
			kS = new int[rS+1];
			d = new int[12];
			d[0] = 0;
			d[1] = 1; d[2] = 2; d[3] = 2;
			d[4] = 4; d[5] = 2; d[6] = 4;
			d[7] = 2; d[8] = 4; d[9] = 6;
			d[10] = 2; d[11] = 6;
			noSieve = false;
			noMoreDiv = false;
			MultiID.ReSizeSI(sizeO2, start);	// init objects
			MultiID.ReSizeSI(size, sieveN);
			MultiID.ReSizeSI(size, remS);
			MultiID.ReSizeSI(sizeO2, root2);
			MultiID.ReSizeSI(sizeO2, delX);
			int sD = 1 + (18 / MultiID.muDMax);
			MultiID.ReSizeSI(sD, cutoff);		// sD is the number of super digits in 19 digit cutoff
			int i;
			MultiID.ValueSI("9223,37203,68547,75807", out i, cutoff); // 2^63-1
		}

		private static void Fini1() // clean up your act
		{
			MultiID.DoneSI(start);
			MultiID.DoneSI(sieveN);
			MultiID.DoneSI(remS);
			MultiID.DoneSI(root2);
			MultiID.DoneSI(delX);
			MultiID.DoneSI(cutoff);
		}

		private static void Fini2() // clean up your act
		{
			MultiID.DoneSI(zIn);
			MultiID.DoneSI(z);
		}

		// old DOS program Help notes:
		// PrimeFA - Multi-precision prime factor algorithm
		// Version 2.50, last revised: 1994-11-20, 1100 hours
		// Copyright (c) 1981-1994 by author: Harry j Smith,
		// 19628 Via Monte Dr., Saratoga CA 95070.  All rights reserved.
		//
		// Prime factors of Integers <= 5001 decimal digits (using 80387 coprocessor)
		// Some may take forever (over 300,000,000 years)
		//
		// Enter q to Quit, E to Echo screen to file: PrimeFA.Out,
		// d for Diagnostics, S to force sieve method, t to toggle find direction,
		// M to force More divide method, n for Normal
		//
		// Can also enter numbers like 5*(2+3)/2-3^(2!)
		// (= => +,   . => *,   , => ignored,   P => ^,   F => !,
		// ? => Help,   -number => find a prime)
		//
		// Or enter numbers like 1234s12
		// to start testing 1234 at start divisor 12
		//
		// Or enter numbers like 9r8,
		// to repeat 9, 8 times i.e. 99999999
		// (use ESC to interrupt if something runs too long)
		//
		// Number to factor: (type q <Enter> to Quit, MemAvail = 442224)
		// end Help

		private static void DelPrimes(ref int dNF, TPrimes dPrimes) // delete prime factors found
		{
			for (int i = 1; i <= dNF; i++)
				dPrimes.a[i] = null;
			dNF = 0;
		}

		private static void ReportFactor(			// part of FactorN, reports div as a factor
						bool report,				// report		= true; to report factors to output screen
						ref int fNF,				// fNF			= Number of prime factors found
						ref TPrimes fPrimes,		// fPrimes		= Array of prime factors found
						Tpowers fPowers,			// fPowers		= Array of powers of the prime factors found
						out bool exitFactor,		// exitFactor	= Exit Factor flag
						MultiSI div,				// div			= Divisor
						int po)						// po			= Power
		{
			int tot = 0;
			exitFactor = false;
			if (MultiID.muAbort) return;
			MultiID.ADebug("Enter ReportFactor");
			if (MultiID.statReq2 || MultiID.debug)
			{
				MultiID.WriteT("Report=", ref tot);
				if (report)
					MultiID.WriteT("True", ref tot);
				else
					MultiID.WriteT("False", ref tot);
				MultiID.WriteT(" NF=" + fNF + 1 + " Po=" + po + " Factor=", ref tot);
				MultiID.WriteSI(div, ref tot, 0, false);
				MultiID.TimeDel(ref tot);
				if (MultiID.statReq2) RunForm.FocusBothUs();
				MultiID.statReq2 = false;
			}
			if (fNF < maxPs)
			{
				fNF++;
				fPrimes.a[fNF] = new MultiSI(div.N);
				MultiID.SetToSI(div, fPrimes.a[fNF]);
				fPowers.a[fNF] = po;
			}
			else
			{
				MultiID.WriteSLn("Error, cannot store more than " + maxPs + " prime factors");
				exitEarly = true;
			}
			if (report && !exitEarly)
			{
				if (MultiID.diagOn)
				{
					MultiID.WriteSI(div, ref tot);
					MultiID.WriteLn(" ^ " + po);
				}
			}
			else
			{
				MultiID.ADebug("Exit FactorN from ReportFactor");
				exitFactor = true;
				return;
			}
			MultiID.ADebug("Exit ReportFactor");
		} // ReportFactor

		private static void SmallFactors(		// part of FactorN, find small factors in pTab
						MultiSI n,
						bool report,
						ref int fNF,
						TPrimes fPrimes,
						Tpowers fPowers,
						out bool exitFactor,
						MultiSI div,
						ref int po)
		{
			int i;								// index
			int sign;							// sign result of ABCompUI Sub
			MultiSI q = new MultiSI(size);		// quotient
			MultiSI r = new MultiSI(size);		// remainder of divide
			bool rZero;							// true if Remainder is zero
			Application.DoEvents();
			MultiID.MuInterrupt();
			MultiID.ADebug("Enter SmallFactors");
			exitFactor = false;
			i = 1;
			do
			{
				MultiID.SetTo1SI(pTab[i], div);
				i++;
				po = 0;
				do
				{
					MultiID.DivQRSI(n, div, r, q); // q = n / div
					rZero = MultiID.ZTestSI(r);
					if (rZero)
					{
						MultiID.SetToSI(q, n);
						po++;
					}
				}
				while (rZero && !MultiID.muAbort);
				if (po > 0)
				{
					ReportFactor(report, ref fNF, ref fPrimes, fPowers, out exitFactor, div, po);
					if ((report && findAPrime) || exitFactor)
					{
						MultiID.ADebug("Quick exit SmallFactors");
						MultiID.DoneSI(r);
						MultiID.DoneSI(q);
						return;
					}
				}
				sign = MultiID.AbCompUI(q, div);
			}
			while (sign < 0 && i <= lenPTab && !MultiID.muAbort && !exitEarly);
			MultiID.DoneSI(r);
			MultiID.DoneSI(q);
			if (!MultiID.EQ1UI(1, n) && !MultiID.GTDbSI(squareMaxP, n))
			{
				MultiID.SetToSI(n, div);
				po = 1;
				MultiID.SetTo1SI(1, n);
				ReportFactor(report, ref fNF, ref fPrimes, fPowers, out exitFactor, div, po);
				if (exitFactor)
				{
					MultiID.ADebug("ReportFactor exit SmallFactors");
					return;
				}
			}
			MultiID.ADebug("Exit SmallFactors");
		} // SmallFactors

		private static void MoreDiv(			// part of FactorN, try the division method some more
					MultiSI n,
					MultiSI start,
					bool report,
					ref int fNF,
					TPrimes fPrimes,
					Tpowers fPowers,
					out bool exitFactor,
					MultiSI div,
					ref int po,
					ref bool found,
					ref bool twoFactors)
		{
			int i;								// index
			int j;								// index
			int sign;							// sign result of ABCompUI Sub
			bool done;							// done flag
			bool rZero;							// true if Remainder is zero
			int tot = 0;
			Application.DoEvents();
			MultiID.MuInterrupt();
			exitFactor = false;
			if (noMoreDiv) return;
			MultiID.ADebug("Enter MoreDiv, Try the division method some more");
			found = false;
			sign = MultiID.AbCompUI(n, start);
			if (sign < 0)
			{
				MultiSI q = new MultiSI(size);		// quotient
				MultiSI r = new MultiSI(size);		// remainder of divide
				MultiSI MaxDi = new MultiSI(sizeO2);// max divisor this try of MoreDiv
				MultiID.Add1SI(start, 10000000, MaxDi);
				done = false;
				j = 0;
				po = 0;
				MultiID.Add1SI(start, -7, div);	// div = ((start - 7) / 1020) * 1020 + 7
				div /= 1020;
				div *= 1020;
				MultiID.Add1SI(div, 7, div);
				i = 4;
				do
				{
					MultiID.Add1SI(div, d[i], div); // div += d[i]
					if (i < 11)
						i++;
					else
					{
						i = 4;
						j++;
						if (j == 34)
						{
							j = 0;
							if (twoFactors)
							{
								sign = MultiID.AbCompUI(MaxDi, div);
								if (sign > 0) done = true;
							}
							else
							{
								MultiID.DivSI(q, div, q);
								sign = MultiID.AbCompUI(q, div); // stop at Cube root
								if (sign > 0)
								{
									done = true;
								}
							}
							if (MultiID.statReq2)
							{
								MultiID.ShortWrSI(40, n, ref tot);
								MultiID.WriteT(" / ", ref tot);
								MultiID.WriteSI(div, ref tot);
								MultiID.WriteT(" MoreDiv", ref tot);
								MultiID.TimeDel(ref tot);
								RunForm.FocusBothUs();
								MultiID.statReq2 = false;
							}
							if (MultiID.muAbort)
							{
								done = true;
								exitFactor = true;
								MultiID.ADebug("MuAbort exit MoreDiv");
								MultiID.DoneSI(MaxDi);
								MultiID.DoneSI(r);
								MultiID.DoneSI(q);
								return;
							}
						}
					}
					if (n <= cutoff)
					{
						MethodX(n, div, ref po, q, ref done);
					}
					else
					{
						do
						{
							MultiID.DivQRSI(n, div, r, q);
							rZero = MultiID.ZTestSI(r);
							if (rZero)
							{
								MultiID.SetToSI(q, n);
								po++;
								done = true;
							}
						}
						while (rZero && !MultiID.ZTestSI(n) && !MultiID.muAbort);
					}
				}
				while (!done && !MultiID.muAbort && !MultiID.muErr);
				MultiID.DoneSI(MaxDi);
				MultiID.DoneSI(r);
				found = (po > 0);
				if (found)
				{
					ReportFactor(report, ref fNF, ref fPrimes, fPowers, out exitFactor, div, po);
					if (exitFactor)
					{
						MultiID.ADebug("ReportFactor exit MoreDiv");
						MultiID.DoneSI(q);
						return;
					}
				}
				sign = MultiID.AbCompUI(q, div);
				if (!MultiID.EQ1UI(1, n) && (sign >= 0))
				{
					MultiID.SetToSI(n, div);
					po = 1;
					MultiID.SetTo1SI(1, n);
					ReportFactor(report, ref fNF, ref fPrimes, fPowers, out exitFactor, div, po);
					if (exitFactor)
					{
						MultiID.ADebug("ReportFactor exit MoreDiv");
						MultiID.DoneSI(q);
						return;
					}
				}
				MultiID.Add1SI(div, 2, start);
				MultiID.DivSI(q, div, q);
				sign = MultiID.AbCompUI(q, div);
				twoFactors = (sign > 0);			// (div > q / div), div > cube root of n
				MultiID.DoneSI(q);
			}
			MultiID.ADebug("Normal exit MoreDiv");
		} // MoreDiv

		private static void MethodX(		// extended Method, part of MoreDiv which is part of FactorN
						MultiSI n,
						MultiSI div,
						ref int po,
						MultiSI q,			// quotient
						ref bool done)
		{
			Int64 dX;
			Int64 qX;
			Int64 nX;
			Int64 rX;
			bool rZero;		// true if Remainder is zero
			int jj;
			int tot = 0;
			qX = 0;
			nX = 0;
			rX = 0;
			jj = 0;
			MultiID.GetExSI(out dX, div);
			MultiID.GetExSI(out nX, n);
			do
			{
				do
				{
					qX = Math.DivRem(nX, dX, out rX);
					rZero = (rX == 0);
					if (rZero)
					{
						nX = qX;
						po++;
						done = true;
					}
					jj++;
					if (jj == 1000)
					{
						jj = 0;
						MultiID.MuInterrupt();
						if (MultiID.muAbort) break; // exit do
						if (MultiID.statReq2)
						{
							tot = 0;
							MultiID.SetToExSI(nX, n);
							MultiID.SetToExSI(dX, div);
							MultiID.ShortWrSI(40, n, ref tot);
							MultiID.WriteT(" / ", ref tot);
							MultiID.WriteSI(div, ref tot);
							MultiID.WriteT(" MethodX", ref tot);
							MultiID.TimeDel(ref tot);
							RunForm.FocusBothUs();
							MultiID.statReq2 = false;
						}
					}
				}
				while (rZero && !MultiID.muAbort);
				if (!done)
				{
					NextP(ref dX);
				}
			}
			while (dX < qX && !MultiID.muAbort && !done);
			MultiID.SetToExSI(qX, q);
			MultiID.SetToExSI(nX, n);
			MultiID.SetToExSI(dX, div);
		} // MethodX

		private static void LargeFactor(	// part of FactorN, determine if large factor is prime
						MultiSI n,
						bool report,
						out bool exitFactor,
						MultiSI div,
						ref int po,
						ref bool found)
		{
			int j;					// index in pTab
			int x;					// base of PowMBSI calculation for prime test
			int maxJ;				// the number of primes that have been tested with PowMBSI
			int lNF;				// number of prime factors found in nP
			TPrimes lPrimes;		// prime factors found in nP
			Tpowers lPowers;		// powers of the prime factors found in nP
			lPrimes.a = new MultiSI[maxPs+1];	// prime factor table
			lPowers.a = new int[maxPs+1];		// factor power table
			int tot = 0;
			Application.DoEvents();
			MultiID.MuInterrupt();
			MultiID.ADebug("Enter LargeFactor");
			exitFactor = false;
			found = false;
			if (MultiID.muAbort)
			{
				exitFactor = true;
				return;
			}
			lNF = 0;
			//lPrimes new
			MultiSI m = new MultiSI(2 * size);		// result of i^(n-1) % n calculation (n-1 power)
			MultiSI nm = new MultiSI(size);			// n - 1
			MultiID.Add1SI(n, -1, nm);				// nm = n - 1
			MultiID.PowMBSI(3, nm, n, m);			// m = 3^nm % n
			if (!MultiID.GT1UI(1, m))
			{
				MultiSI nP = new MultiSI(size);			// n' , a secondary number to be factored
				MultiSI q = new MultiSI(size);			// quotient
				MultiSI pM = new MultiSI(2 * size);		// result of PowMBSI calculation for secondary number
				MultiSI start = new MultiSI(sizeO2);	// minimum factor to start prime test, = 2
				MultiID.PowMBSI(2, nm, n, m);
				if (!MultiID.GT1UI(1, m))
				{
					MultiID.PowMBSI(5, nm, n, m);
					if (!MultiID.GT1UI(1, m))
					{
						MultiID.PowMBSI(7, nm, n, m);
						maxJ = 4;
						MultiID.SetToSI(nm, nP);
						MultiID.SetTo1SI(2, start);
						while (!MultiID.EQ1UI(1, nP) && !MultiID.GT1UI(1, m))	// (nP != 1) and (m < 2)
						{
							DelPrimes(ref lNF, lPrimes);
							FactorN(nP, start, false, out lNF, lPrimes, lPowers);
							if (MultiID.muAbort)
							{
								exitFactor = true;
								DelPrimes(ref lNF, lPrimes);
								MultiID.DoneSI(start);
								MultiID.DoneSI(pM);
								MultiID.DoneSI(q);
								MultiID.DoneSI(nP);
								MultiID.DoneSI(nm);
								MultiID.DoneSI(m);
								return;
							}
							j = 1;
							while (j <= lenPTab)
							{
								x = pTab[j];
								if (j > maxJ)
								{
									MultiID.PowMBSI(x, nm, n, m);
									maxJ = j;
								}
								if (!MultiID.GT1UI(1, m)) // m < 2
								{
									MultiID.DivSI(nm, lPrimes.a[1], q);
									MultiID.PowMBSI(x, q, n, pM);
									if (!MultiID.EQ1UI(1, pM))	// pM != 1
										j = lenPTab;			// abort loop
									else if (j == lenPTab)
									{
										MultiID.SetTo1SI(-1, m);
									}
								}
								else
								{
									j = lenPTab;	// abort loop
								}
								j++;
							}
						}
						if (!MultiID.EQ1UI(1, pM) && !MultiID.GT1UI(1, m))
						{
							// (pM != 1) and (m < 2)
							if (MultiID.debug)
							{
								MultiID.WriteT("Large factor ", ref tot);
								MultiID.WriteSI(n, ref tot);
								MultiID.WriteTLn(" is prime", ref tot);
								MultiID.WriteLn("Should exit LargeFactor");
							}
							MultiID.SetToSI(n, div);
							po = 1;
							MultiID.SetTo1SI(1, n);
							found = true;
							if (report) findAPrime = false;
						}
					}
				}
				MultiID.DoneSI(start);
				MultiID.DoneSI(pM);
				MultiID.DoneSI(q);
				MultiID.DoneSI(nP);
			}
			MultiID.DoneSI(nm);
			if (MultiID.GT1UI(1, m)) // m > 1
			{
				if (MultiID.statReq2 || MultiID.debug)
				{
					MultiID.WriteT("Large factor is not prime", ref tot);
					MultiID.TimeDel(ref tot);
					if (MultiID.statReq2) RunForm.FocusBothUs();
					MultiID.statReq2 = false;
				}
				if (report && findAPrime)
				{
					exitFactor = true;
					DelPrimes(ref lNF, lPrimes);
					MultiID.ADebug("Exit FactorN from LargeFactor");
					MultiID.DoneSI(m);
					return;
				}
			}
			MultiID.DoneSI(m);
			DelPrimes(ref lNF, lPrimes);
			MultiID.ADebug("Exit LargeFactor");
		} // LargeFactor

		private static void SieveMethod(			// part of FactorN
						MultiSI n,
						MultiSI start,
						bool report,
						ref int fNF,
						TPrimes fPrimes,
						Tpowers fPowers,
						out bool exitFactor,
						MultiSI div,
						ref int po,
						ref bool found,
						bool twoFactors)
		{
			int LMax = 2000;		// number of iterations of sieve method per call
			int sign;				// sign result of ABCompUI Sub
			int kJ;					// j'th entry of sieve k table
			int i;					// index
			int j;					// index
			int iL;					// count of the number of iterations of sieve method
			bool rZero;				// true if Remainder of divide is zero
			bool allTrue;			// all match made in sieve table are true
			int sNF;				// number of prime factors found
			double t;				// value of Clock
			TPrimes sPrimes;		// prime factors found
			Tpowers sPowers;		// powers of the prime factors found
			int tot = 0;
			sPrimes.a = new MultiSI[maxPs+1];	// prime factor table
			sPowers.a = new int[maxPs+1];		// factor power table
			Application.DoEvents();
			MultiID.MuInterrupt();
			exitFactor = false;
			if (noSieve) return;
			MultiID.ADebug("Enter SieveMethod");
			sNF = 0;
			MultiSI q = new MultiSI(size);		// quotient
			MultiSI r = new MultiSI(size);		// remainder of divide
			MultiSI x = new MultiSI(size);		// x for search, factor = x - y
			MultiSI y = new MultiSI(size);		// if y**2 = x**2 - n, then x - y is the desired factor
			MultiSI dP = new MultiSI(size);		// double precision temp
			MultiSI root = new MultiSI(size);	// square root temp
			MultiSI remR = new MultiSI(size);	// remainder of square root
			MultiID.WriteLn();
			MultiID.WriteSI(zIn, ref tot);
			MultiID.WriteTLn(" is a hard nut to crack", ref tot);
			found = false;
			sign = MultiID.AbCompUI(sieveN, n);
			if (sign != 0) InitSieve(n);
			iL = 0;
			do
			{
				allTrue = true;
				for (i = 1; i <= rS; i++)
				{
					if (sieve[i, kS[i]])
						{;} // noOp
					else
					{
						allTrue = false;
						i = rS;
					}
				}
				if (allTrue)
				{
					if (MultiID.muAbort)
					{
						exitFactor = true;
						goto ExitSieve;
					}
					iL++;
					// root of (x * x - n) needed but x * x may overflow
					if (MultiID.statReq2)
					{
						MultiID.WriteT("DelX=X-Root(N)=", ref tot);
						MultiID.WriteSI(delX, ref tot, 0, false);
						MultiID.TimeDel(ref tot);
						MultiID.statReq2 = false;
					}
					MultiID.AddSI(delX, root2, remR); // sqRoot(delX * (delX + root2) - remS)
					MultiID.MulSI(remR, delX, remR);
					MultiID.SubSI(remR, remS, remR);
					MultiID.SqRtRemSI(remR, remR, root);
					if (MultiID.statReq2)
					{
						MultiID.SetToSI(root, y);
						MultiID.Div1SI(root2, 2, dP);		// x = delX + root2 / 2
						MultiID.AddSI(dP, delX, x);
						MultiID.SubSI(x, y, dP);			// dP = x - y
						MultiID.ShortWrSI(40, n, ref tot);
						MultiID.WriteT(" / ", ref tot);
						MultiID.ShortWrSI(40, dP, ref tot);
						MultiID.WriteTLn(" SieveMethod", ref tot);
						MultiID.WriteT("Y=Root(X*X - N)=", ref tot);
						MultiID.WriteSI(root, ref tot);
						MultiID.WriteT(" Rem=", ref tot);
						MultiID.WriteSI(remR, ref tot, 0, false);
						MultiID.TimeDel(ref tot);
						RunForm.FocusBothUs();
						MultiID.statReq2 = false;
					}
					if (MultiID.ZTestSI(remR))
					{
						MultiID.SetToSI(root, y);
						MultiID.Div1SI(root2, 2, dP);		// x = delX + root2 / 2
						MultiID.AddSI(dP, delX, x);
						MultiID.SubSI(x, y, dP);			// dP = x - y
						MultiID.WriteSI(dP, ref tot);
						MultiID.Write(" is a big factor of ");
						MultiID.WriteSILn(n, ref tot, 0, false);
						if (twoFactors)
						{
							MultiID.SetToSI(dP, div);
							po = 1;
							MultiID.DivSI(n, div, dP);		// dP = n / div
							sign = MultiID.AbCompUI(dP, div);
							if (sign == 0)
								po++;
							else
							{
								MultiID.SetToSI(dP, n);		// this was missing from Apple )( program
								ReportFactor(report, ref fNF, ref fPrimes, fPowers, out exitFactor, div, po);
								if (exitFactor || MultiID.muAbort)
								{
									exitFactor = true;
									goto ExitSieve;
								}
								MultiID.SetToSI(dP, div);
							}
							MultiID.SetTo1SI(1, n);
							ReportFactor(report, ref fNF, ref fPrimes, fPowers, out exitFactor, div, po);
							if (exitFactor || MultiID.muAbort)
							{
								exitFactor = true;
								goto ExitSieve;
							}
						}
						else
						{
							FactorN(dP, start, false, out sNF, sPrimes, sPowers);
							MultiID.SetToSI(sPrimes.a[1], div);
							po = sPowers.a[1];
							for (j = 1; j <= po; j++)
								MultiID.DivSI(n, div, n);		// n = n / div
							do
							{
								MultiID.DivQRSI(n, div, r, q);	// q = n / div
								rZero = MultiID.ZTestSI(r);
								if (rZero)
								{
									MultiID.SetToSI(q, n);
									po++;
								}
							}
							while (rZero && !MultiID.muAbort);
							ReportFactor(report, ref fNF, ref fPrimes, fPowers, out exitFactor, div, po);
							if (exitFactor || MultiID.muAbort)
							{
								exitFactor = true;
								goto ExitSieve;
							}
						}
						found = true;
						iL = LMax;
					}
				}
				// step x + kS
				MultiID.RAdd1SI(1, delX);
				for (j = 1; j <= rS; j++)
				{
					kJ = kS[j] - 1;
					if (kJ < 0)
						kS[j] = pTab[j + 1] - 1;
					else
					{
						kS[j] = kJ;
					}
				}
			}
			while (iL < LMax && !MultiID.muErr && !MultiID.muAbort);
		ExitSieve:
			MultiID.DoneSI(remR);
			MultiID.DoneSI(root);
			MultiID.DoneSI(dP);
			MultiID.DoneSI(y);
			MultiID.DoneSI(x);
			MultiID.DoneSI(r);
			MultiID.DoneSI(q);
			DelPrimes(ref sNF, sPrimes);
			t = MultiID.TimerF();
			MultiID.ADebug("ReportFactor exit SieveMethod");
		} // SieveMethod

		private static void InitSieve(MultiSI n) // part of SieveMethod
		{
			int m;			// first rS primes
			int mM;			// (-x) % m
			bool match;		// match true or false for sieve table
			int ii;			// index
			int jj;			// index
			int kk;			// index
			int Dig;		// result of one digit Modulo function
			int tot = 0;
			MultiID.SetToSI(n, sieveN);
			MultiSI x = new MultiSI(size);			// x for search, factor = x - y
			MultiSI root = new MultiSI(size);		// square root temp
			MultiSI dP = new MultiSI(size);			// double precision temp
			for (ii = 1; ii <= rS; ii++)
			{
				if (MultiID.muErr || MultiID.muAbort) break;
				m = pTab[ii + 1];
				for (jj = 0; jj < m; jj++)
					squareMod[jj] = (int)((jj * jj) % m);
				for (jj = 0; jj < m; jj++)
				{
					if (MultiID.muErr || MultiID.muAbort) break;
					MultiID.Add1SI(n, -squareMod[jj], dP);	// dP = n - squareMod(jj)
					while (dP.S)
						MultiID.Add1SI(dP, m, dP);			// dP = dP + m
					MultiID.Mod1SI(m, out Dig, dP);			// need -dP % m
					mM = (int)(m - Dig);
					if (mM == m) mM = 0;
					match = false;
					kk = 0;
					while (kk < m)
					{
						if (mM == squareMod[kk])
						{
							match = true;
							kk = (int)m;
						}
						kk++;
					}
					sieve[ii, jj] = match;
					if (MultiID.statReq2)
					{
						MultiID.WriteT("I=" + ii + " J=" + jj + " MATCH=", ref tot);
						if (match)
							MultiID.WriteT("True", ref tot);
						else
						{
							MultiID.WriteT("False", ref tot);
						}
						MultiID.TimeDel(ref tot);
						RunForm.FocusBothUs();
						MultiID.statReq2 = false;
					}
				}
				if (MultiID.statReq2)
				{
					MultiID.WriteLn();
					RunForm.FocusBothUs();
					MultiID.statReq2 = false;
				}
			}
			MultiID.DoneSI(dP);
			MultiID.SqRtRemSI(n, remS, root);
			MultiID.SetToSI(root, x);
			if (MultiID.statReq2 || MultiID.debug)
			{
				MultiID.WriteT("SqRt=", ref tot);
				MultiID.WriteSI(x, ref tot);
				MultiID.WriteT(" Rem=", ref tot);
				MultiID.WriteSI(remS, ref tot, 0, false);
				MultiID.TimeDel(ref tot);
				if (MultiID.statReq2) RunForm.FocusBothUs();
				MultiID.statReq2 = false;
			}
			if (!MultiID.ZTestSI(remS)) MultiID.RAdd1SI(1, x);	// x = x + 1
			MultiID.Mul1SI(root, 2, root2);						// root2 = root * 2
			MultiID.SubSI(x, root, delX);						// delX = x - root(n)
			MultiID.DoneSI(root);
			for (ii = 1; ii <= rS; ii++)
			{
				m = pTab[ii + 1];
				MultiID.Mod1SI(m, out Dig, x);					// need -x % m
				mM = (int)(m - Dig);
				if (mM == m) mM = 0;
				if (MultiID.statReq2)
				{
					MultiID.WriteT("KS(" + ii + ")=" + mM, ref tot);
					MultiID.TimeDel(ref tot);
					MultiID.statReq2 = false;
				}
				kS[ii] = mM;
			}
			MultiID.DoneSI(x);
		} // InitSieve

		private static void FactorN(			// factor n into its prime factors
						MultiSI n,				// n
						MultiSI start,			// start	= Minimum factor to start prime test
						bool report,			// report	= true; to report factors to output screen
						out int fNF,			// fNF		= Number of prime factors found
						TPrimes fPrimes,		// fPrimes	= Array of prime factors found
						Tpowers fPowers)		// fPowers	= Array of powers of the prime factors found
		{
			bool found;				// factor found
			bool twoFactors;		// only two prime factors in number
			bool exitFactor;		// exit Factor flag
			int po = 0;				// power of prime factor
			int tot = 0;
			Application.DoEvents();
			MultiID.MuInterrupt();
			MultiID.ADebug("Enter FactorN");
			MultiSI div = new MultiSI(sizeO2);	// divisor to test
			if (MultiID.debug)
			{
				MultiID.WriteT("N= ", ref tot);
				MultiID.WriteSI(n, ref tot);
				MultiID.WriteT(" Start= ", ref tot);
				MultiID.WriteSILn(start, ref tot, 0, false);
			}
			fNF = 0;
			twoFactors = false;
			if (!MultiID.GT1UI(pTab[lenPTab] + 1, start))	// start < pTab[lenPTab] + 2
			{
				SmallFactors(n, report, ref fNF, fPrimes, fPowers, out exitFactor, div, ref po);
				if (exitFactor)
				{
					MultiID.ADebug("SmallFactors exit FactorN");
					MultiID.DoneSI(div);
					return;
				}
				MultiID.SetTo1SI(pTab[lenPTab] + 2, start);	// start = pTab[lenPTab] + 2
				if (fNF > 0 && report && findAPrime)
				{
					if ( (fNF == 1) && (fPowers.a[1] == 1) && MultiID.EQ1UI(1, n) )
					{
						MultiID.ADiag("N is prime");
						findAPrime = false;
					}
					else
					{
						MultiID.ADiag("N is not prime");
					}
					MultiID.ADebug("Quick exit FactorN");
					MultiID.DoneSI(div);
					return;
				}
				if (MultiID.GT1UI(1, n))
				{
					if (MultiID.diagOn)
						MultiID.WriteSILn(n, ref tot);
					MultiID.ADiag("has no small factors");
				}
			}
			MultiID.ClearSI(sieveN);
			found = true;
			while (MultiID.GT1UI(1, n))
			{
				do
				{
					if (found)
					{
						LargeFactor(n, report, out exitFactor, div, ref po, ref found);
						if (exitFactor)
						{
							MultiID.ADebug("LargeFactor exit FactorN");
							MultiID.DoneSI(div);
							return;
						}
					}
					MultiID.ADebug("In FactorN loop");
					if (found)
					{
						ReportFactor(report, ref fNF, ref fPrimes, fPowers, out exitFactor, div, po);
						if (exitFactor)
						{
							MultiID.ADebug("ReportFactor exit FactorN");
							MultiID.DoneSI(div);
							return;
						}
					}
					if (MultiID.GT1UI(1, n))
					{
						MoreDiv(n, start, report, ref fNF, fPrimes, fPowers, out exitFactor, div, ref po, ref found, ref twoFactors);
						if (exitFactor)
						{
							MultiID.ADebug("MoreDiv exit FactorN");
							MultiID.DoneSI(div);
							return;
						}
					}
				}
				while (found && !MultiID.EQ1UI(1, n));
				if (MultiID.GT1UI(1, n))
				{
					SieveMethod(n, start, report, ref fNF, fPrimes, fPowers, out exitFactor, div, ref po, ref found, twoFactors);
					if (exitFactor)
					{
						MultiID.ADebug("SieveMethod exit FactorN");
						MultiID.DoneSI(div);
						return;
					}
				}
			}
			MultiID.DoneSI(div);
			MultiID.ADebug("Exit FactorN");
		} // FactorN

		private static void NextMainNumber() // main operator loop
		{
			MultiID.ReSizeSI(size, z);
			if (MultiID.muAbort) findAPrime = false;
			MultiID.muErr = false;
			if (findAPrime)
			{
				if (!testPrime)
				{
					MultiID.SetTo1SI(2, start);
					if (findUp)
						MultiID.RAdd1SI(3, zIn);	// to increase
					MultiID.RDiv1SI(2, zIn);		// zIn = -((zIn / 2) * 2 - 1)
					MultiID.RMul1SI(2, zIn);
					MultiID.RAdd1SI(-1, zIn);
				}
				zIn.S = true;
			}
			if (!pTabFilled) FillPTab();
			findAPrime = zIn.S && MultiID.GT1UI(1, zIn);
			zIn.S = false;
			if (zIn.N > (MultiID.muDMax * size))
			{
				MultiID.WriteLn("N is too big");
				MultiID.ClearSI(z);
			}
			else
			{
				MultiID.SetToSI(zIn, z);
			}
			if (!MultiID.ZTestSI(z))
			{
				if (MultiID.GT1UI(1, z))
				{
					MultiID.ADebug("Starting to factor");
					FactorN(z, start, true, out nF, primes, powers);
					byECM = false;
					if (MultiID.GT1UI(1, z))
					{
						complete = false;
					}
					else
					{
						complete = true;
					}
				}
			}
		} // NextMainNumber

		private static void SigmaF()		// computing the Sigma(x) function
		{
			MultiSI pT = new MultiSI(MultiID.sdInDbl);	// a temp power
			MultiID.ReSizeSI(size + 4, z);
			MultiID.ReSizeSI(2 * size, zIn);
			MultiID.SetTo1SI(1, zIn);
			for (int i = 1; i <= nF; i++)
			{
				if (primes.a[i] == 2) // EQ1SI(2, primes.a[i])	// if p = 2 { z = z * (2^(a+1) - 1)
				{
					MultiID.SetToDbSI(powers.a[i] + 1, pT);		// pT = a+1
					MultiID.PowSI(primes.a[i], pT, z);			// z = 2^(a+1)
					MultiID.RAdd1SI(-1, z);						// z = z-1
					MultiID.RMulSI(z, zIn);						// zIn= zIn * z
				}
				else if (powers.a[i] == 1)						// if a = 1 { z = z * (p+1)
				{
					MultiID.RAdd1SI(1, primes.a[i]);			// p = p+1
					MultiID.RMulSI(primes.a[i], zIn);			// zIn= zIn * (p+1)
					MultiID.RAdd1SI(-1, primes.a[i]);			// p=p-1
				}
				else											// else z = z * (p^(a+1) - 1) / (p-1)
				{
					MultiID.SetToDbSI(powers.a[i] + 1, pT);		// pT = a+1
					MultiID.PowSI(primes.a[i], pT, z);			// z = p^(a+1)
					MultiID.RAdd1SI(-1, z);						// z=z-1
					MultiID.RAdd1SI(-1, primes.a[i]);			// p=p-1
					MultiID.RDivSI(primes.a[i], z);				// z=(p^(a+1) - 1) / (p-1)
					MultiID.RAdd1SI(1, primes.a[i]);			// p=p+1
					MultiID.RMulSI(z, zIn);						// zIn= zIn * z
				}
			}
			MultiID.DoneSI(pT);
		} // SigmaF

		public static void PrimeSig(int nI)	// n = Sig(n)
		{
			MultiID.ADDiag("Prime Sig function = Sum of the divisors of x:");
			if (MultiID.EQ1UI(1, Common.nList[nI]))				// +/-1 => 1
			{
				Common.nList[nI].S = false;
				DelPrimes(ref nF, primes);
				complete = false;
				return;
			}
			if (MultiID.ZTestSI(Common.nList[nI]) && !complete)
				return;					// 0 => 0
			ECM(nI);					// factor n by ECM
			if (complete)
			{
				SigmaF();
				MultiID.SetToSI(zIn, Common.nList[nI]);		// n = Sig(n) (Including n)
				Calc.NormN(nI);
			}
			Fini2();
		}

		public static void PrimeSig0(int nI) // n = Sig(n) - n
		{
			MultiID.ADDiag("Prime Sig0 function = Sum of the divisors of x (-x):");
			if (MultiID.EQ1UI(1, Common.nList[nI]))			// +/-1 => 0
			{
				Calc.ClearN(nI);
				DelPrimes(ref nF, primes);
				complete = false;
				return;
			}
			if (MultiID.ZTestSI(Common.nList[nI]))
			{
				if (!complete) return;						// 0 => 0
				PrimeFAInv(Common.nList[nI]);				// reconstruct n from factors
				if (MultiID.diagOn)
				{
					int tot = 0;
					MultiID.WriteSI(Common.nList[nI], ref tot);
					MultiID.WriteT(" already factored with " + nF + " prime(s)", ref tot);
				}
			}
			else
			{
				ECM(nI);									// factor n by ECM
			}
			if (complete)
			{
				SigmaF();
				Common.nList[nI].S = false;
				MultiID.RSubSI(Common.nList[nI], zIn);		// n = Sig(n) - n (not including n)
				MultiID.SetToSI(zIn, Common.nList[nI]);
				Calc.NormN(nI);
			}
			Fini2();
		}

		private static void TauF() // computing the Tau(x) function
		{
			MultiID.ReSizeSI(2, z);
			MultiID.ReSizeSI(size, zIn);
			MultiID.SetTo1SI(1, zIn);
			for (int i = 1; i <= nF; i++)						// Tau(n) = Product (a+1)
			{
				MultiID.SetToDbSI(powers.a[i] + 1, z);			// z = a+1
				MultiID.RMulSI(z, zIn);							// zIn= zIn * z
			}
		}

		public static void PrimeTau(int nI) // n = Tau(n)
		{
			MultiID.ADDiag("Prime Tau function = Number of divisors of x:");
			if (MultiID.EQ1UI(1, Common.nList[nI]))	// +/-1 => 1
			{
				Common.nList[nI].S = false;
				DelPrimes(ref nF, primes);
				complete = false;
				return;
			}
			if (MultiID.ZTestSI(Common.nList[nI]) && !complete)
				return;							// 0 => 0
			ECM(nI);							// factor n by ECM
			if (complete)
			{
				TauF();
				MultiID.SetToSI(zIn, Common.nList[nI]);	// n = Tau(n)
				Calc.NormN(nI);
			}
			Fini2();
		}

		public static void MoebiusMu(int nI)
		// Mu(X) = Moebius Mu(X) function
		{
			MultiID.ADDiag("Moebius Mu(X) function:");
			if (MultiID.EQ1UI(1, Common.nList[nI]))	// +/-1 => 1
			{
				Common.nList[nI].S = false;
				DelPrimes(ref nF, primes);
				complete = false;
				return;
			}
			if (MultiID.ZTestSI(Common.nList[nI]))
			{
				if (byECM && !squareFree && !complete)
					MultiID.ADiag("Stored Number not squarefree");
				if (!complete) return;				// 0 => 0
				if (MultiID.diagOn)
				{
					PrimeFAInv(z);
					int tot = 0;
					MultiID.WriteSI(z, ref tot);
					MultiID.WriteT(" already factored with " + nF + " prime(s)", ref tot);
					MultiID.DoneSI(z);
				}
			}
			else
			{
				testSqFree = true;
				squareFree = true;
				ECM(nI);							// factor n by ECM
				testSqFree = false;
			}
			int mu = 0;
			if (!byECM && complete)
			{
				squareFree = true;
				for (int i = 1; i <= nF; i++)
				{
					if (powers.a[i] != 1)
					{
						squareFree = false;
						break;
					}
				}
			}
			if (squareFree && complete)
			{
				mu = 1;
				if ((nF & 1) == 1)				// if an odd number of non-repeating factors
					mu = -1;					// final mu is -1
			}
			MultiID.SetTo1SI(mu, Common.nList[nI]);
			Calc.NormN(nI);
			Fini2();
		} // MoebiusMu

		public static void IsSq(int nI)
		// set n true if n is a perfect square, else false
		{
			MultiID.ADDiag("IsSq(X) function:");
			if (MultiID.EQ1UI(1, Common.nList[nI]))	// +/-1 => 1 (true) (sign not used)
			{
				Common.nList[nI].S = false;
				DelPrimes(ref nF, primes);
				complete = false;
				return;
			}
			if (MultiID.ZTestSI(Common.nList[nI]))
			{
				if (!complete)
				{
					Common.nList[nI].V[0] = 1;	// 0 => 1 (true)
					return;
				}
				if (MultiID.diagOn)
				{
					PrimeFAInv(z);
					int tot = 0;
					MultiID.WriteSI(z, ref tot);
					MultiID.WriteT(" already factored with " + nF + " prime(s)", ref tot);
					MultiID.DoneSI(z);
				}
			}
			else
			{
				testSquare = true;
				isSquare = true;
				ECM(nI);						// factor n by ECM
				testSquare = false;
			}
			Calc.ClearN(nI);					// isSquare = 0 (false)
			if (!byECM) isSquare = true;
			if (isSquare && complete)		 	// isSquare not proven
			{
				Common.nList[nI].V[0] = 1;		// isSquare = 1 (true)
				for (int i = 1; i <= nF; i++)
				{
					if ((powers.a[i] & 1) == 1)	// if any power is odd
					{
						Common.nList[nI].V[0] = 0;	// final isSquare = 0 (false)
						break;
					}
				}
			}
			Calc.NormN(nI);
			Fini2();
		} // IsSq

		public static void SqFree(int nI)
		// set n true if n is a square free, else false
		{
			MultiID.ADDiag("SqFree(X) function:");
			MoebiusMu(nI);
			Common.nList[nI].S = false;
		} // SqFree

		public static void PrimePhi(int nI) // n = Phi(n)
		{
			PrimePhi(Common.nList[nI]);
			Calc.NormN(nI);
		}

		public static void PrimePhi(MultiSI n) // n = Phi(n) = Euler's totient function
		{
			MultiID.ADDiag("Prime Phi function = Number of integers <= n relatively prime to n:");
			if (MultiID.EQ1UI(1, n))			// +/-1 => 1
			{
				n.S = false;
				DelPrimes(ref nF, primes);
				complete = false;
				return;
			}
			if (MultiID.ZTestSI(n) && !complete)
				return;							// 0 => 0
			ECM(n);								// factor n by ECM
			if (complete)								// Phi(n) = Product (p-1)p^(a-1)
			{
				MultiSI pT = new MultiSI(MultiID.sdInDbl);		// a temp power
				MultiID.ReSizeSI(size, z);
				MultiID.ReSizeSI(size, zIn);
				MultiID.SetTo1SI(1, zIn);
				for (int i = 1; i <= nF; i++)
				{
					MultiID.SetToDbSI(powers.a[i] - 1, pT);		// pT = a-1
					MultiID.PowSI(primes.a[i], pT, z);			// z = p^(a-1)
					MultiID.RAdd1SI(-1, primes.a[i]);			// p=p-1
					MultiID.RMulSI(primes.a[i], z);				// z= z * p
					MultiID.RAdd1SI(1, primes.a[i]);			// p=p+1
					MultiID.RMulSI(z, zIn);						// zIn= zIn * z
				}
				MultiID.SetToSI(zIn, n);						// n = Tau(n)
				MultiID.DoneSI(pT);
			}
			Fini2();
		} // PrimePhi

		public static void PrimeAllD(int nI)	// find all divisors of n
		{
			int[] count;		// count of the factors
			int i;
			bool done;
			int baseS;
			int tot = 0;
			if (MultiID.EQ1UI(1, Common.nList[nI]))
			{
				DelPrimes(ref nF, primes);
				complete = false;
				MultiID.WriteDLn2("All Divisors: 1.");
				return;
			}
			if (MultiID.ZTestSI(Common.nList[nI]) && !complete)
			{
				MultiID.WriteDLn2("All Divisors: 0 has no divisors.");
				return;
			}
			ECM(nI);					// factor n by ECM
			if (complete)				// Phi(n) = Product (p-1)p^(a-1)
			{
				baseS = MultiID.baseO;
				MultiID.baseO = 10;
				TauF();
				Common.DoubleSpace();
				MultiID.WriteT("All ", ref tot);
				MultiID.WriteSI(zIn, ref tot, 0, false);	// display number of divisors in decimal
				MultiID.WriteT(" Divisors: ", ref tot);
				MultiID.baseO = baseS;
				MultiSI pT = new MultiSI(MultiID.sdInDbl);	// a temp power
				MultiID.ReSizeSI(size, z);
				MultiID.ReSizeSI(size, zIn);
				if (nF == 0) nF++;
				count = new int[nF+1];						//ReDim count(1 To nF)
				for (i = 2; i <= nF; i++)
					count[i] = 0;
				count[1] = 1;
				done = false;
				MultiID.WriteT("1", ref tot);
				do
				{
					MultiID.SetTo1SI(1, zIn);
					for (i = 1; i <= nF; i++)
					{
						MultiID.SetToDbSI(count[i], pT);	// pT = C
						MultiID.PowSI(primes.a[i], pT, z);	// z = p^C
						MultiID.RMulSI(z, zIn);				// zIn= zIn * z
					}
					AddToCounts(count, 1, ref done);
					MultiID.WriteT("; ", ref tot);
					MultiID.WriteSI(zIn, ref tot, 0, false);
				}
				while (!done && !MultiID.muAbort);
				if (done) MultiID.WriteT(".", ref tot);		// end of divisors
				if (MultiID.baseO != 10) MultiID.WriteT(" (base " + MultiID.baseO + ")", ref tot);
				Common.DoubleSpace();
				MultiID.DoneSI(pT);
			}
			Fini2();
			MultiID.muErr = false; // because no value is returned to nList[nI]
		} // PrimeAllD

		public static void AddToCounts(int[] count, int i, ref bool done) // add 1 to count array
		{
			count[i]++;
			if (count[i] > powers.a[i])
			{
				count[i] = 0;
				if (i == nF)
				{
					done = true;
					return;
				}
				AddToCounts(count, i + 1, ref done);
			}
		}

		public static void PrimeSearch(int nI) // prime Search
		{
			DelPrimes(ref nF, primes);
			complete = false;
			if (findUp)
			{
				MultiID.ADDiag("Prime Greater Than:");
				if (!MultiID.GT1UI(1, Common.nList[nI]))	// 0,1 => 2
				{
					MultiID.SetTo1SI(2, Common.nList[nI]);
					return;
				}
			}
			else
			{
				MultiID.ADDiag("Prime Less Than:");
				if (!MultiID.GT1UI(2, Common.nList[nI]))	// 0,1,2 => 0
				{
					Calc.ClearN(nI);
					return;
				}
				if (!MultiID.GT1UI(3, Common.nList[nI]))	// 3 => 2
				{
					MultiID.SetTo1SI(2, Common.nList[nI]);
					return;
				}
			}
			MultiID.ReSizeSI(Common.nList[nI].N + 1, zIn);	// allow for it getting 1 SD bigger
			MultiID.SetToSI(Common.nList[nI], zIn);
			zIn.S = false;
			MultiSI d = new MultiSI(1);
			d.V[0] = 2;								// d = 2
			if (findUp) MultiID.RAdd1SI(3, zIn);	// to increase
			else d.S = true;						// d = -2
			MultiID.RDiv1SI(2, zIn);				// zIn = -((zIn / 2) * 2 - 1)
			MultiID.RMul1SI(2, zIn);				// next odd number > or < n
			MultiID.RAdd1SI(-1, zIn);
			while (FnAdleman(zIn) == 0)				// while not prime
			{
				if (MultiID.muErr || MultiID.muAbort)
				{
					MultiID.ClearSI(zIn);
					break;
				}
				MultiID.RAddSI(d, zIn);
			}
			MultiID.DoneSI(d);
			MultiID.SetToSI(zIn, Common.nList[nI]);
			Calc.NormN(nI);
		}

		private static void PrimeFAInv(MultiSI z)
		// prime factor algorithm A inverse z=Pi(p^e), ^ is used here to mean p to the power e, not exclusive or
		{
			MultiSI e = new MultiSI(MultiID.sdInDbl);
			MultiSI t = new MultiSI(size);
			MultiID.ReSizeSI(size, z);
			MultiID.SetTo1SI(1, z);
			for (int i = 1; i <= nF; i++)
			{
				MultiID.SetToDbSI(powers.a[i], e);
				MultiID.PowSI(primes.a[i], e, t);		// t = p^e
				MultiID.RMulSI(t, z);					// z= z * p^e
			}
			MultiID.DoneSI(t);
			MultiID.DoneSI(e);
		}

		public static void PrimeFAA(int nI) // prime factor algorithm A
		{
			MultiID.ADDiag("Prime factor algorithm A:");
			PrimeF(nI, "A");
		}

		public static void PrimeFAE(int nI) // prime factor algorithm E, Elliptic Curve Method (ECM)
		{
			MultiID.ADDiag("Prime factor algorithm E (ECM):");
			PrimeF(nI, "E");
		}

		private static void PrimeF(int nI, string st) // prime factor algorithm A or E
		{
			if (MultiID.EQ1UI(1, Common.nList[nI]))
			{
				MultiID.WriteDLn("1 = 1^1.");
				DelPrimes(ref nF, primes);
				complete = false;
				return;
			}
			findAPrime = false;
			testPrime = false;
			if (MultiID.ZTestSI(Common.nList[nI]))
			{
				if (!complete)
				{
					MultiID.WriteDLn("0 = 0.");
					return;
				}
				if (MultiID.diagOn)
				{
					PrimeFAInv(z);
					int tot = 0;
					MultiID.WriteSI(z, ref tot);
					MultiID.WriteT(" already factored with " + nF + " prime(s)", ref tot);
					MultiID.DoneSI(z);
				}
				MultiID.DoneSI(zIn);					// clear zIn
				MultiID.ReSizeSI(size, z);				// size saved from previous PrimeDoIt call
				MultiID.SetTo1SI(1, z);
			}
			else
			{
				if (st == "A") PrimeDoIt(nI);
				else if (st == "E") ECM(nI);
			}
			if (complete || exitEarly) DisplayFactors();
			Fini2();
			MultiID.muErr = false; // because no value is returned to nList[nI]
		} // PrimeF

		private static void DisplayFactors() // display results
		{
			int sign;	// sign result of ABCompUI Sub
			int tot = 0;
			if (MultiID.diagOn && !MultiID.muErr && !MultiID.muAbort && !exitEarly)
			{
				Common.SingleSpace();
				MultiID.WriteT("Completely factored", ref tot);
				if (MultiID.baseO != 10) MultiID.WriteT(" (Exponents are in decimal)", ref tot);
				MultiID.WriteTLn("", ref tot);
			}
			PrimeFAInv(z);
			Common.DoubleSpace();
			if (MultiID.ZTestSI(zIn))
			{
				sign = 0;
				MultiID.WriteSI(z, ref tot);
			}
			else
			{
				sign = MultiID.AbCompUI(zIn, z);
				MultiID.WriteSI(zIn, ref tot);
			}
			// sign = 1; // diag
			if (sign == 0) MultiID.WriteT(" = ", ref tot);
			if (sign != 0) MultiID.WriteT(" != ", ref tot);
			for (int i = 1; i <= nF; i++)
			{
				if (i != 1) MultiID.WriteT(" * ", ref tot);
				MultiID.WriteSI(primes.a[i], ref tot, 0, false);
				if (powers.a[i] != 1)
					MultiID.WriteT("^" + powers.a[i], ref tot);
			}
			if (sign != 0) // z != zIn
			{
				MultiID.WriteTLn(" =", ref tot);
				MultiID.WriteSILn(z, ref tot, 0, false);
				MultiID.Write("Error, factors do not check");
				DelPrimes(ref nF, primes);
				complete = false;
			}
			else
			{
				MultiID.WriteT(".", ref tot); // end of a good factoring
			}
			Common.DoubleSpace();
		} // DisplayFactors

		private static void InsertFactor(MultiSI d)	// Insert factor in factor arrayCurve Method
		{
			if (nF > 0 && primes.a[nF] == d)		// old factor
			{
				powers.a[nF]++;
				squareFree = false;
				return;
			}
			if (nF < maxPs)							// new factor
			{
				if (nF == 0)						// if first factor
				{
					isSquare = true;
					squareFree = true;
					complete = false;
				}
				else
				{
					if ((powers.a[nF] & 1) == 1)	// if found an odd power
					{
						isSquare = false;
					}
				}
				nF++;
				primes.a[nF] = new MultiSI(d.N);
				MultiID.SetToSI(d, primes.a[nF]);
				powers.a[nF] = 1;
			}
			else
			{
				MultiID.WriteSLn("Error, cannot store more than " + maxPs + " prime factors");
				exitEarly = true;
			}
		} // InsertFactor

		private static void ECM(int nI) // Prime Factorization by ECM, Elliptic Curve Method
		// from UBASIC program emc.ub, Prime Factorization by ECM, 1987-1990 by Yuji KIDA
		{
			ECM(Common.nList[nI]);
		}

		public static void ECM(MultiSI n) // Prime Factorization by ECM, Elliptic Curve Method
		// from UBASIC program emc.ub, Prime Factorization by ECM, 1987-1990 by Yuji KIDA
		{
			int tot = 0;
			if (MultiID.ZTestSI(n))
			{
				if (!complete) return;				// 0 => 0
				if (MultiID.diagOn)
				{
					PrimeFAInv(z);
					MultiID.WriteSI(z, ref tot);
					MultiID.WriteT(" already factored with " + nF + " prime(s)", ref tot);
					MultiID.DoneSI(z);
				}
				return;
			}
			int i;
			bool swb;
			if (!pTabFilled) FillPTab();
			Common.SingleSpace();
			complete = false;
			exitEarly = false;
			int nD = 100; // was 21
			MultiSI[] ss = new MultiSI[nD];
			MultiSI[] pd = new MultiSI[nD];
			for (i = 0; i < nD; i++)
			{
				ss[i] = new MultiSI();
				pd[i] = new MultiSI();
			}
			size = n.N;
			size = Math.Min(size, MultiID.mMax);
			size = Math.Max(size, (48 / MultiID.muDMax) + 1);	// size for at least 49 decimal digits
			Init1();
			MultiID.SetToSI(n, zIn); // zIn sized in Init1
			zIn.S = false;
			MultiID.ReSizeSI(size, z);
			int s1i = 1;
			int spi = 0;
			int dpi = 0;
			int exN = MultiID.sdInLong;
			MultiID.SetToSI(zIn, z);		// note: zIn = NN, z = N, zz = Z
			InitP();
			long dX = 2;					// Divide by small primes
			MultiSI dSI = new MultiSI(exN);
			stopDiv = 131072;	// 2^17 for PrimeDiv function (2^17 from UBASIC program)
			while (true)
			{
				if (MultiID.muErr || MultiID.muAbort) goto ExitECM; // L200:
				PrimeDiv(z, ref dX);
				if (dX == 0) break;
				MultiID.SetToExSI(dX, dSI);
				InsertFactor(dSI);				// insert small prime factors
				if (z == dSI) break;
				if (exitEarly) goto ExitECM;
				MultiID.DiagT(dX + " * ", ref tot);
				if ((testSqFree && !squareFree) || (testSquare && !isSquare))
					goto ExitECM;
				MultiID.RDivSI(dSI, z);			// z = z / dX
				s1i = 0;
			}
			MultiID.DoneSI(dSI);
			if (dX != 0) // L240:
			{
				complete = true;				// completely factored
				if (s1i != 0)
				{
					MultiID.DiagT(dX + " is prime", ref tot);
					goto ExitECM;
				}
				MultiID.DiagT("" + dX + MultiID.nL, ref tot);
				goto ExitECM;
			}
			if (nF > 0 && testPrimR)	// if testing for a primitive root and there are two odd prime factors
			{
				nF = 2;
				goto ExitECM;			// the second factor may be null
			}
			tot = 0;
			Common.SingleSpace();		// When N has a big factor // L290:
			MultiSI n1 = new MultiSI(size);
			MultiSI n2 = new MultiSI(size);
			while (true)
			{
				// ADLEMAN TEST // L330:
				if (MultiID.muErr || MultiID.muAbort) goto ExitECM;
				if (FnAdleman(z) == 0)			// returns 1 if z is prime, 0 if composite
				{
					if(MultiID.muErr) goto ExitECM;
					s1i = 0;
				}
				else
				{
					MultiID.SetToSI(z, pd[dpi]);	// Store in PD if prime
					dpi++;
					if (spi != 0)					// if not finished
					{
						if (dpi > 1 && testPrimR && pd[dpi-2] != pd[dpi-1])
						{
							nF = 2;
							goto ExitECM; // if testing for a primitive root and there are two odd prime factors
						}
						spi--;
						MultiID.SetToSI(ss[spi], z);
						continue;
					}
					break; // end while
				}
				// Check if N is a power of an integer // L430:
				FnPowerCheck(z, n1);
				MultiID.DivSI(z, n1, n2);		// n2 = z / n1
				if (n1 != 1)
				{
					if (n1 > n2) MultiID.SwapSI(ref n1, ref n2);
					MultiID.SetToSI(n2, ss[spi]);
					spi++;
					MultiID.SetToSI(n1, z);
					continue;
				}
				FnECM(z, n1);					// Find a factor by ECM // L500:
				MultiID.DivSI(z, n1, n2);		// n2 = z / n1
				if (n1 > n2) MultiID.SwapSI(ref n1, ref n2);
				MultiID.SetToSI(n2, ss[spi]);
				spi++;
				MultiID.SetToSI(n1, z);
			}
			MultiID.DoneSI(n2);
			MultiID.DoneSI(n1);
			if (s1i != 0)					// Display the result // L560:
			{
				complete = true;			// completely factored
				InsertFactor(z);			// insert the one and only prime factor
				if (MultiID.diagOn) MultiID.WriteSI(z, ref tot, 0, false);
				MultiID.DiagT(" is prime", ref tot);
				complete = true;
				goto ExitECM;
			}
			tot = 0;
			Common.SingleSpace();
			if (dpi == 1)
			{
				InsertFactor(pd[0]);		// insert the only one large prime factor
				if (exitEarly) goto ExitECM;
				if (MultiID.diagOn) MultiID.WriteSILn(pd[0], ref tot, 0, false);
				complete = true;			// completely factored
				goto ExitECM;
			}
			do								// Sort and display the larger factors
			{
				swb = false;
				for (i = 0; i <= dpi - 2; i++)
				{
					if (pd[i] > pd[i + 1])
					{
						MultiID.SwapSI(ref pd[i], ref pd[i+1]); // only pointers are swapped
						swb = true;
					}
				}
			} while (swb);
			for (i = 0; i <= dpi - 2; i++)
			{
				InsertFactor(pd[i]);		// insert large prime factors
				if (exitEarly) goto ExitECM;
				if (MultiID.diagOn) MultiID.WriteSI(pd[i], ref tot, 0, false);
				MultiID.DiagT(" * ", ref tot);
			}
			InsertFactor(pd[dpi-1]);		// insert final large prime factor
			if (exitEarly) goto ExitECM;
			if (MultiID.diagOn) MultiID.WriteSI(pd[dpi-1], ref tot, 0, false);
			complete = true;				// completely factored
		ExitECM:
			pd = new MultiSI[0];
			ss = new MultiSI[0];
			if (!complete && !Common.quiet)
				MultiID.WriteLn("Number not completely factored");
			byECM = true;
			return;
		} // ECM

		private static int FnAdleman(MultiSI z) // Adleman function, returns 1 if z is prime, 0 if composite
		// from UBASIC program emc.ub, Prime Factorization by ECM, 1987-1990 by Yuji KIDA
		{
			int returnV;
			if (MultiID.muErr) return 0;
			if (!pTabFilled) FillPTab();
			int tot = 0;
			MultiID.DiagT("Adleman Test for ", ref tot);
			if (MultiID.diagOn) MultiID.WriteSILn(z, ref tot);
			if (MultiID.ZTestSI(z)) return 0;
			long L = 1;			// reset NextP()
			NextP(ref L);		// NextP will call InitP() if needed and set xI = 1
			int i;
			int pWords = 2 * z.N + 2;
			int ei,eti,gi,hi,ii,ji,mwi,npi,nqi,pi,qi,wi,wki,xi;
			bool ok, swb;
			qi = eti = 0;	// to init ref items
			int[] indexa = new int[9283];		// index mod Q with base G
			int[] phia = new int[9283];
			int[] psia = new int[9283];
			int aSize = 17;
			MultiSI[] tau  = new MultiSI[aSize]; // 16=max of P-1
			MultiSI[] tau1 = new MultiSI[aSize];
			MultiSI[] tau2 = new MultiSI[aSize];
			MultiSI[] taue = new MultiSI[aSize];
			MultiSI[] sum  = new MultiSI[aSize];
			MultiSI[] taup = new MultiSI[aSize];
			for (i = 0; i < aSize; i++)
			{
				tau[i]  = new MultiSI(pWords);
				tau1[i] = new MultiSI(pWords);
				tau2[i] = new MultiSI(pWords);
				taue[i] = new MultiSI(pWords);
				sum[i]  = new MultiSI(pWords);
				taup[i] = new MultiSI(pWords);
			}
			int[,] wka = new int[7, 28];		// 6=no of P -1,27=no of Q -1
			int[,] eta = new int[7, 28];
			int[] pa = { 2, 3, 5, 7, 11, 13, 17 };	// 0...6 PrimeData1
			int[] qa = {3, 7, 11, 31, 43, 71, 211, 23, 67, 331,
						463, 2311, 79, 131, 547, 859, 911, 2003, 2731, 6007,
						103, 239, 443, 1123, 1327, 1871, 3571, 9283 };	// 0...27
			int[] ga = {2, 3, 2, 3, 3, 7, 2, 5, 2, 3,
						3, 3, 3, 2, 2, 2, 17, 5, 3, 3,
						5, 7, 2, 2, 3, 14, 2, 2 };		// 0...27
			MultiSI s  = new MultiSI(pWords);
			MultiID.SetTo1SI(1, s);
			MultiSI ss = new MultiSI(pWords);	// to hold s^2
			for (nqi = 0; nqi < 28; nqi++)
			{
				MultiID.Mul1SI(s, qa[nqi], s);
				MultiID.SqSI(s, ss);
				if (ss > z) break;
			}
			if (nqi >= 28)
			{
				Common.DoubleSpace();
				MultiID.WriteLn("FnAdleman: I quit! Too big.");
				MultiID.DiagT("Max with large factors = ", ref tot);
				if (MultiID.diagOn) MultiID.WriteSILn(ss, ref tot);
				MultiID.muErr = true;
				complete = false;
				returnV = 0;
				goto ExitFnAdleman;
			}
			MultiID.DoneSI(ss);
			npi = 2; // L1130:
			if (nqi > 19) npi += 4;
			else if (nqi > 11) npi += 3;
			else if (nqi > 6) npi += 2;
			else if (nqi > 3) npi += 1;
			stopDiv = 131072;	// 2^17 for PrimeDiv function (2^17 from UBASIC program)
			long wL = 0;
			PrimeDiv(z, ref wL);
			MultiSI w  = new MultiSI(pWords);
			w &= SI(wL);
			if (wL != 0)
			{								// FactorQ:
				returnV = 0;
				if (w == z) returnV = 1;	// All OK
				goto ExitFnAdleman;
			}
			MultiSI up = new MultiSI(pWords);
			MultiSI wp = new MultiSI(pWords);
			Test4_2(out hi, up, wp, z);	// TEST FOR P=2
			if (hi == 1) swb = true;
			else
			{
				Test3_2(nqi, out ok, qa, ref qi, z);
				swb = ok;
			}
				mwi = 1;
			MultiSI t  = new MultiSI(pWords);
			MultiSI u  = new MultiSI(pWords);
			for (ji = 0; ji <= nqi; ji++)
			{
				qi = qa[ji];
				gi = ga[ji];
				Test2_2(out eti, out ok, qi, u, w, z);
				if (!ok)
				{
					returnV = 0; // Failed
					goto ExitFnAdleman;
				}
				if (swb) wki = hi;
				else
				{
					Test5_2(out eti, out ok, qi, t, up, w, out wki, z);
					if (!ok)
					{
						returnV = 0; // Get Factor
						goto ExitFnAdleman;
					}
				}
				wka[0, ji] = wki;
				eta[0, ji] = eti;
				if (wki > mwi) mwi = wki;
			}
			MultiID.DoneSI(u);
			if (mwi != 1)						// decide max of W and eta(chi)
			{
				for (ji = 0; ji <= nqi; ji++)
				{
					if (mwi > wka[0, ji]) eta[0, ji] = 0;
				}
			}
			MultiSI e  = new MultiSI(pWords);
			for (ii = 1; ii <= npi; ii++)		// TEST for P>2
			{
				pi = pa[ii];
				Test4_P(out hi, pi, up, wp, z);
				if (hi == 1) swb = true;
				else
				{
					Test3_P(out ei, nqi, out ok, pi, qa, ref qi, z);
					swb = ok;
				}
				mwi = 1;
				for (ji = 0; ji <= nqi; ji++)
				{
					qi = qa[ji];
					gi = ga[ji];
					do	// do once
					{
						if ((qi - 1) % pi != 0)
						{
							wki = 0;
							eti = 0;
							break;
						}
						Test2_P(e, out ei, gi, indexa, out ok, phia, psia, pi, qi, sum, swb, tau, tau1, tau2, taue, taup, w, out wi, z);
						if (!ok)
						{
							returnV = 0; // Failed
							goto ExitFnAdleman;
						}
						if (swb)
						{
							wki = hi;
							eti = indexa[z % qi] % pi;
							break;
						}
						Test5_P(ref eti, out ok, pi, t, tau, tau1, tau2, taue, taup, up, w, ref wi, out wki, z);
						if (!ok)
						{
							returnV = 0; // Failed
							goto ExitFnAdleman;
						}
					} while (false);
					wka[ii, ji] = wki; // L1550:
					eta[ii, ji] = eti;
					if (wki > mwi) mwi = wki;
				}
				if (mwi != 1)					// decide max of W and eta(chi)
				{
					for (ji = 0; ji <= nqi; ji++)
					{
						wki = wka[ii, ji];
						if(wki != 0 && mwi > wki)
						eta[ii, ji] = 0;
					}
				}
			}
			// final test
			MultiID.DoneSI(up);
			MultiID.DoneSI(e);
			MultiID.DoneSI(wp);
			MultiSI iN = new MultiSI(pWords);
			MultiSI m1 = new MultiSI(pWords);
			MultiSI p1 = new MultiSI(pWords);
			MultiSI u1 = new MultiSI(pWords);
			MultiSI v1 = new MultiSI(pWords);
			MultiSI v = new MultiSI(pWords);
			GET_V(iN, eta, m1, npi, nqi, p1, ga, pa, out pi, qa, ref qi, u1, v, v1, wka, out xi);
			MultiID.DoneSI(v1);
			MultiID.DoneSI(u1);
			MultiID.DoneSI(p1);
			MultiID.DoneSI(m1);
			MultiID.DoneSI(iN);
			w &= 1;
			t &= 1;
			for (ii = 0; ii <= npi; ii++) t *= pa[ii];
			returnV = 1;			// AllOK
			while (t != 0)
			{
				w &= (w * v) % s;
				if (w == 1)
				{
					break;			// AllOK
				}
				if (z % w == 0)
				{
					returnV = 0;	// Got Factor
					break;
				}
				t--;
			}
			MultiID.DoneSI(w);
			MultiID.DoneSI(v);
			MultiID.DoneSI(t);
			MultiID.DoneSI(s);
		ExitFnAdleman:
			tau  = new MultiSI[0];
			tau1 = new MultiSI[0];
			tau2 = new MultiSI[0];
			taue = new MultiSI[0];
			sum  = new MultiSI[0];
			taup = new MultiSI[0];
			return returnV;
		} // FnAdleman

		private static void Test3_2( // Subroutines for P=2, Check condition 3 for P=2
			int nqi,
			out bool ok,
			int[] qa,
			ref int qi,
			MultiSI z
			)
		{
			if (MultiID.debug)
			{
				int tot = 0;
				string st = " In Test3_2: " + "nqi=" + nqi + " ok=" + ".." + " qi=" + qi + " z =";
				MultiID.WriteT(st, ref tot);
				MultiID.WriteSILn(z, ref tot, 0, false);
			}
			int ji;
			ok = false;
			for (ji = 0; ji <= nqi; ji++)
			{
				qi = qa[ji];
					if (Kronecker(z, qi) != 1)
				{
					ok = true;
					break;
				}
			}
			if (MultiID.debug) // L2000:
			{
				int tot = 0;
				string st = "Out Test3_2: " + "nqi=" + nqi + " ok=" + ok + " qi=" + qi + " z =";
				MultiID.WriteT(st, ref tot);
				MultiID.WriteSILn(z, ref tot, 0, false);
			}
		} // Test3_2

		private static void Test4_2( // Check condition 4 for P=2
			out int hi,
			MultiSI up,
			MultiSI wp,
			MultiSI z
			)
		{
			MultiID.Add1SI(z, -1, wp);
			hi = -1;
			do
			{
				up &= wp;
				wp /= 2;
				hi++;
			}
			while (MultiID.rem1 == 0);
		} // Test4_2

		private static void Test2_2( // Check condition 2 for P=2
			out int eti,
			out bool ok,
			int qi,
			MultiSI u,
  			MultiSI w,
			MultiSI z
			)
		{
			if (MultiID.debug)
			{
				int tot = 0;
				string st = " In Test2_2: " + "eti=" + ".." + " ok=" + ".." + " qi=" + qi;
				MultiID.WriteT(st, ref tot);
				MultiID.WriteT(" u=", ref tot);
				MultiID.WriteSI(u, ref tot, 0, false);
				MultiID.WriteT(" w=", ref tot);
				MultiID.WriteSI(w, ref tot, 0, false);
				MultiID.WriteT(" z=", ref tot);
				MultiID.WriteSILn(z, ref tot, 0, false);
			}
			MultiID.SetTo1SI(qi, w);
			if ((qi - 1) % 4 != 0) MultiID.SubSI(z, w, w);
			u &= MultiID.ModPow(w, (z - 1) / 2, z);
			if (Kronecker(z, qi) < 0)
			{
				eti = 1;
				ok = (u == (z - 1));
			}
			else
			{
				eti = 0;
				ok = (u == 1);
			}
			if (MultiID.debug)
			{
				int tot = 0;
				string st = "Out Test2_2: " + "eti=" + eti + " ok=" + ok + " qi=" + qi;
				MultiID.WriteT(st, ref tot);
				MultiID.WriteT(" u=", ref tot);
				MultiID.WriteSI(u, ref tot, 0, false);
				MultiID.WriteT(" w=", ref tot);
				MultiID.WriteSI(w, ref tot, 0, false);
				MultiID.WriteT(" z=", ref tot);
				MultiID.WriteSILn(z, ref tot, 0, false);
			}
		} // Test2_2

		private static void Test5_2( // Check condition 5 for P=2
			out int eti,
			out bool ok,
			int qi,
			MultiSI t,
			MultiSI up,
  			MultiSI w,
			out int wki,
			MultiSI z
			)
		{
			MultiSI temp = new MultiSI(z.N);
			MultiID.SetTo1SI(Kronecker(-1, qi) * qi, t);
			t &= MultiID.ModPow(t, up, z);
			wki = 1;
			ok = true;
			eti = 0;
			if (t == 1) return;
			MultiID.Add1SI(z, -1, temp);
			if (t == temp)
			{
				eti =1;
				return;
			}
			do
			{
				wki++;
				MultiID.SetToSI(t, w);
				MultiID.PowMBSI(t, MultiID.sTwo, z, t);
				MultiID.Add1SI(z, -1, temp);
				if (t == temp)
				{
					eti = 1;
					return;
				}
			}
			while (t != 1 && !MultiID.muErr && !MultiID.muAbort); // prevent looping
			MultiID.Add1SI(w, -1, temp);
			MultiID.GCDSI(temp, z, w);
			ok = (w == 1);
			MultiID.DoneSI(temp);
		} // Test5_2

		private static void Test3_P( // Subroutines for P>2, Check condition 3 for P>2
			out int ei,
			int nqi,
			out bool ok,
			int pi,
			int[] qa,
			ref int qi,
			MultiSI z
			)
		{
			if (MultiID.debug)
			{
				int tot = 0;
				string st = " In Test3_P: " + "ei=" + ".." + " nqi=" + nqi + " ok=" + ".." +
					" pi=" + pi + " qi=" + qi + "qa[]=";
				MultiID.WriteT(st, ref tot);
				MultiID.WriteT(" z=", ref tot);
				MultiID.WriteSILn(z, ref tot, 0, false);
			}
			int ji;
			int rem;
			ei = 0; // to guarantee output item
			ok = false;
			for (ji = 0; ji <= nqi; ji++)
			{
				qi = qa[ji];
				ei = Math.DivRem(qi - 1, pi, out rem);
				if (rem == 0)
				{
					if (MultiID.ModPow(z, SI(ei), SI(qi)) != MultiID.sOne)
					{
						ok = true;
						break;
					}
				}
			}
			if (MultiID.debug) // L2450:
			{
				int tot = 0;
				string st = "Out Test3_P: " + "ei=" + ei + " nqi=" + nqi + " ok=" + ok +
					" pi=" + pi + " qi=" + qi + "qa[]=";
				MultiID.WriteT(st, ref tot);
				MultiID.WriteT(" z=", ref tot);
				MultiID.WriteSILn(z, ref tot, 0, false);
			}
		} // Test3_P

		private static void Test4_P( // Check condition 4 for P>2, decompose N^(P%-1)-1 to P%^H%*U
			out int hi,
			int pi,
			MultiSI up,
			MultiSI wp,
			MultiSI z
			)
		{
			MultiID.useFMB = false;
			wp = (z ^ (pi - 1)) - 1;
			hi = -1;
			do
			{
				up &= wp;
				wp /= pi;
				hi++;
			}
			while (MultiID.rem1 == 0 && !MultiID.muErr && !MultiID.muAbort);
		} // Test4_P

		private static void Test2_P( // Check condition 4 for P>2, Calculation of tau(chi)^n/tau(chi^n)
			MultiSI e,
			out int ei,
			int gi,
			int[] indexa,
			out bool ok,
			int[] phia,
			int[] psia,
			int pi,
			int qi,
			MultiSI[] sum,
			bool swb,
			MultiSI[] tau,
			MultiSI[] tau1,
			MultiSI[] tau2,
			MultiSI[] taue,
			MultiSI[] taup,
			MultiSI w,
			out int wi,
			MultiSI z
			)
		{
			int ii;
			wi = 1;
			for (ii = 1; ii < qi; ii++) // calc index
			{
				wi = (wi * gi) % qi;
				indexa[wi] = ii % pi;
			}
			TAU_P(indexa, phia, pi, psia, qi, sum, tau, tau1, tau2, w, ref wi, z);
			if (!swb)
				for (ib = 0; ib <= pi - 2; ib++) MultiID.SetToSI(tau[ib], taup[ib]);
			MultiID.Div1SI(z, pi, e);
			T_E(e, pi, tau, tau1, tau2, w, ref wi, z);
			for (ib = 0; ib <= pi - 2; ib++) MultiID.SetToSI(tau[ib], taue[ib]);
			ei = z % pi;
			JACOBI(ei, indexa, phia, pi, psia, qi, sum, tau, tau1, tau2, w, ref wi, z);
			T_TM(pi, tau, tau2, taue, w, ref wi, z);
			// MAIN OF TEST2
			ei = (indexa[z % qi] * (-z)) % pi;
			if (ei == pi - 1)
				w &= -(z - 1) * (pi - 1);
			else if (tau[ei] != 1)
			{
				ok = false;
				return;
			}
			else
			{
				tau[ei]--;
				w &= 0;
			}
			for (ii = 0; ii <= pi - 2; ii++)
				w += tau[ii];
			ok = (w == 0);
		} // Test2_P

		private static void Test5_P( // Check condition 4 for P>2, calculation of tan(chi)^p^u
			ref int eti,
			out bool ok,
			int pi,
			MultiSI t,
			MultiSI[] tau,
			MultiSI[] tau1,
			MultiSI[] tau2,
			MultiSI[] taue,
			MultiSI[] taup,
			MultiSI up,
			MultiSI w,
			ref int wi,
			out int wki,
			MultiSI z
			)
		{
			int tot = 0;
			int ii, ji;
			for (ib = 0; ib <= pi - 2; ib++) MultiID.SetToSI(taup[ib], tau1[ib]);
			MultiID.SetTo1SI(1, tau[0]);
			for (ib = 1; ib <= pi - 2; ib++) MultiID.ClearSI(tau[ib]);
			int upLenM = Len(up) - 2;
			for (ii = 0; ii <= upLenM; ii++)
			{
				if (Bit(ii) != 0)
					T_T1(pi, tau, tau1, tau2, w, ref wi, z);
				T1_T1(pi, tau1, tau2, w, ref wi, z);
			}
			T_T1(pi, tau, tau1, tau2, w, ref wi, z);
			wki = 1;
			W_chi(ref eti, out ok, pi, t, tau, w, z);
			if (ok) return;
			do
			{
				for (ib = 0; ib <= pi - 2; ib++) MultiID.SetToSI(tau[ib], tau1[ib]);
				for (ib = 0; ib <= pi - 2; ib++) MultiID.SetToSI(tau[ib], taue[ib]);
				int piLenM = Len(pi) - 2;
				for (ii = 1; ii <= piLenM; ii++)
				{
	  				T1_T1(pi, tau1, tau2, w, ref wi, z);
					if (Bit(ii, pi) != 0)
						T_T1(pi, tau, tau1, tau2, w, ref wi, z);
				}
				T1_T1(pi, tau1, tau2, w, ref wi, z);
				T_T1(pi, tau, tau1, tau2, w, ref wi, z);
				wki++;
				W_chi(ref eti, out ok, pi, t, tau, w, z);
				MultiID.MuInterrupt();
				if (MultiID.muAbort) return; // exit function FnAdleman
				if (MultiID.statReq2)
				{
					MultiID.WriteT(" FnAdleman: Test5_P", ref tot);
					MultiID.TimeDel(ref tot);
					MultiID.statReq2 = false;
				}
			}
			while (!ok);
			if (eti > 0)
			{
				ok = true;
				return;
			}
			ii = 0;		// MAIN OF TEST 5 // never executed
			while (true)
			{
				t &= taue[ii];
				if (GCD(t, z) == 1) break;
				if (t != 0)
				{
					w &= GCD(t, z);
					ok = false;
					return;
				}
				ii++;
			}
			do // do once // L3140:
			{
				if (t != 1)
				{
					if (GCD(t - 1, z) == 1) break;
					w &= GCD(t - 1, z);
					ok = false;
					return;
				}
				else
				{
					ji = ii + 1;
					while (true)
					{
						t &= taue[ji];
						if (GCD(t, z) == 1) break;
						if (t != 0)
						{
							w &= GCD(t, z);
							ok = false;
							return;
						}
						ji++;
					}
				}
			} while (false);
			ii = 0; // L3240:
			while (true)
			{
				t &= taue[ii];
				if (GCD(t + 1, z) == 1)
				{
					ok = true;
					return;
				}
				if (t < (z - 1))
				{
					w &= GCD(t + 1, z);
					ok = false;
					return;
				}
				ii++;
			}
		} // Test5_P

		private static void W_chi( // decide w(chi)
			ref int eti,
			out bool ok,
			int pi,
			MultiSI t,
			MultiSI[] tau,
			MultiSI w,
			MultiSI z
			)
		{
			MultiSI temp = new MultiSI(z.N);
			int ii;
			MultiID.SetToSI(tau[0], t);
			if (t > MultiID.sOne)
			{
				MultiID.SetToSI(z, temp);
				MultiID.RSubSI(MultiID.sOne, temp);		// temp is z - 1
			}
			if (t <= MultiID.sOne)
			{
				MultiID.SetToSI(t, w);
				eti = 0;
				for (ii = 1; ii <= pi - 2; ii++)
				{
					MultiID.SetToSI(tau[ii], t);
					if (t != 0)
					{
						MultiID.RAddSI(t, w);
						eti = ii;
					}
				}
				ok = (w == 1);
			}
			else if (t < temp) ok = false;				// temp is z - 1
			else
			{
				for (ii = 1; ii <= pi - 2; ii++)
				{
					if (tau[ii] < t)
					{
						ok = false;
						return;
					}
				}
				ok = true;
				eti = pi - 1;
			}
			MultiID.DoneSI(temp);
		} // W_chi

		private static void TAU_P( // calculation of tau(chi)^p
			int[] indexa,
			int[] phia,
			int pi,
			int[] psia,
			int qi,
			MultiSI[] sum,
			MultiSI[] tau,
			MultiSI[] tau1,
			MultiSI[] tau2,
			MultiSI w,
			ref int wi,
			MultiSI z
			)
		{
			int ii, ji;
			for (ib = 1; ib < qi; ib++) psia[ib] = indexa[ib];
			for (ib = 1; ib < qi; ib++) phia[ib] = indexa[ib];
			MultiID.SetTo1SI(1, tau[0]);
			for (ib = 1; ib <= pi - 2; ib++) MultiID.ClearSI(tau[ib]); // set 1
			MultiID.SetTo1SI(1, tau1[0]);
			for (ib = 1; ib <= pi - 2; ib++) MultiID.ClearSI(tau1[ib]); // set 1
			int piLenM = Len(pi) - 2;
			for (ii = 1; ii <= piLenM; ii++)
			{
				PSI_PSI(pi, psia, qi, sum, tau1, tau2, w, ref wi, z);
				if (Bit(ii, pi) != 0)
				{
					PHY_PSI(phia, pi, psia, qi, sum, tau, tau1, tau2, w, ref wi, z);					}
			}
			PSI_PSI(pi, psia, qi, sum, tau1, tau2, w, ref wi, z);
			T_T1(pi, tau, tau1, tau2, w, ref wi, z);
			for (ib = 0; ib <= pi - 2; ib++) MultiID.SetToSI(tau[ib], tau1[ib]);
			MultiID.ClearSI(tau1[pi - 1]);
			ji = phia[qi - 1];
			for (ii = 0; ii < pi; ii++)
			{
				wi = (ii + ji) % pi;
				tau[wi] &= tau1[ii];
			}
			w &= tau[pi - 1];
			for (ii = 0; ii <= pi - 2; ii++)
				tau[ii] &= (tau[ii] - w) * qi % z;
		} // TAU_P

		private static void T_E( // calculation of tau(chi)^p
			MultiSI e,
			int pi,
			MultiSI[] tau,
			MultiSI[] tau1,
			MultiSI[] tau2,
			MultiSI w,
			ref int wi,
			MultiSI z
			)
		{
					int ii, ji;
			for (ib = 0; ib <= pi - 2; ib++) MultiID.SetToSI(tau[ib], tau1[ib]);
			ji = 0;
			int eLen = Len(e);
			if (Even(e))
			{
				while (Bit(ji) == 0 && !MultiID.muErr && !MultiID.muAbort) // prevent looping)
				{
					if (MultiID.muErr) return;
					T1_T1(pi, tau1, tau2, w, ref wi, z);
					ji++;
				}
				for (ib = 0; ib <= pi - 2; ib++) MultiID.SetToSI(tau1[ib], tau[ib]);
			}
			for (ii = ji + 1; ii < eLen; ii++)
			{
				T1_T1(pi, tau1, tau2, w, ref wi, z);
				if (Bit(ii) == 1)
					T_T1(pi, tau, tau1, tau2, w, ref wi, z);
			}
		}

		private static void JACOBI( // calculation of  tau(chi)^e/tau(chi^e)
			int ei,
			int[] indexa,
			int[] phia,
			int pi,
			int[] psia,
			int qi,
			MultiSI[] sum,
			MultiSI[] tau,
			MultiSI[] tau1,
			MultiSI[] tau2,
			MultiSI w,
			ref int wi,
			MultiSI z
			)
		{
			int ii, ji;
			for (ib = 0; ib < qi; ib++) psia[ib] = indexa[ib];
			MultiID.SetTo1SI(1, tau1[0]);
			for (ib = 1; ib <= pi - 2; ib++) MultiID.ClearSI(tau1[ib]);
			ji = 0;
			while (Bit(ji, ei) == 0 && !MultiID.muErr && !MultiID.muAbort) // prevent looping
			{
				if (MultiID.muErr) return;
				PSI_PSI(pi, psia, qi, sum, tau1, tau2, w, ref wi, z);
				ji++;
			}
			for (ib = 0; ib <= pi - 2; ib++) MultiID.SetToSI(tau1[ib], tau[ib]);
			for (ib = 0; ib < qi; ib++) phia[ib] = psia[ib];
			int eiLen = Len(ei);
			for (ii = ji + 1; ii < eiLen; ii++)
			{
				PSI_PSI(pi, psia, qi, sum, tau1, tau2, w, ref wi, z);
				if (Bit(ii, ei) == 1)
					PHY_PSI(phia, pi, psia, qi, sum, tau, tau1, tau2, w, ref wi, z);
			}
		} // JACOBI

		private static void T_TM( // calculation of  T=T*TM
			int pi,
			MultiSI[] tau,
			MultiSI[] tau2,
			MultiSI[] taue,
			MultiSI w,
			ref int wi,
			MultiSI z
			)
		{
			int ii, ji;
			for (
				ji = 0; ji <= pi - 2; ji++)
			{
				tau2[ji] &= tau[0] * taue[ji] % z;
			}
			tau2[pi - 1] &= 0;
			for (ii = 1; ii <= pi - 2; ii++)
			{
				for (ji = 0; ji <= pi - 2; ji++)
				{
					wi = (ii + ji) % pi;
					tau2[wi] &= tau2[wi] + tau[ii] * taue[ji] % z;
				}
			}
			w &= tau2[pi - 1];
			for (ii = 0; ii <= pi - 2; ii++)
				tau[ii] &= (tau2[ii] - w) % z;
			for (ib = 0; ib <= pi - 2; ib++) MultiID.SetToSI(tau[ib], tau2[ib]);
		} // T_TM

		private static void PSI_PSI( // calculation of  the sum of PSI(A)*PSI(1-A)
			int pi,
			int[] psia,
			int qi,
			MultiSI[] sum,
			MultiSI[] tau1,
			MultiSI[] tau2,
			MultiSI w,
			ref int wi,
			MultiSI z
			)
		{
			int pWords = 2 * z.N + 2;
			MultiSI temp = new MultiSI(pWords);
			MultiSI tempR = new MultiSI(pWords);
			int ii, ji;
			for (ib = 0; ib < pi; ib++) MultiID.ClearSI(sum[ib]);
			for (ii = 2; ii < qi; ii++)
			{
				MultiID.RAdd1SI( 1, sum[ (psia[ii] + psia[qi + 1 - ii]) % pi ] );
			}
			// calculation of  T1=T1*T1*U
			for (ii = 0; ii <= (pi - 1) / 2; ii++)
			{
				MultiID.SqSI(tau1[ii], temp);
				MultiID.DivQRSI(temp, z, tau2[2 * ii], temp);
			}
			for (ii = (pi + 1) / 2; ii <= pi - 2; ii++)
			{
				MultiID.SqSI(tau1[ii], temp);
				MultiID.DivQRSI(temp, z, tau2[2 * ii - pi], temp);
			}
			MultiID.ClearSI(tau2[pi - 2]);
			for (ii = 0; ii <= pi - 3; ii++)
			{
				for (ji = ii + 1; ji <= pi - 2; ji++)
				{
					wi = (ii + ji) % pi;
					MultiID.MulSI(tau1[ii], tau1[ji], temp);
					MultiID.RMul1SI(2, temp);
					MultiID.DivQRSI(temp, z, tempR, temp);
					MultiID.RAddSI(tempR, tau2[wi]);
				}
			}
			for (ji = 0; ji < pi; ji++)
			{
				MultiID.MulSI(tau2[0], sum[ji], temp);
				MultiID.DivQRSI(temp, z, tau1[ji], temp);
			}
			for (ii = 1; ii < pi; ii++)
			{
				for (ji = 0; ji < pi; ji++)
				{
					wi = (ii + ji) % pi;
					MultiID.MulSI(tau2[ii], sum[ji], temp);
					MultiID.DivQRSI(temp, z, tempR, temp);
					MultiID.RAddSI(tempR, tau1[wi]);
				}
			}
			MultiID.DoneSI(tempR);
			MultiID.SetToSI(tau1[pi - 1], w);
			for (ii = 0; ii <= pi - 2; ii++)
			{
				MultiID.SubSI(tau1[ii], w, temp);
				MultiID.DivQRSI(temp, z, tau1[ii], temp);
				if (tau1[ii].S) tau1[ii] += z;
			}
			for (ii = 0; ii < qi; ii++)
			{
				psia[ii] = 2 * psia[ii] % pi;
			}
			MultiID.DoneSI(temp);
		} // PSI_PSI

		private static void PHY_PSI( // calculation of  the sum of PHY(A)*PSI(1-A)
			int[] phia,
			int pi,
			int[] psia,
			int qi,
			MultiSI[] sum,
			MultiSI[] tau,
			MultiSI[] tau1,
			MultiSI[] tau2,
			MultiSI w,
			ref int wi,
			MultiSI z
			)
		{
			int pWords = 2 * z.N + 2;
			MultiSI temp = new MultiSI(pWords);
			MultiSI tempR = new MultiSI(pWords);
			int ii, ji;
			for (ib = 0; ib < pi; ib++) MultiID.ClearSI(sum[ib]);
			for (ii = 2; ii < qi; ii++)
			{
				MultiID.RAdd1SI(1, sum[ (phia[ii] + psia[qi + 1 - ii]) % pi] );
			}
			for (ji = 0; ji <= pi - 2; ji++)
			{
				MultiID.MulSI(tau[0], tau1[ji], temp);
				MultiID.DivQRSI(temp, z, tau2[ji], temp);
			}
			MultiID.ClearSI(tau2[pi - 1]);
			for (ii = 1; ii <= pi - 2; ii++)
			{
				for (ji = 0; ji <= pi - 2; ji++)
				{
					wi = (ii + ji) % pi;
					MultiID.MulSI(tau[ii], tau1[ji], temp);
					MultiID.DivQRSI(temp, z, tempR, temp);
					MultiID.RAddSI(tempR, tau2[wi]);
				}
			}
			for (ji = 0; ji < pi; ji++)
			{
				MultiID.MulSI(tau2[0], sum[ji], temp);
				MultiID.DivQRSI(temp, z, tau[ji], temp);
			}
			for (ii = 1; ii < pi; ii++)
			{
				for (ji = 0; ji < pi; ji++)
				{
					wi = (ii + ji) % pi;
					MultiID.MulSI(tau2[ii], sum[ji], temp);
					MultiID.DivQRSI(temp, z, tempR, temp);
					MultiID.RAddSI(tempR, tau[wi]);
				}
			}
			MultiID.DoneSI(tempR);
			MultiID.SetToSI(tau[pi - 1], w);
			for (ii = 0; ii <= pi - 2; ii++)
			{
				MultiID.SubSI(tau[ii], w, temp);
				MultiID.DivQRSI(temp, z, tau[ii], temp);
				if (tau[ii].S) tau[ii] += z;
			}
			MultiID.DoneSI(temp);
			// calculation of  PHY=PHY*PSY
			for (ii = 0; ii < qi; ii++)
			{
				phia[ii] = (phia[ii] + psia[ii]) % pi;
			}
		} // PHY_PSI

		private static void T_T1( // calculation of  T=T*T1
			int pi,
			MultiSI[] tau,
			MultiSI[] tau1,
			MultiSI[] tau2,
			MultiSI w,
			ref int wi,
			MultiSI z
			)
		{
			int pWords = 2 * z.N + 2;
			MultiSI temp = new MultiSI(pWords);
			MultiSI tempR = new MultiSI(pWords);
			if (MultiID.debug)
			{
				int tot = 0;
				string st = " In T_T1: " + "pi=" + pi + " wi=" + wi;
				MultiID.WriteT(st, ref tot);
				MultiID.WriteT(" w=", ref tot);
				MultiID.WriteSI(w, ref tot, 0, false);
				MultiID.WriteT(" z=", ref tot);
				MultiID.WriteSILn(z, ref tot, 0, false);
			}
			int ii, ji;
			for (ji = 0; ji <= pi - 2; ji++)
			{
				MultiID.MulSI(tau[0], tau1[ji], temp);
				MultiID.DivQRSI(temp, z, tau2[ji], temp);
			}
			MultiID.ClearSI(tau2[pi - 1]);
			for (ii = 1; ii <= pi - 2; ii++)
			{
				for (ji = 0; ji <= pi - 2; ji++)
				{
					wi = (ii + ji) % pi;
					MultiID.MulSI(tau[ii], tau1[ji], temp);
					MultiID.DivQRSI(temp, z, tempR, temp);
					MultiID.RAddSI(tempR, tau2[wi]);
				}
			}
			MultiID.DoneSI(tempR);
			MultiID.SetToSI(tau2[pi - 1], w);
			for (ii = 0; ii <= pi - 2; ii++)
			{
				MultiID.SubSI(tau2[ii], w, temp);
				MultiID.DivQRSI(temp, z, tau[ii], temp);
				if (tau[ii].S) tau[ii] += z;
			}
			MultiID.DoneSI(temp);
			if (MultiID.debug)
			{
				int tot = 0;
				string st = "Out T_T1: " + "pi=" + pi + " wi=" + wi;
				MultiID.WriteT(st, ref tot);
				MultiID.WriteT(" w=", ref tot);
				MultiID.WriteSI(w, ref tot, 0, false);
				MultiID.WriteT(" z=", ref tot);
				MultiID.WriteSILn(z, ref tot, 0, false);
			}
		} // T_T1

		private static void T1_T1(		// calculation of  T1=T1*T1
			int pi,
			MultiSI[] tau1,
			MultiSI[] tau2,
			MultiSI w,
			ref int wi,
			MultiSI z
			)
		{
			int pWords = 2 * z.N + 2;
			MultiSI temp = new MultiSI(pWords);
			MultiSI tempR = new MultiSI(pWords);
			if (MultiID.debug)
			{
				int tot = 0;
				string st = " In T1_T1: " + "pi=" + pi + " wi=" + wi;
				MultiID.WriteT(st, ref tot);
				MultiID.WriteT(" w=", ref tot);
				MultiID.WriteSI(w, ref tot, 0, false);
				MultiID.WriteT(" z=", ref tot);
				MultiID.WriteSILn(z, ref tot, 0, false);
			}
			int ii, ji;
			for (ii = 0; ii <= (pi - 1) / 2; ii++)
			{
				MultiID.SqSI(tau1[ii], temp);
				MultiID.DivQRSI(temp, z, tau2[2 * ii], temp);
			}
			for (ii = (pi + 1) / 2; ii <= pi - 2; ii++)
			{
				MultiID.SqSI(tau1[ii], temp);
				MultiID.DivQRSI(temp, z, tau2[2 * ii - pi], temp);
			}
			MultiID.ClearSI(tau2[pi - 2]);
			for (ii = 0; ii <= pi - 3; ii++)
			{
				for (ji = ii + 1; ji <= pi - 2; ji++)
				{
					wi = (ii + ji) % pi;
					MultiID.MulSI(tau1[ii], tau1[ji], temp);
					MultiID.RMul1SI(2, temp);
					MultiID.DivQRSI(temp, z, tempR, temp);
					MultiID.RAddSI(tempR, tau2[wi]);
				}
			}
			MultiID.SetToSI(tau2[pi - 1], w);
			for (ii = 0; ii <= pi - 2; ii++)
			{
				MultiID.SubSI(tau2[ii], w, temp);
				MultiID.DivQRSI(temp, z, tau1[ii], temp);
				if (tau1[ii].S) tau1[ii] += z;
			}
			if (MultiID.debug)
			{
				int tot = 0;
				string st = "Out T1_T1: " + "pi=" + pi + " wi=" + wi;
				MultiID.WriteT(st, ref tot);
				MultiID.WriteT(" w=", ref tot);
				MultiID.WriteSI(w, ref tot, 0, false);
				MultiID.WriteT(" z=", ref tot);
				MultiID.WriteSILn(z, ref tot, 0, false);
			}
		} // T1_T1

		private static void GET_V( // subroutine for the final test, calculation of  V mod S
			MultiSI iN,
			int[,] eta,
			MultiSI m1,
			int npi,
			int nqi,
			MultiSI p1,
			int[] ga,
			int[] pa,
			out int pi,
			int[] qa,
			ref int qi,
			MultiSI u1,
			MultiSI v,
			MultiSI v1,
			int[,] wka,
			out int xi
			)
		{
			xi = pi = 0; // to guarantee output items
			int ii, ji;
			for(ji = 0; ji <= nqi; ji++)
			{
				iN &= eta[0, ji];
				p1 &= 2;
				for (ii = 1; ii <= npi; ii++)
				{
					if (wka[ii, ji] != 0)
					{
						pi = pa[ii];
						u1 &= p1 * ModInv(p1, SI(pi));
						v1 &= 1 - u1;
						p1 *= pi;
						iN &= (iN * v1 + eta[ii, ji] * u1) % p1;
					}
				}
				qi = qa[ji];
				xi = Int(MultiID.ModPow(SI(ga[ji]), iN, SI(qi)));
				if (ji != 0)
				{
					u1 &= m1 * ModInv(m1, SI(qi));
					v1 &= 1 - u1;
					m1 *= qi;
					v &= (v * v1 + xi * u1) % m1;
				}
				else
				{
					m1 &= qi;
					v &= xi;
				}
			}
		} // GET_V

		private static void FnPowerCheck(MultiSI z, MultiSI n1) // POWERCHK V1.1
		// POWERCHK V1.1, 1987,89 by Yuji KIDA
		{
			int tot = 0;
			int ei;
			int pWords = 2 * z.N + 2;
			MultiSI dp = new MultiSI(pWords);
			MultiSI e = new MultiSI(pWords);				// ei in SI format
			MultiSI wp = new MultiSI(pWords);
			MultiSI epp = new MultiSI(pWords);
			MultiSI temp = new MultiSI(pWords);
			MultiID.useFMB = false;
			MultiID.DiagT("Power Check ", ref tot);
			for (ei = (Len(z) - 1) / 17; ei > 1; ei--)
			{
				MultiID.DiagT(ei + " ", ref tot);
				double zD;	// z in double format
				MultiID.GetDbSI(out zD, z);
				double dD = Math.Exp(Math.Log(zD) / ei);	// dp in double format
				MultiID.SetToDbSI(dD, dp);
				MultiID.SetTo1SI(ei, e);
				MultiID.PowSI(dp, e, wp);
				MultiID.Add1SI(dp, 1, temp);
				MultiID.PowSI(temp, e, temp);
				MultiID.SubSI(temp, wp, epp);
				MultiID.useFMB = false;
				do
				{
					if (MultiID.muErr) return;
					MultiID.Mul1SI(wp, ei - 1, temp);		// main of Newton’s method
					MultiID.RAddSI(z, temp);
					MultiID.RMulSI(temp, dp);
					MultiID.Mul1SI(wp, ei, temp);
					MultiID.RDivSI(temp, dp);
					wp &= dp ^ ei;
					MultiID.SubSI(wp, z, temp);
					temp.S = false;
				} while (temp >= epp);
				MultiID.PowSI(dp, e, wp);		// test for each integer
				if (wp <= z)
				{
					if (wp == z) goto ExitPowerCheck;
					do	// if smaller
					{
						if (MultiID.muErr || MultiID.muAbort) return; // L5450:
						MultiID.RAdd1SI(1, dp);
						MultiID.PowSI(dp, e, wp);
						if (wp == z) goto ExitPowerCheck;
					} while (wp < z);
					continue; // next ei
				}
				do // L5500:
				{
					MultiID.RAdd1SI(-1, dp);		// never executed
					MultiID.PowSI(dp, e, wp);
					if	(wp == z) goto ExitPowerCheck;
				} while (wp > z);
			}
			MultiID.DoneSI(temp);
			MultiID.DoneSI(epp);
			MultiID.DoneSI(wp);
			MultiID.DoneSI(e);
			MultiID.SetTo1SI(1, dp);		// not a power
		ExitPowerCheck: // L5550:
			if (MultiID.diagOn) Common.SingleSpace();
			MultiID.SetToSI(dp, n1);
			MultiID.DoneSI(dp);
			return;
		} // fnPowerCheck

		private static void FnECM(MultiSI z, MultiSI n3) // ELLIPTIC CURVE METHOD ver 4.2
		// originally   by H. W. LENSTRA
		// fast version by P. L. Montgomery
		// fast version by H. Suyama
		// implemented  by Y. Kida
		{
			int i;
			int tot = 0;
			MultiID.useFMB = false;
			int eci, ii, ipi, ji;
			int pWords = 2 * z.N + 2;
			MultiSI gd = new MultiSI(pWords);
  			MultiSI a  = new MultiSI(pWords);
			MultiSI a0 = new MultiSI(pWords);
			MultiSI a0s = new MultiSI(2 * pWords); // a0^2
			MultiSI aa = new MultiSI(pWords);
			MultiSI dx = new MultiSI(pWords);
			MultiSI dz = new MultiSI(pWords);
			int l1i;
			int l2i;
			double lgd;
			int lsi;
			MultiSI n1 = new MultiSI(pWords);
			MultiSI n2 = new MultiSI(pWords);
			int pi;
			int qi;
  			MultiSI tx = new MultiSI(pWords);
			MultiSI tz = new MultiSI(pWords);
			MultiSI ux = new MultiSI(pWords);
			MultiSI uz = new MultiSI(pWords);
			MultiSI w1 = new MultiSI(pWords);
			MultiSI w2 = new MultiSI(pWords);
			MultiSI w3 = new MultiSI(pWords);
			MultiSI w4 = new MultiSI(pWords);
			MultiSI wx = new MultiSI(pWords);
			MultiSI wz = new MultiSI(pWords);
			MultiSI x  = new MultiSI(pWords);
			MultiSI zz  = new MultiSI(pWords);
			MultiSI mp = new MultiSI(1);
			MultiSI[] xa = new MultiSI[48];
			for (i = 0; i < 48; i++) xa[i] = new MultiSI();
			MultiSI[] v0a = new MultiSI[12];
			MultiSI[] v1a = new MultiSI[12];
			MultiSI[] v2a = new MultiSI[12];
			MultiSI[] v3a = new MultiSI[12];
			for (i = 0; i < 12; i++)
			{
				v0a[i] = new MultiSI();
				v1a[i] = new MultiSI();
				v2a[i] = new MultiSI();
				v3a[i] = new MultiSI();
			}
			if (MultiID.diagOn) MultiID.WriteT("EC Method with limit ", ref tot);
			double zd;
			MultiID.GetDbSI(out zd, z);
			l1i = (int)Math.Floor(Math.Pow(Math.Log(zd), 2.65)/10);
			l1i = Math.Min(130000, Math.Max(l1i, 500));
			l2i = 40 * l1i;
			if (MultiID.diagOn) MultiID.WriteT("(B1 = " + l1i + " B2 = " + l2i + ") ", ref tot);
			lgd = Math.Log(l1i);
			lsi = (int)Math.Floor(Math.Sqrt(l1i));
			eci = 0;
			if (MultiID.diagOn) MultiID.WriteT("Curve ", ref tot);
			do // SetCurve:
			{
				do
				{
					do
					{
						MultiID.MuInterrupt();
						eci++;
						if (MultiID.diagOn) MultiID.WriteT(eci + " ", ref tot);
						a0 &= 2 * (eci + 1) * ModInv(SI(3 * (eci + 1) * (eci + 1) - 1), z) % z;
						if (MultiID.muAbort) return;
						MultiID.useFMB = false;
						a0s &= a0 * a0;
					} while (a0 * (a0s - 1) * (9 * a0s - 1) % z == 0);
					a &= (-3 * (a0s * a0s) - 6 * a0s + 1) * ModInv(4 * (a0s * a0), z) % z;
					aa &= (a + 2) * ModInv(SI(4), z) % z;
					x &= (3 * a0s + 1) % z;
					zz &= 4 * a0 % z;
					int l1iLen = Len(l1i);
					for (ii = 1; ii <= l1iLen; ii++)// 1st step, for power of 2
					{
						w1 &= ((x + zz)^2) % z;
						w2 &= ((x - zz)^2) % z;
						x &= w1 * w2 % z;
						zz &= (w1 - w2) * (w2 + aa * (w1 - w2) % z) % z;
					}
					MultiID.GCDSI(zz, z, gd);
					if (gd > MultiID.sOne) break;
					ipi = 2;		// for powers of odd primes
					do
					{
						pi = pTab[ipi];
						mp = SI(Math.Pow(pi, Math.Floor(lgd / Math.Log(pi)))); // can mp overflow a double
						ipi++;
						ECM_SUB(aa, mp, tx, tz, ux, uz, w1, w2, w3, w4, x, zz, z);
						MultiID.GCDSI(zz, z, gd);
						if (gd > MultiID.sOne) break;
					} while (pi <= lsi);
					if (gd > MultiID.sOne) break;
					do
					{
						mp &= 1;
						for (i = 0; i <= 9; i++, ipi++)
						{
							mp *= pTab[ipi];
						}
						ECM_SUB(aa, mp, tx, tz, ux, uz, w1, w2, w3, w4, x, zz, z);
						MultiID.GCDSI(zz, z, gd);
						if (gd > MultiID.sOne) break;
					} while (pTab[ipi + 9] <= l1i);						// * 2nd STEP
  					if (gd > MultiID.sOne) break;
					ux &= x;											// =Q
					uz &= zz;
					xa[0] &= x * ModInv(zz, z) % z;
					ji = 0;
					w1 &= ((x + zz)^2) % z;
					w2 &= ((x - zz)^2) % z;
					tx &= w1 * w2 % z;									// =2Q
					tz &= (w1 - w2) * (w2 + aa * (w1 - w2) % z) % z;
					w1 &= (x - zz) * (tx + tz) % z;
					w2 &= (x + zz) * (tx - tz) % z;
					x &= ((w1 + w2)^2) % z * uz % z;					// =3Q
					zz &= ((w1 - w2)^2) % z * ux % z;
					for (ii = 5; ii <= 209; ii += 2)
					{
						wx &= x;
						wz &= zz;
						w1 &= (x - zz) * (tx + tz) % z;
						w2 &= (x + zz) * (tx - tz) % z;
						x &= ((w1 + w2)^2) % z * uz % z;				// =7Q,9Q,...
						zz &= ((w1 - w2)^2) % z * ux % z;
						if (ii == 105)									// =105Q
						{
							dx &= x;
							dz &= zz;
						}
						if (GCD(ii, 210) == 1)
						{
							ji++;
							xa[ji] &= x * ModInv(zz, z) % z;
						}
						ux &= wx;
						uz &= wz;
					}
					w1 &= ((dx + dz)^2) % z;
					w2 &= ((dx - dz)^2) % z;
					x &= w1 * w2 % z;								// =210Q
					zz &= (w1 - w2) * (w2 + aa * (w1 - w2) % z) % z;
					for (ji = 0; ji <= 11; ji++)
					{
						ii = ji * 4;
						w1 &= xa[ii] + xa[ii + 1];
						w2 &= xa[ii] * xa[ii + 1] % z;
						w3 &= xa[ii + 2] + xa[ii + 3];
						w4 &= xa[ii + 2] * xa[ii + 3] % z;
						v3a[ji] &= w1 + w3;
						v2a[ji] &= w1 * w3 % z + w2 + w4;
						v1a[ji] &= (w2 * w3 + w1 * w4) % z;
						v0a[ji] &= w2 * w4 % z;
					}
					ux &= x;										// =210Q
					uz &= zz;
					w1 &= ((x + zz)^2) % z;
					w2 &= ((x - zz)^2) % z;
					tx &= w1 * w2 % z;								// =420Q
					tz &= (w1 - w2) * (w2 + aa * (w1 - w2) % z) % z;
					w1 &= (x - zz) * (tx + tz) % z;
					w2 &= (x + zz) * (tx - tz) % z;
					x &= ((w1 + w2)^2) % z * uz % z;				// =630Q
					zz &= ((w1 - w2)^2) % z * ux % z;
					for (qi = 2; qi <= l1i / 420; qi++)
					{
						wx &= x;
						wz &= zz;
						w1 &= (x - zz) * (tx + tz) % z;
						w2 &= (x + zz) * (tx - tz) % z;
						x &= ((w1 + w2)^2) % z * uz % z;
						zz &= ((w1 - w2)^2) % z * ux % z;
						ux &= wx;
						uz &= wz;
					}
					for ( ; qi <= l2i / 420; qi++)
					{
						w1 &= x * ModInv(zz, z) % z;
						w2 &= w1 * w1 % z;
						w3 &= w1 * w2 % z;
						w4 &= w2 * w2 % z;
						mp &= (w4 - v3a[0] * w3 + v2a[0] * w2 - v1a[0] * w1 + v0a[0]) % z;
						for (ji = 1; ji <= 11; ji++)
						{
							mp &= ((w4 - v3a[ji] * w3 + v2a[ji] * w2 - v1a[ji] * w1 + v0a[ji]) % z) * mp % z;
						}
						MultiID.GCDSI(mp, z, gd);
						if (gd > MultiID.sOne) break;
						wx &= x;
						wz &= zz;
						w1 &= (x - zz) * (tx + tz) % z;
						w2 &= (x + zz) * (tx - tz) % z;
						x &= ((w1 + w2)^2) % z * uz % z;
						zz &= ((w1 - w2)^2) % z * ux % z;
						ux &= wx;
						uz &= wz;
					}
					if (gd > MultiID.sOne) break;
					qi = 0;
				} while (true); // change curve
			} while (Len(z) < 150 && gd == z); // change curve // FindIt:
			v3a = new MultiSI[0];
			v2a = new MultiSI[0];
			v1a = new MultiSI[0];
			v0a = new MultiSI[0];
			MultiID.DoneSI(mp);
			MultiID.DoneSI(zz);
			MultiID.DoneSI(x);
			MultiID.DoneSI(wz);
			MultiID.DoneSI(wx);
			MultiID.DoneSI(w4);
			MultiID.DoneSI(w3);
			MultiID.DoneSI(w2);
			MultiID.DoneSI(w1);
			MultiID.DoneSI(uz);
			MultiID.DoneSI(ux);
			MultiID.DoneSI(tz);
			MultiID.DoneSI(tx);
			MultiID.DoneSI(n2);
			MultiID.DoneSI(n1);
			MultiID.DoneSI(dz);
			MultiID.DoneSI(dx);
			MultiID.DoneSI(aa);
			MultiID.DoneSI(a0s);
			MultiID.DoneSI(a0);
			MultiID.DoneSI(a);
			if (MultiID.diagOn) MultiID.WriteLn();
			n3 &= gd;
			MultiID.DoneSI(gd);
		} // FnECM

		private static void ECM_SUB( // main subroutine
			MultiSI aa,
  			MultiSI mp,
			MultiSI tx,
			MultiSI tz,
			MultiSI ux,
			MultiSI uz,
			MultiSI w1,
			MultiSI w2,
			MultiSI w3,
			MultiSI w4,
			MultiSI x,
			MultiSI zz,
			MultiSI z
			)
		{
			int ii;
			tx &= x;	// &= is the = sign for MultiSI
			tz &= zz;
			ux &= x;
			uz &= zz;
			w3 &= tx + tz;
			w4 &= tx - tz;
			MultiID.useFMB = false;
			int mpLen = Len(mp);
			for (ii = 1; ii < mpLen; ii++)		// time consuming loop so did it the long way
			// this shows an attempt at speedup, but only gained about 8% and a lot harder to read
			// still about 11 times longer running time than UBASIC
			{
				MultiID.MuInterrupt();
				if (MultiID.muAbort || MultiID.muErr) break; // exit for
				MultiID.MulSI(w3, w3, w1);		// w1 &= w3 * w3 % z;
				MultiID.ModSI(w1, z, w1);
				MultiID.MulSI(w4, w4, w2);		// w2 &= w4 * w4 % z;
				MultiID.ModSI(w2, z, w2);
				MultiID.MulSI(w1, w2, tx);		// tx &= w1 * w2 % z;
				MultiID.ModSI(tx, z, tx);
				MultiID.SubSI(w1, w2, tz);		// tz &= w1 - w2;
				MultiID.MulSI(aa, tz, w3);		// tz &= tz * (w2 + aa * tz % z) % z;
				MultiID.RAddSI(w2, w3);
				MultiID.ModSI(w3, z, w3);
				MultiID.MulSI(tz, w3, tz);
				MultiID.ModSI(tz, z, tz);
				MultiID.AddSI(tx, tz, w3);		// w3 &= (tx + tz);
				MultiID.SubSI(tx, tz, w4);		// w4 &= (tx - tz);
				if (Bit(ii) != 0) // uses saveBits from Len(mp) above
				{
					MultiID.AddSI(x, zz, tz);	// w1 &= w4 * (x + zz) % z;
					MultiID.MulSI(w4, tz, w1);
					MultiID.ModSI(w1, z, w1);
					MultiID.SubSI(x, zz, tz);	// w1 &= w4 * (x + zz) % z;
					MultiID.MulSI(w3, tz, w2);
					MultiID.ModSI(w2, z, w2);
					MultiID.AddSI(w1, w2, tz);	// x &= ((w1 + w2)^2) % z * uz % z;
					MultiID.MulSI(tz, tz, tx);
					MultiID.ModSI(tx, z, tx);
					MultiID.MulSI(tx, uz, x);
					MultiID.ModSI(x, z, x);
					MultiID.SubSI(w1, w2, tz);
					MultiID.MulSI(tz, tz, tx);	// zz &= ((w1 - w2)^2) % z * ux % z;
					MultiID.ModSI(tx, z, tx);
					MultiID.MulSI(tx, ux, zz);
					MultiID.ModSI(zz, z, zz);
				}
				else
				{
					MultiID.AddSI(ux, uz, tz);	// w1 &= w4 * (ux + uz) % z;
					MultiID.MulSI(w4, tz, w1);
					MultiID.ModSI(w1, z, w1);
					MultiID.SubSI(ux, uz, tz);	// w2 &= w3 * (ux - uz) % z;
					MultiID.MulSI(w3, tz, w2);
					MultiID.ModSI(w2, z, w2);
					MultiID.AddSI(w1, w2, tz);	// ux &= ((w1 + w2)^2) % z * zz % z;
					MultiID.MulSI(tz, tz, tx);
					MultiID.ModSI(tx, z, tx);
					MultiID.MulSI(tx, zz, ux);
					MultiID.ModSI(ux, z, ux);
					MultiID.SubSI(w1, w2, tz);	// uz &= ((w1 - w2)^2) % z * x % z;
					MultiID.MulSI(tz, tz, tx);
					MultiID.ModSI(tx, z, tx);
					MultiID.MulSI(tx, x, uz);
					MultiID.ModSI(uz, z, uz);
				}
				if (MultiID.statReq2)
				{
					MultiID.nDone = ii;
					MultiID.nToGo = mpLen - ii - 1;
					MultiID.WritenDone("ECM:");
				}
			}
			savedBits = new int[0]; // delete savedBits
		} // ECM_SUB

		private static int KroM1(int b)	// KroM1(b) == Kronecker(-1, b)
		// got help from Benoit Cloitre abcloitre(AT)wanadoo.fr. He pointed me to
		// http://www.research.att.com/cgi-bin/access.cgi/as/njas/sequences/eisA.cgi?Anum=A014577
		// and suggested the sequences are the sane if you change their 0 to -1
		// See my sequence A097402
		{
			int f = 1;							// f = final factor
			if (b < 0) { b = -b; f = -1; }
			if (b < 3) return f;
			do
			{
				if ((b & 1) == 1)				// if b is odd
				{
					if ((b & 2) == 0) return f;	// if b mod 4 = 1
					return -f;					// if b mod 4 = 3
				}
				b /= 2;
			} while (true);
		} // KroM1 int

		private static int KroM1(MultiSI bIn)	// KroM1(b) == Kronecker(-1, b)
		// got help from Benoit Cloitre abcloitre(AT)wanadoo.fr. He pointed me to
		// http://www.research.att.com/cgi-bin/access.cgi/as/njas/sequences/eisA.cgi?Anum=A014577
		// and suggested the sequences are the sane if you change their 0 to -1
		// See my sequence A097402
		{
			MultiSI b = new MultiSI(bIn.N);
			b &= bIn;
			int f = 1;									// f = final factor
			if (MultiID.CompSI(MultiID.sZero, b) < 0)	// if (b < 0) { b = -b; f = -1; }
			{
				b.S = false;
				f = -1;
			}
			if (MultiID.CompSI(MultiID.sThree, b) < 0)	// if (b < 3) return f;
			{
				return f;
			}
			do
			{
				if ((b.V[0] & 1) == 1)					// if b is odd
				{
					if ((b.V[0] & 2) == 0) return f;	// if b mod 4 = 1
					return -f;							// if b mod 4 = 3
				}
				b /= 2;
			} while (true);
		} // KroM1 SI

		private static int Kronecker(int a, int b)
		//
		// cf. Cohen p.29
		// return Kronecker symbol (a/b)
		// which, if b is an odd prime, == legendre symbol (a/b)
		//
		{
			if (a == -1) return KroM1(b);
			int f = 1;		// final factor
			if (a < 0)
			{
				f = KroM1(b);
				a = -a;
			}
			if (b < 0) b = -b;
			int[] tab2 = {0, 1, 0, -1, 0, -1, 0, 1};
			// (-1)^((a^2-1)/8) == tab2[a & 7]
			if (b == 0)
			{
				if (a != 1) return 0;
				else		return f;
			}
			if (a == 0)
			{
				if (b != 1)	return 0;
				else		return f;
			}
			// if a and b are both even
			if ((a & 1) == 0 && (b & 1) == 0) return 0;
			int v = 0;
			while ((b & 1) == 0)
			{
				v++;
				b >>= 1;
			}
			int k;
			if ((v & 1) == 0) k = 1;
			else			  k = tab2[ a & 7 ];
			//step3:
			do
			{
				if (a == 0)
				{
					if (b > 1) return 0;
					else	   return f * k;
				}

				v = 0;
				while ((a & 1) == 0)
				{
					v++;
					a >>= 1;
				}
				if ((v & 1) == 1) k *= tab2[b & 7];	// k *= (-1)**((b*b-1)/8)
				// step4:
				if ((a & b & 2) != 0) k = -k;		// k = k*(-1)**((a-1)*(b-1)/4)
				int r = a;
				a = b % r;
				b = r;
				// goto step3;
			}
			while (true);
		} // Kronecker int int

		private static int Kronecker(MultiSI a, int b) // return Kronecker symbol (a/b)
		//		C++ author = Joerg Arndt		email: arndt@jjj.de
		//		the C++ software is online at	http://www.jjj.de/
		// int kronecker(umod_t a, umod_t b)
		//
		// cf. Cohen p.29
		// return Kronecker symbol (a/b)
		// which, if b is an odd prime, == legendre symbol (a/b)
		//
		// From UBASIC:
		//	Extended Kronecker's symbol for m and n.  Equivalent to
		//	Legendre's symbol (for quadratic residue) for odd
		//	prime n.
		//	That is, 1 if m is a quadratic residue modulo n, -1 if
		//	m is a non-quadratic residue, 0 if m can be divided by
		//	n.
		//	Also equivalent to Jacobi's symbol for odd composite n.
		//
		//	In:   integer.
		//	Out:  integer (0, 1, or -1).
		//	Ex:   print kron(3,5)	result:-1.
		{
			int returnV;
			if (a == -1) return KroM1(b);
			int f = 1;		// final factor
			MultiSI aS = new MultiSI(a.N);
			aS &= a;
			if (aS < MultiID.sZero)
			{
				f = KroM1(b);
				aS = -aS;
			}
			if (b < 0) b = -b;
			int[] tab2 = {0, 1, 0, -1, 0, -1, 0, 1};
			// (-1)^((a^2-1)/8) == tab2[a & 7]
			if (b == 0)
			{
				if (aS != 1) returnV = 0;
				else		 returnV = f;
				goto ExitKronecker2;
			}
			if (aS == 0)
			{
				if (b != 1)	returnV = 0;
				else		returnV = f;
				goto ExitKronecker2;
			}
			// if a and b are both even
			int ai = a.V[0];
			if ((ai & 1) == 0 && (b & 1) == 0) return 0;
			int v = 0;
			while ((b & 1) == 0)
			{
				v++;
				b >>= 1;
			}
			int k;
			if ((v & 1) == 0) k = 1;
			else			  k = tab2[ai & 7];
			MultiSI r = new MultiSI();
			MultiSI bS = new MultiSI(2);
			bS &= b;
			int bi = b;
			//step3:
			do
			{
				if (aS == 0)
				{
					if (bS > MultiID.sOne)	returnV = 0;
					else					returnV = f * k;
					goto ExitKronecker1;
				}
				v = 0;
				while ((ai & 1) == 0)
				{
					if (MultiID.muErr || MultiID.muAbort) return 0; // prevent looping
					v++;
					aS /= 2;
					ai = aS.V[0];
				}
				if ((v & 1) == 1) k *= tab2[bi & 7]; // k *= (-1)**((b*b-1)/8)
				// step4:
				if ((ai & bi & 2) != 0) k = -k;		 // k = k*(-1)**((a-1)*(b-1)/4)
				r  &= aS;			// int r = a;
				aS &= (bS % r);		// a = b % r;
				bS &= r;			// b = r;
				bi = bS.V[0];
				ai = aS.V[0];
				// goto step3;
			}
			while (true);
		ExitKronecker1:
			MultiID.DoneSI(bS);
			MultiID.DoneSI(r);
		ExitKronecker2:
			MultiID.DoneSI(aS);
			return returnV;
		} // Kronecker SI int

		private static int Kronecker(MultiSI a, MultiSI bIn) // return Kronecker symbol (a/b)
		{
			MultiSI bS = new MultiSI(bIn.N);
			bS &= bIn;
			int returnV;
			if (a == -1) return KroM1(bS);
			int f = 1;		// final factor
			MultiSI aS = new MultiSI(a.N);
			aS &= a;
			if (aS < MultiID.sZero)
			{
				f = KroM1(bS);
				aS = -aS;
			}
			bS.S = false;	// if (b < 0) b = -b;
			int[] tab2 = {0, 1, 0, -1, 0, -1, 0, 1};
			// (-1)^((a^2-1)/8) == tab2[a & 7]
			if (bS == 0)
			{
				if (aS != 1) returnV = 0;
				else		 returnV = f;
				goto ExitKronecker2;
			}
			if (aS == 0)
			{
				if (bS != 1) returnV = 0;
				else		 returnV = f;
				goto ExitKronecker2;
			}
			// if a and b are both even
			int ai = a.V[0];
			if ((ai & 1) == 0 && (bS.V[0] & 1) == 0) return 0;
			int v = 0;
			while ((bS.V[0] & 1) == 0)	// b is even
			{
				v++;
				bS /= 2;		// b >>= 1;
			}
			int k;
			if ((v & 1) == 0) k = 1;
			else			  k = tab2[ai & 7];
			MultiSI r = new MultiSI();
			int bi = bS.V[0];		// int bi = b;
			//step3:
			do
			{
				if (aS == 0)
				{
					if (bS > MultiID.sOne)	returnV = 0;
					else					returnV = f * k;
					goto ExitKronecker1;
				}
				v = 0;
				while ((ai & 1) == 0)
				{
					if (MultiID.muErr || MultiID.muAbort) return 0; // prevent looping
					v++;
					aS /= 2;
					ai = aS.V[0];
				}
				if ((v & 1) == 1) k *= tab2[bi & 7]; // k *= (-1)**((b*b-1)/8)
				// step4:
				if ((ai & bi & 2) != 0) k = -k;		 // k = k*(-1)**((a-1)*(b-1)/4)
				r  &= aS;			// int r = a;
				aS &= (bS % r);		// a = b % r;
				bS &= r;			// b = r;
				bi = bS.V[0];
				ai = aS.V[0];
				// goto step3;
			}
			while (true);
		ExitKronecker1:
			MultiID.DoneSI(bS);
			MultiID.DoneSI(r);
		ExitKronecker2:
			MultiID.DoneSI(aS);
			return returnV;
		} // Kronecker SI SI

		private static int GCD(int a, int b)
		// returns the Greatest Common Divisor of a and b
		{
			a = Math.Abs(a);
			b = Math.Abs(b);
			while (true)
			{
				if (b == 0) return a;
				a %= b;
				if (a == 0) return b;
				b %= a;
			}
		} // GCD int

		private static long GCD(long a, long b)
		// returns the Greatest Common Divisor of a and b
		{
			a = Math.Abs(a);
			b = Math.Abs(b);
			while (true)
			{
				if (b == 0) return a;
				a %= b;
				if (a == 0) return b;
				b %= a;
			}
		} // GCD long

		private static MultiSI GCD(MultiSI a, MultiSI b)
		// returns the Greatest Common Divisor of a and b
		{
			MultiSI z = new MultiSI(Math.Max(a.M, b.M));
			MultiID.GCDSI(a, b, z);
			return z;
		} // GCD SI

		private static int ExEuclidGCD(int u, int v, out int u1, out int u2)
		// extended Euclid algorithm
		// see Knuth vol. 2, page 302.
		// returns gcd of u and v
		// computes u1, u2 such that u1 * u + u2 * v = gcd
		{
			int u3, v1, v2, v3, t1, t2, t3, q;
			u1 = 1; u2 = 0; u3 = u;
			v1 = 0; v2 = 1; v3 = v;
			while (v3 != 0)
			{
				q = u3 / v3;
				t1 = u1 - v1 * q;
				t2 = u2 - v2 * q;
				t3 = u3 - v3 * q;
				u1 = v1; u2 = v2; u3 = v3;
				v1 = t1; v2 = t2; v3 = t3;
			}
			return u3;
		} // ExEuclidGCD int

		private static MultiSI ExEuclidGCD(MultiSI u, MultiSI v, out MultiSI u1, out MultiSI u2)
		// extended Euclid algorithm
		{
			int n = Math.Max(u.N, v.N);
			u1 = new MultiSI(n);
			u2 = new MultiSI(n);
			MultiSI u3 = new MultiSI(n);
			MultiSI v1 = new MultiSI(n);
			MultiSI v2 = new MultiSI(n);
			MultiSI v3 = new MultiSI(n);
			MultiSI t1 = new MultiSI(n);
			MultiSI t2 = new MultiSI(n);
			MultiSI t3 = new MultiSI(n);
			MultiSI q = new MultiSI(n);
			u1 &= 1;
			u2 &= 0;
			u3 &= u;
			v1 &= 0;
			v2 &= 1;
			v3 &= v;
			while (v3 != 0)
			{
				q &= u3 / v3;
				t1 &= u1 - v1 * q;
				t2 &= u2 - v2 * q;
				t3 &= u3 - v3 * q;
				u1 &= v1;
				u2 &= v2;
				u3 &= v3;
				v1 &= t1;
				v2 &= t2;
				v3 &= t3;
			}
			MultiID.DoneSI(q);
			MultiID.DoneSI(t3);
			MultiID.DoneSI(t2);
			MultiID.DoneSI(t1);
			MultiID.DoneSI(v3);
			MultiID.DoneSI(v2);
			MultiID.DoneSI(v1);
			return u3;
		} // ExEuclidGCD SI

		private static int ModInv(int a, int n)
		//	Linear Congruences, a inverse mod n
		//	returns x such that (a * x - 1) % n = 0, i.e. a * x = 1 modulo n, x = 1/a modulo n
		//	It returns 0 if no such x exists, i.e. if a and n are
		//	not coprime.
		//	In:		integer (n > 0).
		//	Out:	integer.
		//	Ex:		print modinv(23,45)		result: 2
		{
			int gcd, x, u2;
			gcd = ExEuclidGCD(a, n, out x, out u2);
			if (gcd != 1) return 0;
			if (x < 0) x += n;
			return x;
		} // ModInv int

		private static MultiSI ModInv(MultiSI a, MultiSI n)
		//	Gives x such that a * x mod n = 1,  0 < x < n. a = 1234567901296543209877321, n = 1111111111139
		{
			MultiSI gcd, x, u2;
			gcd = ExEuclidGCD(a, n, out x, out u2);
			if (gcd != 1) return MultiID.sZero;
			if (x < MultiID.sZero) x += n;
			return x;
		} // ModInv SI

		private static int Int(MultiSI z)
		{
			int i;
			MultiID.Get1SI(out i, z);
			return i;
		}

		private static long Long(MultiSI z)
		{
			long l;
			MultiID.GetExSI(out l, z);
			return l;
		}

		private static double Dbl(MultiSI z)
		{
			double d;
			MultiID.GetDbSI(out d, z);
			return d;
		}

		private static MultiSI SI(int i)
		{
			MultiSI z = new MultiSI(MultiID.sdInInt);
			MultiID.SetTo1SI(i, z);
			return z;
		}

		private static MultiSI SI(long ex)
		{
			MultiSI z = new MultiSI(MultiID.sdInLong);
			MultiID.SetToExSI(ex, z);
			return z;
		}

		private static MultiSI SI(double d)
		{
			MultiSI z = new MultiSI(3 * MultiID.sdInDbl);
			MultiID.SetToDbSI(d, z);
			return z;
		}

		private static bool Even(MultiSI z)
		{
			return ((z.V[0] & 1) == 0);
		}

		private static bool Odd(MultiSI z)
		{
			return ((z.V[0] & 1) == 1);
		}

		private static int Bit(int n, int x) // bit n of x, 0 or 1 (0 <= n <= 63)
		{
			if (n > 30) return Bit(n, (long)x);
			int b = 1 << n;
			if ((b & x) == 0) return 0;
			return 1;
		}

		private static int Bit(int n) // bit n of saveBits[], 0 or 1
		{
			int r;
			int i = Math.DivRem(n, sdBits, out r);
			int b = 1 << r;
			if ((b & savedBits[i]) == 0) return 0;
			return 1;
		}

		private static int Bit(int n, long x) // bit n of x, 0 or 1 (0 <= n <= 63)
		{
			long b = (long)1 << n;
			if ((b & x) == 0) return 0;
			return 1;
		}

		private static int Bit(int n, MultiSI z) // bit n of x, 0 or 1
		{
			int returnV;
			long x;
			if (z <= SI(Int64.MaxValue))
			{
				MultiID.GetExSI(out x, z);
				return Bit(n, x);
			}
			int m = Math.Max(1 + n / sdBits, z.N);
			MultiSI t = new MultiSI(m);
			MultiID.PowSI(MultiID.sTwo, SI(n), t);
			MultiID.DivSI(z, t, t);
			returnV = Bit(0, t.V[0]);
			MultiID.DoneSI(t);
			return returnV;
		}

		private static int Len(int n) // number of significant bits in abs(n)
		{
			if (n == 0) return 0;
			n = Math.Abs(n);
			double d = Math.Log(n, 2.0);
			int i = 1 + (int)Math.Floor(d); // Floor could be one too small if log2(2^p) = p - epsilon
			int max = (1 << i) - 1;
			if (n > max) i++;
			return i;
		}

		private static int Len(MultiSI z) // number of significant bits in z
		{
			savedBits = new int[(int)Math.Ceiling(1.03 * z.N)];
			int m;
			int len = 0;
			int i = 0;
			if (z.N == 1) m = z.V[0];
			else
			{
				MultiSI q = new MultiSI();
				MultiID.SetToSI(z, q);
				q.S = false;
				for (i = 0; q.N > 1; i++)
				{
					MultiID.Div1QRSI(q, twoToSDB, out savedBits[i], q);
					len += sdBits;
				}
				m = q.V[0];
				MultiID.DoneSI(q);
			}
			len += Len(m);
			if (m < twoToSDB) savedBits[i] = m;
			else			  savedBits[i + 1] = Math.DivRem(m, twoToSDB, out savedBits[i]);
			return len;
		}

		private static void PrimeDiv(MultiSI n, ref long d) // d = the least prime divisor of n >= d
		// returns 0 if n is greater than stopDiv^2 and has no divisor less than stopDiv
		{
			if (MultiID.ZTestSI(n))
			{
				d = 0;
				return;
			}
			if ( ((n.V[0] & 1) == 0) && (d <= 2))
			{
				d = 2;
				return;
			}
			long dX = d - 1;
			if (dX < 2) dX = 2;				// 2 has already been tested
			MultiSI q = new MultiSI(n.N);
			int rN = 1 + 10 / MultiID.muDMax;
			MultiSI r = new MultiSI(rN);
			MultiSI dSI = new MultiSI(rN);
			do
			{
				if (MultiID.statReq)
				{
					MultiID.WriteLn("PrimeDiv: Testing divisor " + dX);
					MultiID.statReq = false;
				}
				NextP(ref dX);
				MultiID.SetToExSI(dX, dSI);
				MultiID.DivQRSI(n, dSI, r, q);
				if (MultiID.ZTestSI(r))
				{
					d = dX;
					dX = 0;
					break;						// return prime d as a divisor
				}
				if (q < dSI)
				{
					MultiID.GetExSI(out d, n);	// n is prime so return d = n
					dX = 0;
					break;
				}
			} while (dX < stopDiv);
			if (dX >= stopDiv) d = 0;			// divisor not found
			MultiID.DoneSI(dSI);
			MultiID.DoneSI(r);
			MultiID.DoneSI(q);
		} // PrimeDiv

		public static void Prime(int nI) // return n = 1 if |n| is prime, else 0
		{
			MultiID.ADDiag("Is x a prime number?:");
			int i = Prime(Common.nList[nI]);
			Calc.ClearN(nI);
			Common.nList[nI].V[0] = i;
		}

		public static int Prime(MultiSI z) // return n = 1 if |z| is prime, else 0
		{
			if (!MultiID.GT1UI(1, z))	return 0;	// -1, 0, 1 => 0
			if (!MultiID.GT1UI(3, z))	return 1;	// -3, -2, 2, 3 => 1
			if ((z.V[0] & 1) == 0)		return 0;	// remaining even => 0
			bool sign = z.S;
			z.S = false;
			int i = FnAdleman(z);
			z.S = sign; // restore sign
			if (MultiID.muErr || MultiID.muAbort) return 0;
			return i;
		}

		public static void PrimePiLT(int nI) // n = PiL(n) =
		// Pi(n) = number of primes <= n by Lehmer's formula (Top level called by calculator)
		{
			if (!pTabFilled) FillPTab();
			PrimePiL(nI);
		}

		public static void PrimePiL(int nI) // n = PiL(n) =
		// Pi(n) = number of primes <= n by Lehmer's formula
		{
			PiLDepth++;
			if (MultiID.CompSI(MultiID.sTwo, Common.nList[nI]) < 0)	// if n < 2
			{
				Calc.ClearN(nI);
				return;
			}
			int aI;
			Common.InitNAdd(1, "PrimPiL-aI-" + PiLDepth, out aI);
			int bI;
			Common.InitNAdd(1, "PrimPiL-bI-" + PiLDepth, out bI);
			XICommon.SqRtN(nI, bI);			// b = SqRt(n)
			XICommon.SqRtN(bI, aI);			// a = SqRt(SqRt(n))
			XICommon.PrimePi(aI);			// a = Pi(SqRt(SqRt(n)))
			long a;
			MultiID.GetExSI(out a, Common.nList[aI]);
			XICommon.PrimePi(bI);			// b = Pi(SqRt(n))
			long b;
			MultiID.GetExSI(out b, Common.nList[bI]);
			int cI;
			Common.InitNAdd(1, "PrimPiL-cI-" + PiLDepth, out cI);
			XICommon.CubeRtN(nI, cI);		// c = CuRt(n)
			XICommon.PrimePi(cI);			// c = Pi(CuRt(n))
			long c;
			MultiID.GetExSI(out c, Common.nList[cI]);
			int vI;
			Common.InitNAdd(Common.nList[nI].N, "PrimPiL-vI-" + PiLDepth, out vI);
			if (a == 0)
				Calc.SetToN(Common.nList[nI], vI);					// v = n
			else
				PhiLSI(Common.nList[nI], (int)a, Common.nList[vI]);	// v = PhiL(n, a)
			int tI;
			Common.InitNAdd(1, "PrimPiL-tI-" + PiLDepth, out tI);
			Calc.SetToN(MultiID.sTwo, tI);
			Common.nList[tI].S = true;			// t = -2
			XICommon.RAddN(bI, tI);					// t = b - 2
			XICommon.RAddN(aI, tI);					// t = b + a - 2
			int uI;
			Common.InitNAdd(1, "PrimPiL-uI-" + PiLDepth, out uI);
			Calc.SetToN(MultiID.sOne, uI);		// u = 1
			XICommon.RAddN(bI, uI);					// u = b + 1
			XICommon.RSubN(aI, uI);					// u = b - a + 1
			XICommon.MulN(tI, uI, tI);				// t = (b + a - 2)*(b - a + 1)
			MultiID.ReSizeSI(MultiID.sdInLong, Common.nList[uI]);
			MultiID.DivSI(Common.nList[tI], MultiID.sTwo, Common.nList[tI]); // t = t / 2
			XICommon.RAddN(tI, vI);					// v = PhiL(n, a) + t
			MultiSI nd = new MultiSI(Common.nList[nI].N);
			for (long i = a+1; i <= b; i++)		// for PiL(10^10), b got larger than lenPTab
			{
				if (MultiID.muAbort || MultiID.muErr) break; // exit for
				if (i <= lenPTab)
					MultiID.Div1QRSI(Common.nList[nI], pTab[i], out MultiID.rem1, nd);				// nd =n/p(i)
				else if (i <= NPrimes.jLimM)
					MultiID.Div1QRSI(Common.nList[nI], (int)NPrimes.iP[i], out MultiID.rem1, nd);	// nd =n/p(i)
				else
				{
					MultiID.SetToExSI(i, Common.nList[uI]);
					XICommon.PrimeP(uI);
					MultiID.DivSI(Common.nList[nI], Common.nList[uI], nd);	// nd =n/p(i)
				}
				Calc.SetToN(nd, tI);					// t = n/p(i)
				XICommon.PrimePi(tI);						// t = Pi(n/p(i))
				XICommon.RSubN(tI, vI);						// subtract t from Pi(nI)
				if (MultiID.statReq2)
				{
					MultiID.nDone = i - a;
					MultiID.nToGo = b - i;
					MultiID.WritenDone("PrimePiL:");
				}
			}
			MultiSI n2 = new MultiSI(Common.nList[nI].N);
			MultiSI jm = new MultiSI(MultiID.sdInLong);
			for (long i = a+1; i <= c; i++)
			{
				if (MultiID.muAbort || MultiID.muErr) break; // exit for
				if (i <= lenPTab)
					MultiID.Div1QRSI(Common.nList[nI], pTab[i], out MultiID.rem1, nd);				// nd =n/p(i)
				else if (i <= NPrimes.jLimM)
					MultiID.Div1QRSI(Common.nList[nI], (int)NPrimes.iP[i], out MultiID.rem1, nd);	// nd =n/p(i)
				else
				{
					MultiID.SetToExSI(i, Common.nList[uI]);
					XICommon.PrimeP(uI);
					MultiID.DivSI(Common.nList[nI], Common.nList[uI], nd);	// nd =n/p(i)
				}
				Calc.SetToN(nd, bI);					// bi = n/p(i)
				XICommon.SqRtN(bI, bI);						// bi = SqRt(n/p(i))
				XICommon.PrimePi(bI);						// bi = Pi(n/p(i))
				long bi;
				MultiID.GetExSI(out bi, Common.nList[bI]);
				for (long j = i; j <= bi; j++)
				{
					if (j <= lenPTab)
						MultiID.Div1QRSI(nd, pTab[j], out MultiID.rem1, n2);			// n2 = nd/p(j)
					else if (j <= NPrimes.jLimM)
						MultiID.Div1QRSI(nd, (int)NPrimes.iP[j], out MultiID.rem1, n2);	// n2 = nd/p(j)
					else
					{
						MultiID.SetToExSI(j, Common.nList[uI]);
						XICommon.PrimeP(uI);
						MultiID.DivSI(nd, Common.nList[uI], n2);	// n2 = nd/p(j)
					}
					Calc.SetToN(n2, tI);							// t = n/(p(i)*p(j))
					XICommon.PrimePi(tI);								// t = Pi(n/(p(i)*p(j))
					XICommon.RSubN(tI, vI);								// subtract t from Pi(nI)
					MultiID.SetToExSI(j - 1, jm);					// jm = j-1
					XICommon.RAddSIN(jm, vI);							// + (j-1)
				}
				if (MultiID.statReq2)
				{
					MultiID.nDone = i - a;
					MultiID.nToGo = c - i;
					MultiID.WritenDone("PrimePiL:");
				}
			}
			MultiID.DoneSI(jm);
			MultiID.DoneSI(n2);
			MultiID.DoneSI(nd);
			Common.LDelete(uI);
			Common.LDelete(tI);
			Calc.SetToN(Common.nList[vI], nI);
			Common.LDelete(vI);
			Common.LDelete(cI);
			Common.LDelete(bI);
			Common.LDelete(aI);
			PiLDepth--;
		} // PrimePiL

		public static void PrimePiM(int nI) // n = PiM(n) =
		// Pi(n) = number of primes <= n by Meissel's formula
		{
			if (MultiID.CompSI(MultiID.sTwo, Common.nList[nI]) < 0)	// if n < 2
			{
				Calc.ClearN(nI);
				return;
			}
			int bI;
			Common.InitNAdd(1, "PrimPiM-bI", out bI);
			XICommon.SqRtN(nI, bI);			// b = SqRt(n)
			XICommon.PrimePi(bI);			// b = Pi(SqRt(n))
			long b;
			MultiID.GetExSI(out b, Common.nList[bI]);
			int cI;
			Common.InitNAdd(1, "PrimPiM-cI", out cI);
			XICommon.CubeRtN(nI, cI);		// c = CuRt(n)
			XICommon.PrimePi(cI);			// c = Pi(CuRt(n))
			long c;
			MultiID.GetExSI(out c, Common.nList[cI]);
			int vI;
			Common.InitNAdd(Common.nList[nI].N, "PrimPiM-vI", out vI);
			PhiLSI(Common.nList[nI], (int)c, Common.nList[vI]);	// v = PhiL(n, c)
			int tI;
			Common.InitNAdd(1, "PrimPiM-tI", out tI);
			Calc.SetToN(MultiID.sTwo, tI);
			Common.nList[tI].S = true;			// t = -2
			XICommon.RAddN(bI, tI);					// t = b - 2
			XICommon.RAddN(cI, tI);					// t = b + c - 2
			int uI;
			Common.InitNAdd(1, "PrimPiM-uI", out uI);
			Calc.SetToN(MultiID.sOne, uI);		// u = 1
			XICommon.RAddN(bI, uI);					// u = b + 1
			XICommon.RSubN(cI, uI);					// u = b - c + 1
			XICommon.MulN(tI, uI, tI);				// t = (b + c - 2)*(b - c + 1)
			MultiID.ReSizeSI(MultiID.sdInLong, Common.nList[uI]);
			MultiID.DivSI(Common.nList[tI], MultiID.sTwo, Common.nList[tI]); // t = t / 2
			XICommon.RAddN(tI, vI);					// v = PhiL(n, c) + t
			MultiSI nd = new MultiSI(Common.nList[nI].N);
			for (long i = c+1; i <= b; i++)
			{
				if (MultiID.muAbort || MultiID.muErr) break; // exit for
				if (i <= lenPTab)
					MultiID.Div1QRSI(Common.nList[nI], pTab[i], out MultiID.rem1, nd);				// nd =n/p(i)
				else if (i <= NPrimes.jLimM)
					MultiID.Div1QRSI(Common.nList[nI], (int)NPrimes.iP[i], out MultiID.rem1, nd);	// nd =n/p(i)
				else
				{
					MultiID.SetToExSI(i, Common.nList[uI]);
					XICommon.PrimeP(uI);
					MultiID.DivSI(Common.nList[nI], Common.nList[uI], nd);	// nd =n/p(i)
				}
				Calc.SetToN(nd, tI);					// t = n/p(i)
				XICommon.PrimePi(tI);						// t = Pi(n/p(i))
				XICommon.RSubN(tI, vI);						// subtract t from Pi(nI)
				if (MultiID.statReq2)
				{
					MultiID.nDone = i - c;
					MultiID.nToGo = b - i;
					MultiID.WritenDone("PrimePiM:");
				}
			}
			MultiID.DoneSI(nd);
			Common.LDelete(uI);
			Common.LDelete(tI);
			Calc.SetToN(Common.nList[vI], nI);
			Common.LDelete(vI);
			Common.LDelete(cI);
			Common.LDelete(bI);
		} // PrimePiM

		public static void PrimePiL1(int nI) // n = PiL1(n) =
		// Pi(n) = number of primes <= n by Legendre's formula
		{
			if (MultiID.CompSI(MultiID.sTwo, Common.nList[nI]) < 0)	// if n < 2
			{
				Calc.ClearN(nI);
				return;
			}
			int pI;
			Common.InitNAdd(1, "PrimPiL-pI", out pI);
			XICommon.SqRtN(nI, pI);			// p = Sqrt(n)
			XICommon.PrimePi(pI);			// p = Pi(Sqrt(n))
			PhiLN(nI, pI, nI);
			XICommon.RAddN(pI, nI);
			Calc.SetToN(MultiID.sOne, pI);
			XICommon.RSubN(pI, nI);
			Common.LDelete(pI);
		}

		public static void PhiLN(int nI, int aI, int pI) // p = PhiL(n, a) =
		// Legendre's formula i.e., the number of integers in [1, n]
		// not divisible by any of the first a primes.
		{
			if (Common.nList[aI].S)
			{
				Common.nList[pI] &= Common.nList[nI];
				Calc.NormN(pI);
				return;
			}
			if (!pTabFilled) FillPTab();
			long a = Int64.MaxValue;
			if (Common.nList[aI] <= MultiID.max64)
				MultiID.GetExSI(out a, Common.nList[aI]);
			if (a > Int32.MaxValue)
			{
				MultiID.MuWriteErr("a > " + Int32.MaxValue + " too large for Legendre's formula");
				return;
			}
			if (a == 0)
				Calc.SetToN(Common.nList[nI], pI);
			else
				PhiLSI(Common.nList[nI], (int)a, Common.nList[pI]);
			if (MultiID.muAbort || MultiID.muErr) Calc.ClearN(pI);
			Calc.NormN(pI);
		} // PhiLN

		public static void PhiLInt(int n, int a, out int p) // p = PhiL(n, a) =
		// Legendre's formula i.e., the number of integers in [1,n]
		// not divisible by any of the first a primes.
		{
			if (MultiID.muAbort) { p = 0; return;}
			int nd;
			if (MultiID.statReq2)
			{
				int tot = 0;
				MultiID.WriteT("PhiLInt: n = " + n + ", a = " + a, ref tot);
				MultiID.TimeDel(ref tot);
				MultiID.statReq2 = false;
			}
			if (a < 2)
			{
				MultiID.MuInterrupt();
				p = n - n/2; // assume a >= 1
			}
			else
			{
				if (a <= lenPTab)
					nd = n / pTab[a];
				else if (a <= NPrimes.jLimM)
					nd = n / (int)NPrimes.iP[a];
				else
				{
					int uI;
					Common.InitNAdd(MultiID.sdInLong, "PhiLEx-uI", out uI);
					MultiID.SetToExSI(a, Common.nList[uI]);
					XICommon.PrimeP(uI);
					MultiID.GetIntSI(out nd, Common.nList[uI]);
					nd = n / nd;							// nd = n / p(a)
					Common.LDelete(uI);
				}
				int am = a - 1;
				int p1;
				PhiLInt(n, am, out p1);						// p1 = PhiL(n, a-1)
				if (nd < 1)									// if nd < 1 i.e. zero
					p = p1;
				else
				{
					int p2;
					PhiLInt(nd, am, out p2);				// p2 = PhiL(n/p(a), a-1)
					p = p1 - p2;
				}
			}
		} // PhiLInt

		public static void PhiLEx(long n, int a, out long p) // p = PhiL(n, a) =
		// Legendre's formula i.e., the number of integers in [1,n]
		// not divisible by any of the first a primes.
		{
			if (n <= int.MaxValue)
			{
				int pi;
				PhiLInt((int)n, a, out pi);
				p = pi;
				return;
			}
			long nd;
			if (MultiID.muAbort) { p = 0; return; }
			if (MultiID.statReq2)
			{
				int tot = 0;
				MultiID.WriteT("PhiLEx: n = " + n + ", a = " + a, ref tot);
				MultiID.TimeDel(ref tot);
				MultiID.statReq2 = false;
			}
			if (a < 2)
			{
				MultiID.MuInterrupt();
				p = n - n/2; // assume a >= 1
			}
			else
			{
				if (a <= lenPTab)
					nd = n / pTab[a];
				else if (a <= NPrimes.jLimM)
					nd = n / (long)NPrimes.iP[a];
				else
				{
					int uI;
					Common.InitNAdd(MultiID.sdInLong, "PhiLEx-uI", out uI);
					MultiID.SetToExSI(a, Common.nList[uI]);
					XICommon.PrimeP(uI);
					MultiID.GetExSI(out nd, Common.nList[uI]);
					nd = n / nd;							// nd = n / p(a)
					Common.LDelete(uI);
				}
				int am = a - 1;
				long p1;
				PhiLEx(n, am, out p1);						// p1 = PhiL(n, a-1)
				if (nd < 1)									// if nd < 1 i.e. zero
					p = p1;
				else
				{
					long p2;
					PhiLEx(nd, am, out p2);					// p2 = PhiL(n/p(a), a-1)
					p = p1 - p2;
				}
			}
		} // PhiLEx

		public static void PhiLSI(MultiSI n, int a, MultiSI p) // p = PhiL(n, a) =
		// Legendre's formula i.e., the number of integers in [1,n]
		// not divisible by any of the first a primes.
		{
			if (n <= MultiID.max64)
			{
				long ni, pi;
				MultiID.GetExSI(out ni, n);
				PhiLEx(ni, a, out pi);
				MultiID.SetToExSI(pi, p);
				return;
			}
			MultiSI nd;
			if (MultiID.muAbort || MultiID.muErr) return;
			PhiLSIDepth++;
			if (MultiID.statReq2)
			{
				int tot = 0;
				MultiID.WriteT("PhiLSI: n = ", ref tot);
				MultiID.WriteSIF(n, ref tot);
				MultiID.WriteT(", a = " + a + ", depth = " + PhiLSIDepth, ref tot);
				MultiID.TimeDel(ref tot);
				MultiID.statReq2 = false;
			}
			if (a < 2)
			{												// p = n - n/2, assume a >= 1
				if ((n.V[0] & 1) == 0)						// if n is even, assume a >= 1
					MultiID.DivSI(n, MultiID.sTwo, p);		// p = n/2
				else
				{
					MultiID.DivSI(n, MultiID.sTwo, p);		// p = n/2
					MultiID.RAddSI(MultiID.sOne, p);		// p = 1 + n/2
				}
			}
			else
			{
				nd = new MultiSI(n.N);
				int uI;
				Common.InitNAdd(MultiID.sdInLong, "PhiLSI-uI", out uI);

				if (a <= lenPTab)
					MultiID.Div1QRSI(n, pTab[a], out MultiID.rem1, nd);				// nd = n/p(a)
				else if (a <= NPrimes.jLimM)
					MultiID.Div1QRSI(n, (int)NPrimes.iP[a], out MultiID.rem1, nd);	// nd = n/p(a)
				else
				{
					MultiID.SetToExSI(a, Common.nList[uI]);
					XICommon.PrimeP(uI);
					MultiID.DivSI(n, Common.nList[uI], nd);	// nd = n/p(a)
				}
				int am = a - 1;
				MultiSI p1 = new MultiSI(n.N);
				PhiLSI(n, am, p1);							// p1 = PhiL(n, a-1)
				if (nd.V[0] == 0 && nd.N == 1)				// if nd < 1 i.e. zero
					MultiID.SetToSI(p1, p);
				else
				{
					MultiSI p2 = new MultiSI(nd.N);
					PhiLSI(nd, am, p2);						// p2 = PhiL(n/p(a), a-1)
					MultiID.SubSI(p1, p2, p);				// p = p1 - p2
					MultiID.DoneSI(p2);
				}
				MultiID.DoneSI(p1);
				Common.LDelete(uI);
				MultiID.DoneSI(nd);
			}
			PhiLSIDepth--;
		} // PhiLSI

		public static void SorderSI(MultiSI a, MultiSI n, MultiSI e) // e = Sord(a, n) =
		// multiplicative suborder of base a (mod n) or zero if it does not exist.
		// From Stephen Wolfram, Algebraic Properties of Cellular Automata (1984).
		{
			MultiSI ee = new MultiSI(e.N);	// Sord = Mord or MORD / 2
			MorderSI(a, n, ee);				// = Mord / 2 iff a^e mod n = -1
			if ((ee.V[0] & 1) == 1 || ee == MultiID.sZero)
			{
				e &= ee;					// if ee is odd or zero, return ee
				MultiID.DoneSI(ee);
				return;
			}
			MultiSI g1 = new MultiSI(2 * Math.Max(a.N, n.N));
			ee /= 2;
			MultiID.PowMBSI(a, ee, n, g1);			// g1 = a^ee mod n
			if (g1 != (n-1)) ee *= 2;				// test a^(Mord(a,n)/2) mod n = -1
			e &= ee;
			MultiID.DoneSI(g1);
			MultiID.DoneSI(ee);
		}

		public static void MorderSI(MultiSI a, MultiSI n, MultiSI e) // e = Mord(a, n) =
		// multiplicative order of base a (mod n) or zero if it does not exist. From Henri Cohen's book,
		// A Course in Computational Algebraic Number Theory, Springer, 1996, Algorithm 1.4.3, page 25.
		{
			if (a <= MultiID.sOne || n <= MultiID.sOne)
			{
				e &= 0;
				return;
			}
			if (GCD(a, n) != MultiID.sOne)
			{
				e &= 0;
				return;
			}
			MultiSI f = new MultiSI(n.N);
			f &= n;
			PrimePhi(f);						// f = Phi(n)
			ECM(f);								// factor f by ECM, f = Prod(pi^vi) i = 1, ..., nF
			if (!complete)
			{
				e &= 0;
				MultiID.DoneSI(f);
				return;
			}
			MultiSI ee = new MultiSI(f.N);
			MultiSI g1 = new MultiSI(2 * Math.Max(a.N, n.N));
			MultiSI vi = new MultiSI(MultiID.sdInDbl);
			ee &= f;								// ee = Phi(n) 
			for (int i = 1; i <= nF; i++)
			{
				MultiID.SetToDbSI(powers.a[i], vi);
				MultiID.PowSI(primes.a[i], vi, f);	// f = pi^vi
				ee &= (ee / f);						// ee = ee / (pi^vi)
				MultiID.PowMBSI(a, ee, n, g1);		// g1 = a^ee mod n
				while (g1 != MultiID.sOne)
				{
					MultiID.PowMBSI(g1, primes.a[i], n, g1);	// g1 = g1^pi mod n
					MultiID.MulSI(ee, primes.a[i], ee);			// ee = ee * pi
					if (MultiID.muErr || MultiID.muAbort) break;
				}
				if (MultiID.muErr || MultiID.muAbort) break;
			}
			if (MultiID.muErr || MultiID.muAbort) ee &= 0;
			e &= ee;
			MultiID.DoneSI(vi);
			MultiID.DoneSI(g1);
			MultiID.DoneSI(ee);
			MultiID.DoneSI(f);
		} // MorderSI

		public static void PrimRQ(int xI, int yI, int zI) // return z = 1 (true) if x is a primitive root of y, else 0 (false)
		{
			MultiID.ADDiag("PrimRQ function, is x a primitive root of y?:");
			bool so = false;								// set up to return 0 (false)
			bool signY = Common.nList[yI].S;				// save sign of y
			Common.nList[yI].S = false;								// don't use sign of y
			MultiSI g = new MultiSI(Common.nList[yI].N);
			MultiID.ModSI(Common.nList[xI], Common.nList[yI], g);	// g = x mod y
			MorderSI(g, Common.nList[yI], g);						// g = Mord(x, y) or zero
			if (!MultiID.ZTestSI(g))
			{
				int phI;
				Common.InitNAdd(1, "PrimRQ-phI", out phI);
				Common.nList[phI] &= Common.nList[yI];
				PrimePhi(phI);									// ph = Phi(y)
				so = (g == Common.nList[phI]);					// true Mord(x, y) == Phi(y)
				Common.LDelete(phI);
			}
			Common.nList[yI].S = signY;							// restore sign of y
			Calc.ClearN(zI);									// set up to return 0 (false)
			if (so) Common.nList[zI].V[0] = 1;					// return 1 (true)
			MultiID.DoneSI(g);
			Calc.NormN(zI);
		} // PrimRQ

		public static void PrimR(int xI, bool rPrime) // return x = gen = least primitive root of x modulo x or zero
		// rRrime true => return least primitive prime root of x modulo x or zero
		{
			int rI;							// put largest root on "Re"
			Common.LFindAdd("Re", out rI);
			if (xI == rI || rI == Common.iTop)
			{
				MultiID.MuWriteErr("Cannot set same location to both max and min primitive root, continuing...");
				return;
			}
			Calc.ClearN(rI);
			MultiID.ADDiag("Find the least primitive root of x:");
			if (Common.nList[xI] == MultiID.sZero)
			{
				MultiID.WriteLn("The non-prime number 0 has 0 primitive roots");
				return;	// 0 => 0
			}
			Common.nList[xI].S = false;						// don't use sign of x
			if (rPrime)
			{
				if (Common.nList[xI] == MultiID.sOne)
				{
					MultiID.WriteLn("The non-prime number 1 has 0 prime primitive roots");
					Common.nList[xI].V[0] = 0;
					return;	// 1 => 0 prime primitive roots
				}
				if (Common.nList[xI] == MultiID.sTwo)
				{
					MultiID.WriteLn("The prime number 2 has 0 prime primitive roots");
					Common.nList[xI].V[0] = 0;
					return;	// 2 => 0 prime primitive roots
				}
			}
			if (Common.nList[xI] <= MultiID.sFour)
			{
				if (Common.nList[xI] == MultiID.sOne || Common.nList[xI] == MultiID.sFour)
					MultiID.Write("The non-prime number " + Common.nList[xI].V[0]);
				else
					MultiID.Write("The prime number " + Common.nList[xI].V[0]);
				Common.nList[xI].V[0]--;	// 1 => 0,  2 => 1, 3 => 2, 4 => 3
				MultiID.Write(" has 1 ");
				if (rPrime) MultiID.Write("prime ");
				MultiID.WriteLn("primitive root equal to " + Common.nList[xI].V[0]);
				Common.nList[rI].V[0] = Common.nList[xI].V[0];
				return;
			}
			int tot = 0;
			if ((Common.nList[xI].V[0] & 3) == 0) // if x = 4 * a number
			{
					MultiID.WriteT("The non-prime number ", ref tot);
					MultiID.WriteSIFLn(Common.nList[xI], ref tot);
					MultiID.WriteTLn("has 0 primitive roots", ref tot);
					Calc.ClearN(xI);
					return;			
			}
			MultiSI ord = new MultiSI(Common.nList[xI].N);		// modulo order
			MultiSI nw = new MultiSI(Common.nList[xI].N);
			MultiSI rem = new MultiSI(Common.nList[xI].N);		// remainder
			MultiSI gen = new MultiSI(Common.nList[xI].N);		// generator
			MultiSI div = new MultiSI(2 * Common.nList[xI].N);	// divisor
			MultiSI p = new MultiSI(Common.nList[xI].N);		// a prime
			long d = 0;
			int phI;											// totient function Phi
			bool zEven = false;
			nw &= Common.nList[xI];
			if ((nw.V[0] & 1) == 0)		// if x = 2 * a number
			{
				zEven = true;
				MultiID.DivSI(nw, MultiID.sTwo, nw);
			}
			gen &= 1;
			bool onlyOne = true;							// only one primitive root flag
			bool nwPrime = (FnAdleman(nw) == 1);			// is nw prime?
			bool xPrime = (!zEven && nwPrime);				// is x prime
			Common.SingleSpace();
			if (xPrime && !MultiID.muErr && !MultiID.muAbort)
			{
				MultiID.SubSI(Common.nList[xI], MultiID.sOne, ord);	// ord = x - 1 = order = Phi(x)
				Common.InitNAdd(1, "PrimR-phI", out phI);
				Calc.SetToN(ord, phI);
				PrimePhi(phI);										// ph = Phi(Phi(x))
				tot = 0;
				MultiID.WriteT("The prime number ", ref tot);
				MultiID.WriteSILn(Common.nList[xI], ref tot);
				MultiID.WriteT("has ", ref tot);
				onlyOne = (Common.nList[phI] == 1);
				if (onlyOne)
					MultiID.WriteT("1 primitive root equal to ", ref tot);
				else
				{
					MultiID.WriteSI(Common.nList[phI], ref tot);
					MultiID.WriteT(" primitive roots, the largest ", ref tot);
					if (rPrime) MultiID.WriteT("prime ", ref tot);
					MultiID.WriteT("is ", ref tot);
				}
				Common.LDelete(phI);
			}
			else if (!MultiID.muErr && !MultiID.muAbort)
			{
				int k = 1;
				if (nwPrime) p &= nw;
				else
				{
					testPrimR = true;
					ECM(nw);			// factor nw, nw is odd and either x or x/2
					testPrimR = false;
					if (MultiID.muErr || MultiID.muAbort)
					{
						MultiID.ClearSI(gen);
						goto ExitPrimeR;
					}
					if (nF > 1)			// if more than one odd prime factor
					{
						MultiID.WriteT("The non-prime number ", ref tot);
						MultiID.WriteSILn(Common.nList[xI], ref tot);
						MultiID.WriteTLn("has 0 primitive roots", ref tot);
						MultiID.ClearSI(gen);
						goto ExitPrimeR;
					}
					p &= primes.a[1];
					k = powers.a[1];
				}
				MultiID.DivSI(nw, p, ord);
				MultiID.SubSI(nw, ord, ord);		// ord = nw - nw / p = order = Phi(nw = p^k)
				Common.InitNAdd(1, "PrimR-phI", out phI);
				Calc.SetToN(ord, phI);
				PrimePhi(phI);						// ph = Phi(Phi(p^k))
				Common.SingleSpace();
				MultiID.WriteT("The non-prime number ", ref tot);
				MultiID.WriteSI(Common.nList[xI], ref tot);
				MultiID.WriteT(" = ", ref tot);
				if (zEven) MultiID.WriteT("2 * ", ref tot);
				MultiID.WriteSI(p, ref tot);
				if (k != 1)
					MultiID.WriteTLn("^" + k, ref tot);
				else
					MultiID.WriteTLn("", ref tot);
				MultiID.WriteT("has ", ref tot);
				onlyOne = (Common.nList[phI] == 1);
				if (onlyOne)
					MultiID.WriteT("1 primitive root equal to ", ref tot);
				else
				{
					MultiID.WriteSI(Common.nList[phI], ref tot);
					MultiID.WriteT(" primitive roots, the largest ", ref tot);
					if (rPrime) MultiID.WriteT("prime ", ref tot);
					MultiID.WriteT("is ", ref tot);
				}
				Common.LDelete(phI);
			}
			if (MultiID.muErr || MultiID.muAbort) goto ExitPrimeR;
			int getN = 2;
			if (onlyOne) getN = 1;
			for (int i = 0; i < getN; i++)
			{
				int inc;
				if (onlyOne || (i > 0))
				{
					inc = 1;
					gen &= 1;			// search for smallest
				}
				else
				{
					inc = -1;
					gen &= Common.nList[xI];	// search for largest
				}
				if (i > 0)
				{
					MultiID.WriteT(", the smallest ", ref tot);
					if (rPrime) MultiID.WriteT("prime ", ref tot);
					MultiID.WriteT("is ", ref tot);
				}
				int sw = 0;
				while (sw == 0 && !MultiID.muErr && !MultiID.muAbort)
				{
					do								// increment gen
					{
						if (xPrime) gen += inc;
						else
						{
							do
							{
								gen += inc;
								MultiID.GCDSI(gen, Common.nList[xI], nw);
							}
							while (nw != 1);
						}
					}
					while (rPrime && (FnAdleman(gen) == 0));
					sw = 1;
					nw &= ord;
					stopDiv = 131072;	// 2^17 for PrimeDiv function (2^17 from UBASIC program)
					do
					{
						d = 1;
						NextP(ref d);			// NextP will call InitP() if needed
						PrimeDiv(nw, ref d);	// will return 0 if nw too large
						div &= d;
						if (d == 0)
						{
							if (FnAdleman(nw) == 0)		// returns 1 if z is prime, 0 if composite
							{
								FnECM(nw, div);			// find a factor of nw by ECM
								if (MultiID.muErr)
								{
									MultiID.WriteSLn("PrimR: number too large");
									MultiID.ClearSI(gen);
									goto ExitPrimeR;
								}
							}
							else
							{
								div &= nw;
							}
						}
						do
						{
							MultiID.RDivQRSI(div, rem, nw);	// nw = nw / div, rem = remainder
						} 
						while (rem == 0 && !MultiID.muErr && !MultiID.muAbort);
						nw *= div;
						nw += rem;		// nw = nw * div + rem
						MultiID.DivSI(ord, div, rem);
						MultiID.PowMBSI(gen, rem, Common.nList[xI], div); // div = gen^(ord/div) mod x
						if (div == 1)
						{
							sw = 0;
							nw &= 1;
						}
					}
					while (nw != 1 && !MultiID.muErr && !MultiID.muAbort);
				}
				MultiID.WriteSI(gen, ref tot);
				if (!onlyOne && (i == 0))
				{
					MultiID.WriteT(" = n - ", ref tot);
					gen &= Common.nList[xI] - gen;	// display n - largest
					MultiID.WriteSI(gen, ref tot);
					gen &= Common.nList[xI] - gen;	// restore largestPrimR(137842
				}
				if (i == 0)
				{
					Calc.SetToN(gen, rI);
					Calc.NormN(rI);
				}
			}
		ExitPrimeR:
			Calc.SetToN(gen, xI);
			MultiID.DoneSI(p);
			MultiID.DoneSI(div);
			MultiID.DoneSI(gen);
			MultiID.DoneSI(rem);
			MultiID.DoneSI(nw);
			MultiID.DoneSI(ord);
			Calc.NormN(xI);
		} // PrimeR

		private static void PrimeDoIt(int nI) // prime factor Do It
		{
			PrimeDoIt(Common.nList[nI]);
		}

		private static void PrimeDoIt(MultiSI z) // prime factor Do It
		{
			exitEarly = false;
			size = z.N;
			size = Math.Min(size, MultiID.mMax);
			size = Math.Max(size, (48 / MultiID.muDMax) + 1);	// size for 49 decimal digits
			sizeO2 = 1 + size / 2;
			Init();
			MultiID.SetToSI(z, zIn);
			MultiID.SetTo1SI(2, start);
			zIn.S = false;
			do
				NextMainNumber();
			while (!testPrime && findAPrime && !MultiID.muAbort && !exitEarly);
			if (!complete && !testPrime && !Common.quiet)
				MultiID.WriteLn("Number not completely factored");
			Fini1();
		}

		public static void KronSI(MultiSI x, MultiSI y, MultiSI z)
		{
			MultiID.SetTo1SI(Kronecker(x, y), z);
		}

		public struct TPrimes
		{
			public MultiSI[] a;				// prime factor table[maxPs+1]
		};

		public struct Tpowers
		{
			public int[] a;					// factor power table[maxPs+1]
		};

		private static bool byECM;			// flag that the input number was factored by ECM
		private static bool exitEarly;		// flag that the input number has too many factors to store
		private static bool findAPrime;		// find first prime < zIn or > zIn
		private static bool isSquare;		// flag that the factored number is a square like 100 = 2^2 * 5^2
		private static bool noMoreDiv;		// do not use MoreDiv method
		private static bool noSieve;		// do not use sieve method
		private static bool squareFree;		// flag that the factored number is squarefree like 30 = 2*3*5
		private static bool testPrime;		// flag for test if number is prime
		private static bool testPrimR;		// flag for ECM to exit early if testing for a primitive root and 2 odd primes
		private static bool testSqFree;		// flag for ECM to exit early if squareFree = false
		private static bool testSquare;		// flag for ECM to exit early if isSquare = false
		private static bool[,] sieve;		// [rS+1, mRM+1] for SieveMethod, rS tables
		private static int maxPs;			// constant, Max number of prime factors, 24 for 36 digit numbers
											// 54 for 101 digit numbers, about 1400 for 5000 digit #s
		private static int PhiLSIDepth;		// depth of calls to PhiLSI
		private static int rS;				// number of sieve tables (r) for SieveMethod
		private static int size;			// number of super digits in a multi-precision number
		private static int sizeO2;			// 1 + size / 2
		private static int xI;				// belongs to NextP()
		private static int[] kS;			// [rS+1]
		private static int[] squareMod;		// [mRM+1] for SieveMethod
		private static long stopDiv;		// for PrimeDiv function
		private static MultiSI cutoff;		// limit for using MethodX
		private static MultiSI delX;		// for SieveMethod
		private static MultiSI remS;		// for SieveMethod
		private static MultiSI root2;		// for SieveMethod
		private static MultiSI sieveN;		// for SieveMethod
		private static MultiSI start;		// minimum factor to start prime test, normally 27
		private static MultiSI z;			// factor of zIn not yet factored
		private static MultiSI zIn;			// original input number to be factored

		public static bool complete;		// flag that the input number was completely factored
		public static bool findUp;			// true to find first prime > zIn, false for first prime < zIn
		public static bool pTabFilled;		// pTab has been filled
		public static double log2OfBase;	// Math.Log(MultiID.baseM , 2)
		public static double squareMaxP;	// square of max prime in pTab
		public static int ib;				// universal index for a block transfer
		public static int lenPTab;			// 200 low, 3512 at max for Int16 primes
		public static int nF;				// number of prime factors found in zIn
		public static int PiLDepth;			// depth of calls to PrimePiL
		public static int sdBits;			// bits in maxDigit - 1
		public static int twoToSDB;			// 2 ^ sdBits = largest power of 2 < maxDigit
		public static int[] d;				// 12 items at 0 through 11, constant table used in MoreDiv array
		public static int[] pTab;			// small prime table
		public static int[] savedBits;		// bits saved by Len(z) function, sdBits per entry
		public static Tpowers powers;		// powers of the prime factors found in zIn
		public static TPrimes primes;		// prime factors found in zIn
	}
}

// Notes:
// when the command AllD(x), PFA(x), Sig(x), Sig0(x), Tau(x), or Phi(x) completes the
// factorization of x, this factorization is saved. When one of these commands is executed
// with x = 0 and there is a saved factorization, the saved factors and exponents are used
// to complete the command. This can save a lot of time when more than one of these
// commands are needed for the same number.
//
// executing one of these commands with x = 1 will delete the saved factors. The commands
// Prime(x), PGT(x), and PLT(x) will also delete any saved factors. The prime related
// commands PFB(x), Pi(x), P(x), and pTab(x) have no effect on the saved factors.
