using System;

namespace MyNamespace
{
	/// <summary>
	///		XIMult - Extra Precision Integer Multiply using Fast Hartley Transform Multiply
	///		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 XIMult
	{
		public XIMult()
		{
			// no constructor logic
		}

		public static void FHTnLdn(int xN, int yN, out int n, out int ldn) // compute n and ldn
		{
			int nMax = 2 * Math.Max(xN, yN);
			ldn = (int)(Math.Ceiling(Math.Log(nMax, 2)));
			n = 1 << ldn;

			if (n >= nMax + nMax)		// for safety, needed if nMax is a power of 2
			{							// and Log(n, 2) = integer + epsilon
				ldn--;					// the only one I could find was n = 2^29, log(n, 2) gave 29.000000000000004 (14 zeros)
				n = 1 << ldn;			// so this may not be needed but better safe ...
			}
		}

		private static int SumOfDigits(MultiUI x) // Sum of the digits of x mod (radix - 1)
		{
			int s = 0;
			for (int i = 0; i < x.N; i++)
			{
				s += x.V[i];
				if (s >= MultiID.maxDigit) s -= MultiID.maxDigit;
			}
			return s;
		}

		public static void XIMul(MultiUI x, MultiUI y, MultiUI z) // z = x * y
		// x or y may NOT have the same address in memory as z
		// x and y may have the same address, they are not changed
		// sets muErr true if overflow
		{
			MultiUI zL;
			zL = z;
			if (Object.ReferenceEquals(x, z) || Object.ReferenceEquals(y, z))
			{
				zL = new MultiUI(x.N + y.N);
			}
			int xN = x.N;
			int yN = y.N;
			int zN = xN + yN;
			if (zN > zL.M)
			{
				Common.ReportError("XIMul", "z.N " + zN + " > z.M " + zL.M);
				return;
			}
			if (Object.ReferenceEquals(x, y))
			{
				XISquare(x, zL);
				MultiID.SetToUI(zL, z);
				return;
			}
			int sodX;
			int sodY;
			int sodXY;
			int sodZ;
			double maxError;
			double radix;	// radix = max value of a super digit + 1
			int failed = -1;
			do
			{
				failed++;
				if (MultiID.muDMax != 8)
					radix = MultiID.baseM;		// radix = max value of a super digit + 1
				else if (failed < 2)
					radix = MultiID.baseMr;
				else
					radix = 100;
				GC.Collect();	// force garbage collection of all generations
				double timeAtStart = MultiID.TimerF();
				int n, ldn;
				double[] f;
				double[] g;
				if (MultiID.muDMax != 8)
					FHTnLdn(xN, yN, out n, out ldn);
				else if (failed < 2)
					FHTnLdn(2 * xN, 2 * yN, out n, out ldn);
				else
					FHTnLdn(4 * xN, 4 * xN, out n, out ldn);
				if (n > (1 << 25)) // 2^25 = 33,554,432
				{
					int dig = 8 * Math.Max(xN, yN); // convert n to decimal digits
					if (failed < 2)
					{
						MultiID.WriteLn("Numbers too large for FHT multiply, digits = ");
						MultiID.WriteLn(dig + " > 2^26 = 67,108,864");
					}
					else
					{
						MultiID.WriteLn("Numbers too large for radix 100 FHT multiply, digits = ");
						MultiID.WriteLn(dig + " > 2^25 = 33,554,432");
					}
					MultiID.muErr = true;
					MultiID.ClearUI(z);
					return;
				}
				sodX = SumOfDigits(x);
				sodY = SumOfDigits(y);
				sodXY = (int)(Math.BigMul(sodX, sodY) % MultiID.maxDigit);
				try
				{
					f = new double[n];		// initialized to zero, Out-of-Memory exception here
					g = new double[n];
				}
				catch (Exception ex)
				{
					long totalMemory = GC.GetTotalMemory(false);
					string st = ex.Message + ", n = " + n + ", z set to zero, Memory used = " + totalMemory;
					Common.ReportError("XIMul", st);
					z.N = 1;
					z.M = 1;
					z.V = new int[2];
					return;
				}
				if (MultiID.muDMax != 8)
				{
					for (int i = 0; i < xN; i++)	// copy x, int into double
						f[i] = x.V[i];
					for (int i = 0; i < yN; i++)	// copy y, int into double
						g[i] = y.V[i];
				}
				else if (failed < 2)
				{
					for (int i = 0; i < xN; i++)	// copy x, int into double
					{
						int rem;
						int q = Math.DivRem(x.V[i], MultiID.baseMr, out rem);
						f[i + i] = rem;
						f[i + i + 1] = q;
					}
					for (int i = 0; i < yN; i++)	// copy y, int into double
					{
						int rem;
						int q = Math.DivRem(y.V[i], MultiID.baseMr, out rem);
						g[i + i] = rem;
						g[i + i + 1] = q;
					}
				}
				else
				{
					for (int i = 0; i < xN; i++)  // copy x, int into double
					{
						int rem;
						int q = Math.DivRem(x.V[i], 100, out rem);
						f[4 * i] = rem;
						q = Math.DivRem(q, 100, out rem);
						f[4 * i + 1] = rem;
						q = Math.DivRem(q, 100, out rem);
						f[4 * i + 2] = rem;
						q = Math.DivRem(q, 100, out rem);
						f[4 * i + 3] = rem;
					}
					for (int i = 0; i < xN; i++)  // copy y, int into double
					{
						int rem;
						int q = Math.DivRem(y.V[i], 100, out rem);
						g[4 * i] = rem;
						q = Math.DivRem(q, 100, out rem);
						g[4 * i + 1] = rem;
						q = Math.DivRem(q, 100, out rem);
						g[4 * i + 2] = rem;
						q = Math.DivRem(q, 100, out rem);
						g[4 * i + 3] = rem;
					}
 				}
				FHTMult.FHTConvolution0(f, g, ldn);	// do Hartley convolution
				if (MultiID.muAbort)
				{
					MultiID.ClearUI(z);
					return;
				}
				FHTMult.radix = radix;
				FHTMult.ConvolutionRound(g, ldn, out maxError);
				FHTMult.ConvolutionCarry(g, ldn);
				zN = xN + yN;
				if (zN > zL.M)
				{
					zN = zL.M;						// overflow condition, should not happen
					MultiID.WriteLn("FHT overflow, should not happen but it did!!!!!!!");
				}
				if (MultiID.muDMax != 8)
				{
					if (g[zN - 1] == 0.0d) zN--;
					for (int i = 0; i < zN; i++)
						zL.V[i] = (int)g[i];
					zL.N = zN;
				}
				else if (failed < 2)
				{
					for (int i = 0, j = 0; i < zN; i++, j += 2)
						zL.V[i] = (int)g[j + 1] * MultiID.baseMr + (int)g[j];
					if (zL.V[zN - 1] == 0.0d) zN--;
					zL.N = zN;
				}
				else
				{
					for (int i = 0, j = 0; i < zN; i++, j += 4)
						zL.V[i] = (int)g[j + 3] * 1000000 + (int)g[j + 2] * 10000 + (int)g[j + 1] * 100 + (int)g[j];
					if (zL.V[zN - 1] == 0.0d) zN--;
					zL.N = zN;
				}
				sodZ = SumOfDigits(zL);
				bool test = false;
				if (test)
					sodZ = 0;
				if (sodZ != sodXY)
					MultiID.WriteLn("FHT Sum-of-Digits failed, error in FHT multiply!!!!!!!");
				if (Math.Abs(maxError) > 0.45)	//error, this should not happen either I hope
					MultiID.WriteLn("FHT Max Error = " + maxError + ", probable error in FHT multiply!!!!!!!");
			} while ((sodZ != sodXY || Math.Abs(maxError) > 0.45) && failed < 2);
			if (failed != 0)
				if (sodZ == sodXY && Math.Abs(maxError) <= 0.45)
					MultiID.WriteLn("FHT Recovered from failure/probable error in FHT multiply");
				else
				{
					MultiID.WriteLn("FHT cannot recover from failure/probable error in FHT multiply");
					MultiID.muErr = true;
				}
			MultiID.SetToUI(zL, z);
//if (MultiID.diagOn)
//{
//	double delta = MultiID.TimerF() - timeAtStart;
//	if (delta > 0.025)
//	{
//		string st = string.Format("FHT Radix = {0}  ldn = {1}  ", FHTMult.radix, ldn);
//		st += string.Format("Max Error = {0:0.00000}  DT = {1:0.00} sec.  ", maxError, delta);
//		st += string.Format("x.N = {0}  y.N = {1}  z.N = {2}", xN, yN, z.N);
//		MultiID.WriteTDLn(st);
//	}
//}
		} // XIMul

		public static void XISquare(MultiUI x, MultiUI z) // z = x * x
		// x may NOT have the same address in memory as z
		// x is not changed
		// sets muErr true if overflow
		{
			int sodX;
			int sodXY;
			int sodZ;
			double maxError;
			double radix;		// radix = max value of a super digit + 1
			int failed = -1;
			do
			{
				failed++;
				if (MultiID.muDMax != 8)
					radix = MultiID.baseM;      // radix = max value of a super digit + 1
				else if (failed < 2)
					radix = MultiID.baseMr;
				else
					radix = 100;
				GC.Collect();	// force garbage collection of all generations
				double timeAtStart = MultiID.TimerF();
				int xN = x.N;
				int n, ldn;
				double[] g;
				if (MultiID.muDMax != 8)
					FHTnLdn(xN, xN, out n, out ldn);
				else if (failed < 2)
					FHTnLdn(2 * xN, 2 * xN, out n, out ldn);
				else
					FHTnLdn(4 * xN, 4 * xN, out n, out ldn);
				if (n > (1 << 25)) // 2^25 = 33,554,432
				{
					int dig = 8*xN; // convert n to decimal digits
					if (failed < 2)
					{
						MultiID.WriteLn("Numbers too large for FHT multiply, digits = ");
						MultiID.WriteLn(dig + " > 2^26 = 67,108,864");
					}
					else
					{
						MultiID.WriteLn("Numbers too large for radix 100 FHT multiply, digits = ");
						MultiID.WriteLn(dig + " > 2^25 = 33,554,432");
					}
					MultiID.muErr = true;
					MultiID.ClearUI(z);
					return;
				}
				sodX = SumOfDigits(x);
				sodXY = (int)(Math.BigMul(sodX, sodX) % MultiID.maxDigit);
				try
				{
					g = new double[n];		// initialized to zero
				}
				catch (Exception ex)
				{
					long totalMemory = GC.GetTotalMemory(false);
					string st = ex.Message + ", n = " + n + ", z set to zero, Memory used = " + totalMemory;
					Common.ReportError("XIMul", st);
					z.N = 1;
					z.M = 1;
					z.V = new int[2];
					return;
				}
				if (MultiID.muDMax != 8)
				{
					for (int i = 0; i < xN; i++)	// copy x, int into double
						g[i] = x.V[i];
				}
				else if (failed < 2)
				{
					for (int i = 0; i < xN; i++)	// copy x, int into double
					{
						int rem;
						int q = Math.DivRem(x.V[i], MultiID.baseMr, out rem);
						g[i + i] = rem;
						g[i + i + 1] = q;
					}
				}
				else
				{
					for (int i = 0; i < xN; i++)	// copy x, int into double
					{
						int rem;
						int q = Math.DivRem(x.V[i], 100, out rem);
						g[4 * i] = rem;
						q = Math.DivRem(q, 100, out rem);
						g[4 * i + 1] = rem;
						q = Math.DivRem(q, 100, out rem);
						g[4 * i + 2] = rem;
						q = Math.DivRem(q, 100, out rem);
						g[4 * i + 3] = rem;
					}
				}
				FHTMult.FHTConvolution0(g, ldn);	// do Hartley convolution of g with g
				if (MultiID.muAbort)
				{
					MultiID.ClearUI(z);
					return;
				}
				FHTMult.radix = radix;
				FHTMult.ConvolutionRound(g, ldn, out maxError);
				FHTMult.ConvolutionCarry(g, ldn);
				int zN = xN + xN;
				if (zN > z.M)
				{
					zN = z.M;					// overflow condition, should not happen
					MultiID.WriteLn("FHT overflow, should not happen but it did!!!!!!!");
				}
				if (MultiID.muDMax != 8)
				{
					if (g[zN-1] == 0.0d) zN--;
					for (int i = 0; i < zN; i++)
						z.V[i] = (int)g[i];
					z.N = zN;
				}
				else if (failed < 2)
				{
					for (int i = 0, j = 0; i < zN; i++, j += 2)
						z.V[i] = (int)g[j + 1] * MultiID.baseMr + (int)g[j];
					if (z.V[zN - 1] == 0.0d) zN--;
					z.N = zN;
				}
				else
				{
					for (int i = 0, j = 0;  i < zN;  i++, j += 4)
						z.V[i] = (int)g[j + 3] * 1000000 + (int)g[j + 2] * 10000 + (int)g[j + 1] * 100 + (int)g[j];
					if (z.V[zN-1] == 0.0d) zN--;
					z.N = zN;
				}
				sodZ = SumOfDigits(z);
				bool test = false;
				if (test)
					sodZ = 0;
				if (sodZ != sodXY)
					MultiID.WriteLn("XISquare: FHT Sum-of-Digits failed, error in FHT multiply!!!!!!!");
				if (Math.Abs(maxError) > 0.45)	//error, this should not happen either I hope
					MultiID.WriteLn("XISquare: FHT Max Error = " + maxError + ", probable error in FHT multiply!!!!!!!");
			} while ((sodZ != sodXY || Math.Abs(maxError) > 0.45) && failed < 2);
			if (failed != 0)
				if (sodZ == sodXY && Math.Abs(maxError) <= 0.45)
					MultiID.WriteLn("FHT Recovered from failure/probable error in FHT multiply");
				else
				{
					MultiID.WriteLn("FHT cannot recover from failure/probable error in FHT multiply");
					MultiID.muErr = true;
				}
		} // XISquare
	}
}
