using System;
using System.IO; // for StreamReader, StreamWriter i/O to a file
using System.Windows.Forms;

namespace MyNamespace
{
	/// <summary>
	///		NPrimes - Algorithm 357 - Collected Algorithms from ACM
	///		An Efficient Prime Number Generator using the sieve of Eratosthenes
	///
	/// ---------------- *** 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 NPrimes
	{
		public NPrimes()				// constructor logic
		{
			NPrimeDone();				// set first time flag for function NPrime
			mCMaxS = 1000000;			// standard value for maximum mC
			maxT = 2750159;				// maximum integer testable. For reference only.
			maxT *= maxT;				// 756,33745,25281 (13 Digits) (Based on QS = 200000)
			jLimM = 0;
			first = true;				// first time flag for function NPrime
		}

		public NPrimes(int iJLim, int iQSize)	// constructor logic
		{
			first = true;				// first time flag for function NPrime
			jLim = iJLim + 1;			// size of working array, nominal value,
			jjLim = jLim;
			jLimM = iJLim;
			iP = new double[jLim+1];	// public working storage for primes
			qSize = iQSize;				// size of iQ and jQ arrays, must be large enough
										// to hold all primes <= Sqrt(Maximum integer testable)
			iQ = new double[qSize+1];	// secondary table of primes squared
			jQ = new int[qSize+1];		// secondary table of primes < Sqrt(current
										// prime) 2, 3, 5, 7, ... , 2750159 (if Qs = 200000)
		}

		public static void NPrimeDone() // end of prime number generator
		{
			iP = new double[0];			// re-dimention iP[0] to release memory
			iQ = new double[0];
			jQ = new int[0];
			jLim = 0;
			jjLim = 0.0;
			jLimM = 0;
			first = true;				// first time flag for function NPrime
		}

		public static void WriteSnapShotNP(StreamWriter writeSSFile)
		// write a snapshot file so calculations can be restarted after a program termination
		{
			writeSSFile.WriteLine("NPrime Function SnapShot: ij; ik; inc; j; nj; ni; njD");
			writeSSFile.WriteLine(" " + ij + " " + ik + " " + inc + " " + j + " " + nj + " " + ni + " " + njD);
			for (int i = 1; i <= ij; i++)
			{
				writeSSFile.WriteLine(i + " " + iQ[i] + " " + jQ[i]);
			}
		}

		public static bool ReadSnapShotNP(StreamReader readSSFile) // returns an ok flag
		// read the snapshot file so calculations can be restarted
		{
			string line = readSSFile.ReadLine();	// read in next line
			Global.PLn(line);

			line = readSSFile.ReadLine();			// read in next line
			line = Global.DeleteExtraDelimiters(line);
			string[] fields;
			fields	= line.Split(' ');				// convert the string to an array each, number is separated by a space
			ij		= Convert.ToInt32 ( fields[0] );
			ik		= Convert.ToInt32 ( fields[1] );
			inc		= Convert.ToInt32 ( fields[2] );
			j		= Convert.ToInt32 ( fields[3] );
			nj		= Convert.ToInt32 ( fields[4] );
			ni		= Convert.ToDouble( fields[5] );
			njD		= Convert.ToDouble( fields[6] );
			Global.PLn(" " + ij + " " + ik + " " + inc + " " + j + " " + nj + " " + ni + " " + njD);

			for (int i = 1; i <= ij; i++)
			{
				line	= readSSFile.ReadLine();		// read in next line
				fields	= line.Split(' ');				// convert the string to an array, each number is separated by a space
				int k	= Convert.ToInt32 ( fields[0] );
				iQ[i]	= Convert.ToDouble( fields[1] );
				jQ[i]	= Convert.ToInt32 ( fields[2] );

				if (i == 1 || i == ij) Global.PLn(k + " " + iQ[i] + " " + jQ[i]);
				if (i == 1) Global.PLn(" . . .");
			}
			iP[jLim] = 0;
			if (j != jLimM)
			{
				Global.PLn("");
				Global.PLn("   <><><><><> SnapShot File Has Wrong j Value <><><><><>");
				Global.PLn("SnapShot j = " + j + ", Should be " + (jLimM));
				Global.PLn("");
				return false;
			}
			return true;
		}

		public static int NPrime(int mC) // generate mC more primes in iP[1] thru iP[mC]
		{
			return NPrime(mC, false);
		}

		public static int NPrime(int mC, bool flagOnlyATest) // generate mC more primes in iP[1] thru iP[mC]
		{
			// Algorithm 357 - Collected Algorithms from ACM
			// An Efficient Prime Number Generator using the sieve of Eratosthenes

			// mC = Number of new primes desired this call, must be >= 2 on 1st call
			// ij, the current length of the tables in array iQ and jQ is returned
			//   This is approximately the number of primes < Sqrt(current prime)
			//   It must not exceed qSize !!!
			// NB: if SnapShot files are used, mC must be iJLim = jLim - 1
			//   and iJLim cannot be changed between write and read.
			// algorithm 357 was originally written in ALGOL and FORTRAN with indices
			// running from 1 to n. I left it that way.

			int i;		// array index
			int jqi;	// jQ[i] to store in iP[j]
			int k;		// array index

			if (flagOnlyATest)		// test "Too Large"
			{
				iQ[qSize + 1] = 0;	// cause an index out-of-bounds exception
			}
			k = 0;
			if (!first) goto Lf;
			if (mC < 2) return 0;	// mC must be >= 2 on first call
			// first time entered so set initial conditions
			first = false;
			for (i = 1; i <= qSize; i++)
			{
				iQ[i] = 0;
				jQ[i] = 0;
			}
			iP[1] = 2.0;
			jQ[1] = ik = inc = 2;
			jQ[2] = ij = 3;
			jQ[3] = nj = 5;
			iQ[1] = 3.0;
			iQ[2] = 9.0;
			iQ[3] = 25.0;
			k = 1;
			// prepare to delete a sequence of composite numbers
		La:
			j = k + 1;
			ni = iQ[1] - j - j;
			iQ[1] = ni + jjLim + jjLim;
			for (i = j; i <= jLim; i++) iP[i] = 0.0;
			// diagnostic
			// Global.PLn("iQ[1]= " + iQ[1] + " nI= " + nI);
		Lb:
			i = ij;
			if (iQ[ij] >= iQ[1]) goto Le;
			// extend the list of primes in array jQ counting so as to omit
			// multiples of 2 and 3
		Lc:
			nj += inc;
			inc = 6 - inc;
			j = jQ[ik + 1];
			if (j * j <= nj) ik++;
			for (j = 3; j <= ik; j++)
			{
				if (nj % jQ[j] == 0) goto Lc;
			}
			ij++;
			jQ[ij] = nj;
			njD = nj;
			iQ[ij] = njD * njD;
			// Global.PLn("iQ[" + ij + "]=" + iQ[ij]); // diagnostic
			goto Lb;
			// if j + j + ni is composite, enter its smallest prime factor in iP[j].
			// if j + j + ni is prime, then set iP[j] = 0
		Ld:
			iP[j] = jqi;
			j += jqi;		// put prime in all slots divisible by it
			if (j < jLim) goto Ld;
			iQ[i] = (ni + j) + j;
			// diagnostic
			// Global.PLn "iQ[" + i + "]=" + iQ(i) + " j=" + j + " nI=" + nI
		Le:
			i--;
			jj = Math.Floor((iQ[i] - ni) / 2.0);
			if (jj < jjLim) { j = (int)jj; jqi = jQ[i]; goto Ld; }
			if (i != 1) goto Le;
			j = k;
			// pack the next mC primes in iP[1], ..., iP[mC]
		Lf:
			j++;
			if (iP[j] != 0) goto Lf;
			if (j == jLim) goto La;
			k++;
			iP[k] = (ni + j) + j;
			if (k % 5000 == 0) Application.DoEvents();
			if (k != mC) goto Lf;
			return ij; // return the length of the tables in arrays iQ and jQ
		} // NPrime

		//static local variables for function NPrime
		private static double jj;	// j a double
		private static double maxT;	// maximum integer testable. For reference only, not used
		private static double ni;	// test number - 2J for primes
		private static double njD;	// nj in DOUBLE format
		private static int ij;		// current length of the tables in array iQ and jQ
		private static int ik;		// max index into jQ to test nj
		private static int inc;		// = 2, 4, 2, 4, ...
		private static int j;		// array index
		private static int nj;		// test number for secondary table of primes
		private static int qSize;	// size of iQ and jQ arrays, must be large enough
									//   to hold all primes <= Sqrt(Maximum integer testable)
		public static bool first;	// first time flag for function NPrime (true = first)
		public static double jjLim;	// jLim as a double
		public static double[] iP;	// public working storage for primes
		public static double[] iQ;	// secondary table of primes squared
		public static int jLim;		// size of working array, nominal value,
									// don't ask for more than this (mC < jLim)
		public static int jLimM;	// jLim - 1
		public static int mCMaxS;	// standard value for maximum mC
		public static int[] jQ;		// secondary table of primes < Sqrt(current
	}
}
