#define	RBP0_SYMM2		// 2, 4 (default is 2)

#define	RBP0_SYMM4		// 2, 4 (default is 2)
#undef	RBP0_SYMM4		// 2, 4 (default is 2)

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

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

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

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

using System;
using System.IO;		// for TextWriter, StreamWriter I/O to a file
using System.Text;		// for Encoding I/O to a Text file

namespace MyNamespace
{
	/// <summary>
	///		FHTMult - Fast Hartley Transform Multiply
	///		fxt subroutines converted from C++ to C#
	///		including FHT Convolution with zero padded date
	///		Author: Harry J. Smith, Private individual.
	///		Copyright (c) 1981-2005 by author: Harry J. Smith,
	///		19628 Via Monte Dr., Saratoga CA 95070.  All rights reserved.
	///
	///		C++ author = Joerg Arndt		email: arndt@jjj.de
	///		the C++ software is online at	http://www.jjj.de/
	///			email: hjsmithh@sbcglobal.net
	///			web page: http://www.geocities.com/hjsmithh/
	///
	/// ---------------- *** 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 FHTMult
	{
		public FHTMult()	// constructor logic
		{
#if  TEST
#else
			if (MultiID.muDMax != 8)
				radix = MultiID.baseM;			// radix = max value of a super digit + 1
			else
				radix = MultiID.baseMr;
#endif // TEST
			sqrt2 = Math.Sqrt(2);				// 2^0.5
			sqrt1O2 = Math.Sqrt(0.5);			// 0.5^0.5 M_SQRT1_2
			piO4 = Math.PI / 4.0;				// pi Over 4
			cosPiO8 = Math.Cos(Math.PI / 8.0);	// cos(pi/8)
			sinPiO8 = Math.Sin(Math.PI / 8.0);	// sin(pi/8)
			bitsPerInt = 32;					// bits in an type int integer
		}

		public static void Swap0(ref double x, ref double y)
		// swap() for y known to be zero
		{
			y = x;
			x = 0.0;
		}

		public static void Swap(ref double x, ref double y) // using the name Swap for both double Complex
		{
			double temp = y;
			y = x;
			x = temp;
		}

// the following has been replaced with inline code for speedup
//		public static void IdxSwap(double[] f, int k, int r)
//		// was #define  idx_swap(f, k, r)  { ulong kx=(k), rx=(r);  Swap(f[kx], f[rx]); }
//		// indexed swap
//		{
//			double temp = f[r];
//			f[r] = f[k];
//			f[k] = temp;
//		}
//
//		public static void CSqr(double a, double b, out double u, out double v)
//		// was csqr(double a, double b, double &u, double &v);
//		// {u,v} <--| {a*a-b*b, 2*a*b}
//		{
//			u = a*a - b*b;
//			v = 2*a*b;
//		}
//
//		public static void CMult(double c, double s, double x, double y, out double u, out double v)
//		// was cmult(double c, double s, double x, double y, double &u, double &v);
//		// {u,v} <--| {x*c-y*s, x*s+y*c}
//		{
//			u = x*c - y*s;
//			v = x*s + y*c;
//		}
//
//		public static void SumDiff(ref double a, ref double b)
//		{
//			double t = a;
//			a = t + b;
//			b = t - b;
//		}
//
//		public static void SumDiff(double a, double b, out double s, out double d)
//		{
//			s = a + b;
//			d = a - b;
//		}
//
//		public static void SinCos(double a, out double s, out double c)
//		{
//			s = Math.Sin(a);
//			c = Math.Cos(a);
//		}
//
//		public static void FHTMul(double xi, double xj, ref double yi, ref double yj, double v)
//		// was static inline void fht_mul(double xi, double xj, double &yi, double &yj, double v)
//		// Auxiliary routine for fht_convolution_core()
//		{
//			double h1p = xi,  h1m = xj;
//			double s1 = h1p + h1m,  d1 = h1p - h1m;
//			double h2p = yi,  h2m = yj;
//			yi = (h2p * s1 + h2m * d1) * v;
//			yj = (h2m * s1 - h2p * d1) * v;
//		}

		public static void RevbinPermute(double[] f, int n)
		{
			double temp;
			if (n <= 8)
			{
				if (n == 8)
				{
					Swap(ref f[1], ref f[4]);
					Swap(ref f[3], ref f[6]);
				}
				else if (n == 4)  Swap(ref f[1], ref f[2]);

				return;
			}

			int nh = n >> 1;
			int[] x = new int[bitsPerInt];	// was ALLOCA(ulong, x, BITS_PER_LONG);
			x[0] = nh;
			{	// initialize xor-table:
				int i, m = nh;
				for (i = 1;  m != 0;  ++i)
				{
					m >>= 1;
					x[i] = x[i-1] ^ m;
				}
			}
#if  ( RBP0_SYMM2 || RBP0_SYMM4 )
			int rbp0Symm = 2;
			int n1 = n - 1;		// = 11111111
#if  ( RBP0_SYMM4 )
			rbp0Symm = 4;
			int nx1 = nh - 2;		// = 01111110
			int nx2 = n1 - nx1;	// = 10111101
#endif //  ( RBP0_SYMM4 )
#endif //  ( RBP0_SYMM2 || RBP0_SYMM4 )
			int k = 0, r = 0;
			int kk, rr;
			while (k < n/rbp0Symm)	// n>=16, n/2>=8, n/4>=4
			{
				// ----- k%4 == 0:
				if (r > k)
				{
//					Swap(ref f[k], ref f[r]);		// <nh, <nh 11
					temp = f[r];	f[r] = f[k];	f[k] = temp;
#if  ( RBP0_SYMM2 || RBP0_SYMM4 )
//					IdxSwap(f, n1^k, n1^r);		// >nh, >nh 00
					temp = f[rr = n1^r];	f[rr] = f[kk = n1^k];	f[kk] = temp;
#if  ( RBP0_SYMM4 )
					IdxSwap(f, nx1^k, nx1^r);	// <nh, <nh 11
					IdxSwap(f, nx2^k, nx2^r);	// >nh, >nh 00
#endif //  ( RBP0_SYMM4 )
#endif //  ( RBP0_SYMM2 || RBP0_SYMM42 )
				}

				r ^= nh;
				++k;

				// ----- k%4 == 1:
				if (r > k)
				{
//					Swap(ref f[k], ref f[r]);		// <nh, >nh 10
					temp = f[r];	f[r] = f[k];	f[k] = temp;
#if  ( RBP0_SYMM4 )
					IdxSwap(f, n1^k, n1^r);		// >nh, <nh 01
#endif //  ( RBP0_SYMM4 )
				}

				{ // scan for lowest unset bit of k:
#if  RBP_USE_ASM
					int i = LowestBitIdx(~k);
#else
					int m = 2,  i = 1;
					while ((m & k) != 0)  { m <<= 1;  ++i; }
#endif // RBP_USE_ASM
					r ^= x[i];
				}
				++k;

				// ----- k%4 == 2:
				if (r > k)
				{
//					Swap(ref f[k], ref f[r]);	// <nh, <nh 11
					temp = f[r];	f[r] = f[k];	f[k] = temp;
#if  ( RBP0_SYMM2 )
//					IdxSwap(f, n1^k, n1^r); // >nh, >nh 00
					temp = f[rr = n1^r];	f[rr] = f[kk = n1^k];	f[kk] = temp;
#endif //  ( RBP0_SYMM2 )
				}

				r ^= nh;
				++k;

				// ----- k%4 == 3:
				if (r > k)
				{
//					Swap(ref f[k], ref f[r]);	// <nh, >nh 10
					temp = f[r];	f[r] = f[k];	f[k] = temp;
#if  ( RBP0_SYMM4 )
					IdxSwap(f, nx1^k, nx1^r);	// <nh, >nh 10
#endif //  ( RBP0_SYMM4 )
				}

				{ // scan for lowest unset bit of k:
#if  RBP_USE_ASM
					int i = LowestBitIdx(~k);
#else
					int m = 4,  i = 2;
					while ((m & k) != 0)  { m <<= 1;  ++i; }
#endif // RBP_USE_ASM
					r ^= x[i];
				}
				++k;
			}
		} // RevbinPermute(double...

		public static void RevbinPermute0(double[] f, int n)
		// was void revbin_permute0(Type *f, ulong n)
		// put data in revbin order
		// version for zero padded data
		{
			if (n <= 8)
			{
				if (n == 8)
				{
					Swap0(ref f[1], ref f[4]);
					Swap0(ref f[3], ref f[6]);
				}
				else if (n == 4)  Swap0(ref f[1], ref f[2]);

				return;
			}

			int nh = n >> 1;
			int[] x = new int[bitsPerInt];	// was ALLOCA(ulong, x, BITS_PER_LONG);
			x[0] = nh;
			{	// initialize xor-table:
				int i, m = nh;
				for (i = 1;  m != 0;  ++i)
				{
					m >>= 1;
					x[i] = x[i-1] ^ m;
				}
			}
			int rbp0Symm = 2;
#if  ( RBP0_SYMM2 || RBP0_SYMM4 )
#if  ( RBP0_SYMM4 )
			rbp0Symm = 4;
			int n1  = n - 1;		// = 11111111
			int nx1 = nh - 2;		// = 01111110
			// int nx2 = n1 - nx1;	// = 10111101
#endif //  ( RBP0_SYMM4 )
#endif //  ( RBP0_SYMM2 || RBP0_SYMM4 )
			int k = 0, r = 0;
			double temp;
			while (k < n/rbp0Symm)	// n>=16, n/2>=8, n/4>=4
			{
				// ----- k%4 == 0:
				if (r > k)
				{
//					Swap(ref f[k], ref f[r]);	// <nh, <nh 11
					temp = f[r];	f[r] = f[k];	f[k] = temp;
#if  ( RBP0_SYMM2 || RBP0_SYMM4 )
					// IdxSwap(f, n1^k, n1^r);	// >nh, >nh 00
#if  ( RBP0_SYMM4 )
					IdxSwap(f, nx1^k, nx1^r);		// <nh, <nh 11
					// IdxSwap(f, nx2^k, nx2^r);	// >nh, >nh 00
#endif //  ( RBP0_SYMM4 )
#endif //  ( RBP0_SYMM2 || RBP0_SYMM42 )
				}

				r ^= nh;
				++k;

				// ----- k%4 == 1:
				if (r > k)
				{
//					Swap0(ref f[k], ref f[r]);	// <nh, >nh 10
					f[r] = f[k];	f[k]= 0.0d;
#if  ( RBP0_SYMM4 )
					Swap0(ref f[r^n1], ref f[k^n1]);
					// IdxSwap(f, n1^k, n1^r);	// >nh, <nh 01
#endif //  ( RBP0_SYMM4 )
				}

				{ // scan for lowest unset bit of k:
					int m = 2,  i = 1;
					while ((m & k) != 0)  { m <<= 1;  ++i; }
					r ^= x[i];
				}
				++k;

				// ----- k%4 == 2:
				if (r > k)
				{
//					Swap(ref f[k], ref f[r]);	// <nh, <nh 11
					temp = f[r];	f[r] = f[k];	f[k] = temp;
#if  ( RBP0_SYMM2 )
					// IdxSwap(f, n1^k, n1^r); // >nh, >nh 00
#endif //  ( RBP0_SYMM2 )
				}

				r ^= nh;
				++k;

				// ----- k%4 == 3:
				if (r > k)
				{
//					Swap0(ref f[k], ref f[r]);	// <nh, >nh 10
					f[r] = f[k];	f[k]= 0.0d;
#if  ( RBP0_SYMM4 )
					Swap0(ref f[k^nx1], ref f[r^nx1]);	// <nh, >nh 10
					// IdxSwap(f, nx1^k, nx1^r);	// <nh, >nh 10
#endif //  ( RBP0_SYMM4 )
				}

				{ // scan for lowest unset bit of k:
					int m = 4,  i = 2;
					while ((m & k) != 0)  { m <<= 1;  ++i; }
					r ^= x[i];
				}
				++k;
			}
		} // RevbinPermute0(double...

		public static void FHT(double[] f, int ldn)
		//void fht_dit2(double *f, ulong ldn)
		// decimation in time radix-2 fht
		{
			int n = 1 << ldn;
			RevbinPermute(f, n);
			for (int ldm = 1; ldm <= ldn; ++ldm)
			{
				int m = 1 << ldm;
				int mh = m >> 1;
				int m4 = mh >> 1;
				double phi0 = Math.PI / mh;

				for (int r = 0; r < n; r += m)
				{
					{	// j == 0:
						int t1 = r;
						int t2 = t1 + mh;
//						SumDiff(ref f[t1], ref f[t2]);
						double ta = f[t1], tb = f[t2];	f[t1] = ta + tb;	f[t2] = ta - tb;
					}

					if (m4 != 0)
					{
						int t1 = r + m4;
						int t2 = t1 + mh;
//						SumDiff(ref f[t1], ref f[t2]);
						double ta = f[t1], tb = f[t2];	f[t1] = ta + tb;	f[t2] = ta - tb;
					}
				}

				for (int j = 1, k = mh-1;  j < k;  ++j, --k)
				{
					double s, c, a;
//					SinCos(phi0*j, out s, out c);
					s = Math.Sin(a = phi0*j);	c = Math.Cos(a);

					for (int r = 0; r < n; r += m)
					{
						int tj = r + mh + j;
						int tk = r + mh + k;
						double fj = f[tj];
						double fk = f[tk];
						f[tj] = fj * c + fk * s;
						f[tk] = fj * s - fk * c;

						int t1 = r + j;
						int t2 = tj;	// == t1 + mh;
//						SumDiff(ref f[t1], ref f[t2]);
						double ta = f[t1], tb = f[t2];	f[t1] = ta + tb;	f[t2] = ta - tb;

						t1 = r + k;
						t2 = tk;		// == t1 + mh;
//						SumDiff(ref f[t1], ref f[t2]);
						ta = f[t1];	tb = f[t2];	f[t1] = ta + tb;	f[t2] = ta - tb;
					}
#if  TEST
#else
					MultiID.MuInterrupt();
					if (MultiID.muAbort) return;
					if (MultiID.statReq)
					{
						MultiID.Write("FHT: " + Math.Floor(100.0 * ldm / ldn).ToString() + " Percent done");
						MultiID.TimeDel();
						RunForm.FocusBothUs();
						MultiID.statReq = false;
					}
#endif // TEST
				}
			}
		} // FHT(double...

		public static void FHT0(double[] f, int ldn)
		// was void fht0(double *f, ulong ldn)
		// Fast Hartley Transform
		// version for zero padded data:
		//		f[k] == 0.0 for k=n/2 ... n-1
		// ldn := base-2 logarithm of the array length
		// split radix decimation in time algorithm
		{
			if (ldn <= 2)
			{
				if (ldn == 1)		// two point fht
				{
					f[1] = f[0];
				}
				else if (ldn == 2)	// four point fht
				{
					double f0, f2;
//					SumDiff(f[0], f[1], out f0, out f2);
					double ta = f[0], tb = f[1];	f0 = ta + tb;	f2 = ta - tb;

					f[0] = f[1] = f0;
					f[2] = f[3] = f2;
				}
				return;
			}
			int n =  1 << ldn;
			RevbinPermute0(f, n);

			int ldk = ldn & 1;
			if (ldk == 0)		// ldn is multiple of 2  => n is a power of 4
			{
#if  ( INITIAL_RADIX_16 )
				for (int fi = 0;  fi < n;  fi += 16) // radix-16 step
				{
//					SumDiff(ref f[0+fi], ref f[2+fi]);
					double ta = f[fi], tb = f[2+fi];	f[fi] = ta + tb;	f[2+fi] = ta - tb;

					f[3+fi] = f[2+fi];
					f[1+fi] = f[0+fi];
//					SumDiff(ref f[4+fi], ref f[6+fi]);
					ta = f[4+fi]; tb = f[6+fi];	f[4+fi] = ta + tb;	f[6+fi] = ta - tb;

					f[7+fi] = f[6+fi];
					f[5+fi] = f[4+fi];
//					SumDiff(ref f[8+fi], ref f[10+fi]);
					ta = f[8+fi]; tb = f[10+fi];	f[8+fi] = ta + tb;	f[10+fi] = ta - tb;

					f[11+fi] = f[10+fi];
					f[9+fi] = f[8+fi];
//					SumDiff(ref f[12+fi], ref f[14+fi]);
					ta = f[12+fi]; tb = f[14+fi];	f[12+fi] = ta + tb;	f[14+fi] = ta - tb;

					f[15+fi] = f[14+fi];
					f[13+fi] = f[12+fi];

					double f0, f1, f2, f3;
//					SumDiff(f[0+fi], f[4+fi], out f0, out f1);
					ta = f[fi];	tb = f[4+fi];	f0 = ta + tb;	f1 = ta - tb;

//					SumDiff(f[8+fi], f[12+fi], out f2, out f3);
					ta = f[8+fi];	tb = f[12+fi];	f2 = ta + tb;	f3 = ta - tb;

//					SumDiff(f0, f2, out f[0+fi], out f[8+fi]);
					f[fi] = f0 + f2;	f[8+fi] = f0 - f2;

//					SumDiff(f1, f3, out f[4+fi], out f[12+fi]);
					f[4+fi] = f1 + f3;	f[12+fi] = f1 - f3;

//					SumDiff(f[2+fi], f[6+fi], out f0, out f1);
					ta = f[2+fi];	tb = f[6+fi];	f0 = ta + tb;	f1 = ta - tb;

					f3 = sqrt2 * f[14+fi];
					f2 = sqrt2 * f[10+fi];
//					SumDiff(f0, f2, out f[2+fi], out f[10+fi]);
					f[2+fi] = f0 + f2;	f[10+fi] = f0 - f2;

//					SumDiff(f1, f3, out f[6+fi], out f[14+fi]);
					f[6+fi] = f1 + f3;	f[14+fi] = f1 - f3;

					double a, b, g0, g1, g2, g3;
//					SumDiff(f[5+fi], f[7+fi], out a, out b);
					ta = f[5+fi];	tb = f[7+fi];	a = ta + tb;	b = ta - tb;

					a *= sqrt1O2;
					b *= sqrt1O2;
//					SumDiff(f[1+fi], a, out f0, out f1);
					ta = f[1+fi];	f0 = ta + a;	f1 = ta - a;

//					SumDiff(f[3+fi], b, out g0, out g1);
					ta = f[3+fi];	g0 = ta + b;	g1 = ta - b;

//					SumDiff(f[13+fi], f[15+fi], out a, out b);
					ta = f[13+fi];	tb = f[15+fi];	a = ta + tb;	b = ta - tb;

					a *= sqrt1O2;
					b *= sqrt1O2;
//					SumDiff(f[9+fi], a, out f2, out f3);
					ta = f[9+fi];	f2 = ta + a;	f3 = ta - a;

//					SumDiff(f[11+fi], b, out g2, out g3);
					ta = f[11+fi];	g2 = ta + b;	g3 = ta - b;

					double c1 = cosPiO8;	// jjkeep
					double s1 = sinPiO8;	// jjkeep
//					CMult(s1, c1, f2, g3, out b, out a);
					b = f2*s1 - g3*c1;		a = f2*c1 + g3*s1;

//					SumDiff(f0, a, out f[1+fi], out f[9+fi]);
					f[1+fi] = f0 + a;	f[9+fi] = f0 - a;

//					SumDiff(g1, b, out f[7+fi], out f[15+fi]);
					f[7+fi] = g1 + b;	f[15+fi] = g1 - b;

//					CMult(c1, s1, g2, f3, out b, out a);
					b = g2*c1 - f3*s1;		a = g2*s1 + f3*c1;

//					SumDiff(g0, a, out f[3+fi], out f[11+fi]);
					f[3+fi] = g0 + a;	f[11+fi] = g0 - a;

//					SumDiff(f1, b, out f[5+fi], out f[13+fi]);
					f[5+fi] = f1 + b;	f[13+fi] = f1 - b;
				}
				ldk = 4;
#else // INITIAL_RADIX_16
//				for (double *fi=f; fi<fn; fi+=4)	// radix-4 step
				for (int fi = 0;  fi < n;  fi += 4)	// radix-4 step
				{
					SumDiff(ref f[0+fi], ref f[2+fi]);
					f[3+fi] = f[2+fi];
					f[1+fi] = f[0+fi];
				}
				ldk = 2;
#endif // INITIAL_RADIX_16
			}
			else		// ldk==1,  ldn is no multiple of 2  => n is no power of 4
			{
				for (int fi = 0;  fi < n;  fi += 8) // radix-8 step
				{
					double f0, f1;
//					SumDiff(f[0+fi], f[2+fi], out f0, out f1);
					double ta = f[fi], tb = f[2+fi];	f0 = ta + tb;	f1 = ta - tb;

					double s1 = f[4+fi] * sqrt2;
					double s2 = f[6+fi] * sqrt2;
					double f2, f3;
//					SumDiff(f[4+fi], f[6+fi], out f2, out f3);
					ta = f[4+fi];	tb = f[6+fi];	f2 = ta + tb;	f3 = ta - tb;

//					SumDiff(f0, f2, out f[0+fi], out f[4+fi]);
					f[fi] = f0 + f2;	f[4+fi] = f0 - f2;

//					SumDiff(f1, f3, out f[2+fi], out f[6+fi]);
					f[2+fi] = f1 + f3;	f[6+fi] = f1 - f3;

//					SumDiff(f0, s1, out f[1+fi], out f[5+fi]);
					f[1+fi] = f0 + s1;	f[5+fi] = f0 - s1;

//					SumDiff(f1, s2, out f[3+fi], out f[7+fi]);
					f[3+fi] = f1 + s2;	f[7+fi] = f1 - s2;
				}
				ldk = 3;
			}

			for ( ; ldk < ldn; ldk += 2)
			{
				int k	= 1 << ldk;
				int kh  = k >> 1;
				int k2  = k << 1;
				int k3  = k2 + k;
				int k4  = k2 << 1;

				for (int fi = 0, gi = fi + kh;  fi < n;  fi += k4, gi += k4)
				{
					double f0, f1, f2, f3;
//					SumDiff(f[0+fi], f[k+fi], out f0, out f1);
					double ta = f[fi],	tb = f[k+fi];	f0 = ta + tb;	f1 = ta - tb;

//					SumDiff(f[k2+fi], f[k3+fi], out f2, out f3);
					ta = f[k2+fi];	tb = f[k3+fi];	f2 = ta + tb;	f3 = ta - tb;

//					SumDiff(f0, f2, out f[0+fi], out f[k2+fi]);
					f[fi] = f0 + f2;	f[k2+fi] = f0 - f2;

//					SumDiff(f1, f3, out f[k+fi], out f[k3+fi]);
					f[k+fi] = f1 + f3;	f[k3+fi] = f1 - f3;

//					SumDiff(f[0+gi], f[k+gi], out f0, out f1);
					ta = f[gi];	tb = f[k+gi];	f0 = ta + tb;	f1 = ta - tb;

					f3 = sqrt2 * f[k3+gi];
					f2 = sqrt2 * f[k2+gi];
//					SumDiff(f0, f2, out f[0+gi], out f[k2+gi]);
					f[gi] = f0 + f2;	f[k2+gi] = f0 - f2;

//					SumDiff(f1, f3, out f[k+gi], out f[k3+gi]);
					f[k+gi] = f1 + f3;	f[k3+gi] = f1 - f3;
				}

				double tt = (piO4)/kh;	// was M_PI_4/kh;	// jjkeep
#if	TRIG_REC
				double s1 = 0.0,  c1 = 1.0;	// jjkeep
				double al = Math.Sin(0.5*tt);	// jjkeep
				al *= (2.0*al);
				double be = Math.Sin(tt);	// jjkeep
#endif // TRIG_REC
				for (int i = 1;  i < kh;  ++i)
				{
#if	TRIG_REC
					c1 -= (al*(tt=c1)+be*s1);
					s1 -= (al*s1-be*tt);
#else
					double s1, c1, aa;	// jjkeep
//					SinCos(tt*i, out s1, out c1);
					s1 = Math.Sin(aa = tt*i);	c1 = Math.Cos(aa);
#endif // TRIG_REC

					double c2, s2;	// jjkeep
//					CSqr(c1, s1, out c2, out s2);
					c2 = c1*c1 - s1*s1;		s2 = 2*c1*s1;

					for (int fi = i, gi = k-i;  fi < n;  fi += k4, gi += k4)
					{
						double a, b, g0, f0, f1, g1, f2, g2, f3, g3;
						int ff, gg;
//						CMult(s2, c2, f[k+fi], f[k+gi], out b, out a);
						b = f[ff = k+fi]*s2 - f[gg = k+gi]*c2;	a = f[ff]*c2 + f[gg]*s2;

//						SumDiff(f[0+fi], a, out f0, out f1);
						double ta = f[fi];	f0 = ta + a;	f1 = ta - a;

//						SumDiff(f[0+gi], b, out g0, out g1);
						ta = f[gi];	g0 = ta + b;	g1 = ta - b;

//						CMult(s2, c2, f[k3+fi], f[k3+gi], out b, out a);
						b = f[ff = k3+fi]*s2 - f[gg = k3+gi]*c2;	a = f[ff]*c2 + f[gg]*s2;

//						SumDiff(f[k2+fi], a, out f2, out f3);
						ta = f[k2+fi];	f2 = ta + a;	f3 = ta - a;

//						SumDiff(f[k2+gi], b, out g2, out g3);
						ta = f[k2+gi];	g2 = ta + b;	g3 = ta - b;

//						CMult(s1, c1, f2, g3, out b, out a);
						b = f2*s1 - g3*c1;	a = f2*c1 + g3*s1;

//						SumDiff(f0, a, out f[0+fi], out f[k2+fi]);
						f[fi] = f0 + a;	f[k2+fi] = f0 - a;

//						SumDiff(g1, b, out f[k+gi], out f[k3+gi]);
						f[k+gi] = g1 + b;	f[k3+gi] = g1 - b;

//						CMult(c1, s1, g2, f3, out b, out a);
						b = g2*c1 - f3*s1;	a = g2*s1 + f3*c1;

//						SumDiff(g0, a, out f[0+gi], out f[k2+gi]);
						f[gi] = g0 + a;	f[k2+gi] = g0 - a;

//						SumDiff(f1, b, out f[k+fi], out f[k3+fi]);
						f[k+fi] = f1 + b;	f[k3+fi] = f1 - b;
					}
#if  TEST
#else
					MultiID.MuInterrupt();
					if (MultiID.muAbort) return;
					if (MultiID.statReq)
					{
						MultiID.Write("FHT0: " + Math.Floor(100.0 * ldk / ldn).ToString() + " Percent done");
						MultiID.TimeDel();
						RunForm.FocusBothUs();
						MultiID.statReq = false;
					}
#endif // TEST
				}
			}
		} // FHT0(double

		public static void FHTConvolution0(double[] f, double[] g, int ldn)
		// was void	fht_convolution0(double * restrict f, double * restrict g, ulong ldn)
		// (linear, real) convolution:  g[] :=  f[] (*) g[]
		// ldn := base-2 logarithm of the array length
		// input data must be zero padded:
		//		f[n/2] .. f[n-1] == 0 and g[n/2] .. g[n-1] == 0
		// n = 2**ldn  must be >=2
		{
#if  TEST
#else
			MultiID.statReq = false;
#endif // TEST
			FHT0(f, ldn);
#if  TEST
#else
			if (MultiID.muAbort) return;
#endif // Test

			FHT0(g, ldn);
#if  TEST
#else
			if (MultiID.muAbort) return;
#endif // TEST

			FHTConvolutionCore(f, g, ldn, 0.0);

			FHT(g, ldn);
		}

		public static void FHTConvolution0(double[] g, int ldn)
		// (linear, real) convolution:  g[] :=  g[] (*) g[]
		// ldn := base-2 logarithm of the array length
		// input data must be zero padded:
		//		g[n/2] .. g[n-1] == 0
		// n = 2**ldn  must be >=2
		{
#if  TEST
#else
			MultiID.statReq = false;
#endif // TEST
			FHT0(g, ldn);
#if  TEST
#else
			if (MultiID.muAbort) return;
#endif // TEST

			FHTConvolutionCore(g, g, ldn, 0.0);

			FHT(g, ldn);
		}

		public static void FHTConvolutionCore(double[] f, double[] g, int ldn, double v)
		// was void	fht_convolution_core(const double * restrict f, double * restrict g, ulong ldn,
		//				double v/*=0.0*/) // jjkeep
		// Auxiliary routine for the computation of convolutions
		// via Fast Hartley Transforms.
		// ldn := base-2 logarithm of the array length.
		// v!=0.0 chooses alternative normalization.
		{
			int n = 1 << ldn;
			if (v == 0.0) v = 1.0 / n;
			g[0] *= (v * f[0]);
			int nh = n/2;

			if (nh > 0)
			{
				g[nh] *= (v * f[nh]);
				v *= 0.5;
				for (int i = 1, j = n-1; i < j; ++i, j--)
				{
//					FHTMul(f[i], f[j], ref g[i], ref g[j], v);
					double h1p = f[i],  h1m = f[j];
					double s1 = h1p + h1m,  d1 = h1p - h1m;
					double h2p = g[i],  h2m = g[j];
					g[i] = (h2p * s1 + h2m * d1) * v;
					g[j] = (h2m * s1 - h2p * d1) * v;
				}
			}
		} // FHTConvolutionCore(double...

		public static void ConvolutionRound(double[] g, int ldn, out double maxError)
		{
			int n = 1 << ldn;
			maxError = 0.0;
			for (int i = 0; i < n; ++i)
			{
				double gD;
				gD = Math.Round(g[i]);
				if (gD < 0) gD = 0;			// a zero digit with negative error can end up -1
				double error = g[i] - gD;
				g[i] = gD;
				if (Math.Abs(error) > Math.Abs(maxError))
					maxError = error;
			}
		}

		public static void ConvolutionCarry(double[] g, int ldn)
		{
			try
			{
				int n = 1 << ldn;
				long radixL = (long)radix;
				for (int i = 0; i < n; ++i)
				{

					if (g[i] >= radix)		// g[n-1] > radix is an overflow condition, should never happen
					{
						long rem;
						double carry = (double)Math.DivRem((long)g[i], radixL, out rem);
						g[i] = (double)rem;
						g[i+1] += carry;
					}
				}
			}
			catch (Exception ex)
			{
				Common.ReportError("ConvolutionCarry", ex.Message);
			}
		} // ConvolutionCarry

		public static double cosPiO8;			// cos(pi/8) COS_1_PI_8
		public static double piO4;				// pi over 4 M_PI_4
		public static double radix;				// radix = max value of a super digit + 1
		public static double sinPiO8;			// sin(pi/8) SIN_1_PI_8
		public static double sqrt1O2;			// 0.5**0.5 M_SQRT1_2
		public static double sqrt2;				// 2**0.5 M_SQRT2;
		public static int bitsPerInt;			// bits in an type int integer
		public static string dateTimeString;	// string for date and time
	}
}
