/*
 * test4100.redc - test REDC operations for test 41dd
 *
 * Copyright (C) 1999,2023  Ernest Bowen and Landon Curt Noll
 *
 * Primary author:  Ernest Bowen
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 *
 * Under source code control:	1996/03/13 03:53:22
 * File existed as early as:	1996
 *
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */

/*
 * Some severe tests and timing functions for REDC functions and pmod.
 *
 * testall(str,n,N,M,verbose)
 *	performs n tests using arguments x, y, ...
 *	randomly selected from [-N, N) or when nonnegative values are
 *	required, [0, N), and m an odd positive integer in [1,N],
 *	and where a "small" (say less than 10000) exponent k is to be
 *	used (when computing x^k % m directly) k is random in [0,M).
 *	Default values for N and M are 1e20 and 100.
 *
 * times(str,N,n,verbose)
 *	gives times for n evaluations of rcin, etc. with
 *	N-word arguments; default n is ceil(K1/power(N,1.585).
 *
 * powtimes(str, N1,N2,n, verbose)
 *	gives times for n evaluations of pmod(x,k,m)
 *	for the three algorithms "small", "normal", "bignum" that
 *	pmod may use, and equivalent functions rcpow(xr,k,m) for
 *	"small" or "bignum" cases, where xr = rcin(x,m).  The
 *	modulus m is a random positive odd N1-word number; x has
 *	random integer values in [0, m-1]; k has random N2-word values.
 *	N2 defaults to 1; n defaults to ceil(K2/power(N1,1.585)/N2).
 *
 * inittimes(str, N, n, verbose)
 *	displays the times and tests n evaluations of rcin(x,m)
 *	and rcout(x,m) where m is a random positive odd N-word integer,
 *	x is a random integer in [0, m-1]; n defaults to ceil(K1/N^2).
 *
 * rlen_4100(N)
 *	generates a random positive N-word integer.  The global
 *	BASEB should be set to the word-size for the computer being
 *	used.  The parameters K1, K2 which control the default n
 *	should be adjusted to give reasonable runtimes.
 *
 * olen(N)
 *	generates a random odd positive N-word number.
 *
 */


defaultverbose = 1;	/* default verbose value */

/*
 * test defaults
 */
global test4100_K1 = 2^17;
global test4100_K2 = 2^12;
global test4100_BASE = 2^config("baseb");

define rlen_4100(N) = rand(test4100_BASE^(N-1), test4100_BASE^N);

define olen(N)
{
	local v;

	v = rlen_4100(N);
	if (iseven(v))
		v++;
	return v;
}

define test4101(x,y,m,k,z1,z2,verbose)
{
	local xr, yr, v, w, oneone;

	if (isnull(verbose))
		verbose = defaultverbose;
	xr = rcin(x,m);
	yr = rcin(y,m);
	oneone = rcin(rcin(1,m),m);

	if (xr >= m || xr < 0) {
		if (verbose > 1)
			printf("Failure 1 for x = %d, m = %d\n", x, m);
		return 1;
	}
	if (rcin(x * y, m) != mod(xr * y, m, 0)) {
		if (verbose > 1) {
			printf("Failure 2 for x = %d, y = %d, m = %d\n",
				x, y, m);
		}
		return 2;
	}
	if (rcout(xr, m) != x % m) {
		if (verbose > 1)
			printf("Failure 3 for x = %d, m = %d\n", x, m);
		return 3;
	}
	if (rcout(rcmul(xr,yr,m),m) != mod(x * y, m, 0)) {
		if (verbose > 1)
			printf("Failure 4 for x = %d, y = %d, m = %d\n",
				 x, y, m);
		return 4;
	}
	if (rcmul(x,yr,m) != mod(x * y, m, 0)) {
		if (verbose > 1)
			printf("Failure 5 for x = %d, y = %d, m = %d\n",
				 x, y, m);
		return 5;
	}
	if (rcin(rcmul(x,y,m),m) != mod(x * y, m, 0)) {
		if (verbose > 1)
			printf("Failure 6 for x = %d, y = %d, m = %d\n",
				x, y, m);
		return 6;
	}
	if (rcout(rcsq(xr,m),m) != mod(x^2, m, 0)) {
		if (verbose > 1)
			printf("Failure 7 for x = %d, m = %d\n", x, m);
		return 7;
	}
	if (rcin(rcsq(x,m),m) != mod(x^2, m, 0)) {
		if (verbose > 1)
			printf("Failure 8 for x = %d, m = %d\n",
				x, y, m);
		return 8;
	}
	if (rcout(rcpow(xr,k,m),m) != mod(x^k, m, 0)) {
		if (verbose > 1)
			printf("Failure 9 for x = %d, m = %d, k = %d\n",
				x, m, k);
		return 9;
	}
	if (rcpow(x,k,m) != rcin(rcout(x,m)^k, m)) {
		if (verbose > 1)
			printf("Failure 10: x = %d, k = %d, m = %d\n",
				x, k, m);
		return 10;
	}
	if (rcpow(x, z1 * z2, m) != rcpow(rcpow(x,z1,m), z2, m)) {
		if (verbose > 1)
			printf("Failure 11: x = %d, z1 = %d, z2 = %d, m = %d\n",
			x, z1, z2, m);
		return 11;
	}
	if (xr != rcmul(oneone, x, m)) {
		if (verbose > 1)
			printf("Failure 12: x = %d, m = %d\n", x, m);
		return 12;
	}

	return 0;
}

define testall(str,n,N,M,verbose)
{
	local i, p, x, y, z1, z2, c, k, m;

	if (isnull(verbose))
		verbose = defaultverbose;
	if (verbose > 0) {
		print str:":",:;
	}
	m = 0;
	if (isnull(N))
		N = 1e20;
	if (isnull(M))
		M = 100;
	c = 0;
	for (i = 0; i < n; i++) {
		x = rand(-N, N);
		y = rand(-N, N);
		z1 = rand(N);
		z2 = rand(N);
		c = rand(N);
		if (iseven(c))
			c++;
		k = rand(M);
		if (verbose > 1)
			printf("x = %d, y = %d, c = %d, k = %d\n", x, y, c, k);
		p = test4101(x,y,c,k,z1,z2);
		if (p) {
			m++;
			if (verbose > 0) {
				printf("*** Failure %d in test %d\n", p, i);
			}
		}
	}
	if (verbose > 0) {
		if (m) {
			printf("*** %d error(s)\n", m);
		} else {
			printf("passed %d tests\n", n);
		}
	}
	return m;
}

define times(str,N,n,verbose)
{
	local m, m2, A, B, C, x, y, t, i, z;
	local trcin, trcout, trcmul, trcsq;
	local tmul, tsq, tmod, tquomod;

	if (isnull(verbose))
		verbose = defaultverbose;
	if (verbose > 0) {
		print str:":",:;
	}
	m = olen(N);
	m2 = m^2;
	if (isnull(n)) {
		n = ceil(test4100_K1/power(N,1.585));
		if (verbose > 1)
			printf("n = %d\n", n);
	}
	mat A[n];
	mat B[n];
	mat C[n];
	for (i = 0; i < n; i++) {
		A[i] = rand(m);
		B[i] = rand(m);
		C[i] = rand(m2);
	}
	z = rcin(0,m);	/* to initialize redc and maybe lastmod information */
	t = usertime();
	for (i = 0; i < n; i++)
		z = rcin(A[i],m);
	trcin = round(usertime() - t, 3);
	t = usertime();
	for (i = 0; i < n; i++)
		z = rcout(A[i],m);
	trcout = round(usertime() - t, 3);
	t = usertime();
	for (i = 0; i < n; i++)
		z = rcmul(A[i],B[i],m);
	trcmul = round(usertime() - t, 3);
	t = usertime();
	for (i = 0; i < n; i++)
		z = rcsq(A[i],m);
	trcsq = round(usertime() - t, 3);
	t = usertime();
	for (i = 0; i < n; i++)
		z = A[i] * B[i];
	tmul = round(usertime() - t, 3);
	t = usertime();
	for (i = 0; i < n; i++)
		z = A[i]^2;
	tsq = round(usertime() - t, 3);
	t = usertime();
	for (i = 0; i < n; i++)
		z = C[i] % A[i];
	tmod = round(usertime() - t, 3);
	t = usertime();
	for (i = 0; i < n; i++)
		quomod(C[i], A[i], x, y);
	tquomod = round(usertime() - t,3);

	if (verbose > 1) {
		printf("rcin: %d, rcout: %d, rcmul: %d, rcsq: %d\n",
			trcin, trcout, trcmul, trcsq);
		printf("%s: mul: %d, sq: %d, mod: %d, quomod: %d\n",
			str, tmul, tsq, tmod, tquomod);
	} else if (verbose > 0) {
		printf("no error(s)\n");
	}
	return 0;
}


define powtimes(str, N1, N2, n, verbose)
{
	local A, Ar, B, v, i, t, z1, z2, z3, z4, z5, cp, crc;
	local tsmall, tnormal, tbignum, trcsmall, trcbig, m;

	if (isnull(verbose))
		verbose = defaultverbose;
	if (verbose > 0) {
		print str:":",:;
	}
	m = 0;

	if (isnull(N2))
		N2 = 1;

	if (isnull(n)) {
		n = ceil(test4100_K2/power(N1, 1.585)/N2);
		printf ("n  = %d\n", n);
	}
	mat A[n];
	mat Ar[n];
	mat B[n];
	v = olen(N1);

	cp = config("pow2", 2);
	crc = config("redc2", 2);

	/* initialize redc and lastmod info */

	z1 = z2 = z3 = z4 = z5 = rcin(0,v);

	for (i = 0; i < n; i++) {
		A[i] = rand(v);
		Ar[i] = rcin(A[i], v);
		B[i] = rlen_4100(N2);
	}
	t = usertime();
	for (i = 0; i < n; i++)
		z1 += pmod(A[i], B[i], v);
	tbignum = round(usertime() - t, 4);
	config("pow2", 1e6);
	t = usertime();
	for (i = 0; i < n; i++)
		z2 += pmod(A[i], B[i], v);
	tnormal = round(usertime() - t, 4);
	config("redc2",1e6);
	t = usertime();
	for (i = 0; i < n; i++)
		z3 += pmod(A[i], B[i], v);
	tsmall = round(usertime() - t, 4);
	t = usertime();
	for (i = 0; i < n; i++)
		z4 += rcpow(Ar[i], B[i], v);
	trcsmall = round(usertime() - t, 4);
	config("redc2", 2);
	t = usertime();
	for (i = 0; i < n; i++)
		z5 += rcpow(Ar[i], B[i], v);
	trcbig = round(usertime() - t, 4);

	if (z1 != z2) {
		++m;
		if (verbose > 0) {
			print "*** z1 != z2";
		}
	} else if (z1 != z3) {
		++m;
		if (verbose > 0) {
			print "*** z1 != z3";
		}
	} else if (z2 != z3) {
		++m;
		if (verbose > 0) {
			print "*** z2 != z3";
		}
	} else if (rcout(z4, v) != z1 % v) {
		++m;
		if (verbose > 0) {
			print "*** z4 != z1";
		}
	} else if (z4 != z5) {
		++m;
		if (verbose > 0) {
			print "*** z4 != z5";
		}
	}
	config("pow2", cp);
	config("redc2", crc);
	if (verbose > 1) {
	}
	if (verbose > 1) {
		printf("Small: %d, normal: %d, bignum: %d\n",
			tsmall, tnormal, tbignum);
		printf("%s: rcsmall: %d, rcbig: %d\n",
			str, trcsmall, trcbig);
	} else if (verbose > 0) {
		if (m) {
			printf("*** %d error(s)\n", m);
		} else {
			printf("passed\n");
		}
	}
	return 0;
}

define inittimes(str,N,n,verbose)
{
	local A, M, B, R, i, t, trcin, trcout, m;

	if (isnull(verbose))
		verbose = defaultverbose;
	if (verbose > 0) {
		print str:":",:;
	}
	m = 0;
	if (isnull(n)) {
		n = ceil(test4100_K1/N^2);
		if (verbose > 1) {
			printf ("n = %d\n", n);
		}
	}
	mat A[n];
	mat M[n];
	mat B[n];
	mat R[n];
	for (i = 0; i < n; i++) {
		M[i] = olen(N);
		A[i] = rand(M[i]);
	}
	t = usertime();
	for (i = 0; i < n; i++)
		R[i] = rcin(A[i], M[i]);
	trcin = round(usertime() - t, 4);
	for (i = 0; i < n; i++)
		B[i] = rcout(R[i], M[i]);
	trcout = round(usertime() - t, 4);
	for (i = 0; i < n; i++) {
		if (B[i] != A[i]) {
			++m;
			if (verbose > 0) {
				print "*** Error!!!";
			}
			break;
		}
	}
	if (verbose > 0) {
		if (m) {
			printf("*** %d error(s)?\n", m);
		} else {
			if (verbose > 1) {
			    printf("%d successful tests: rcin: %d, rcout: %d\n",
			      n, trcin, trcout);
			} else {
			    printf("%d successful tests: passed\n", n);
			}
		}
	}
	return m;
}

/*
 * test4100 - perform all of the above tests a bunch of times
 */
define test4100(v, tnum)
{
	local n;	/* test parameter */

	/*
	 * set test parameters
	 */
	srand(4100e4100);

	/*
	 * test a lot of stuff
	 */
	err += testall(strcat(str(tnum++),": testall(10,,500)"), 10,, 500, v);
	err += testall(strcat(str(tnum++),": testall(200)"), 200,,, v);

	err += times(strcat(str(tnum++),": times(3,3000)"), 3, 3000, v);
	err += times(strcat(str(tnum++),": times(30,300)"), 30, 300, v);
	err += times(strcat(str(tnum++),": times(300,30)"), 300, 30, v);
	err += times(strcat(str(tnum++),": times(1000,3)"), 1000, 3, v);

	err += powtimes(strcat(str(tnum++),": powtimes(100)"),100,,v);
	err += powtimes(strcat(str(tnum++),": powtimes(250)"),250,,v);

	err += inittimes(strcat(str(tnum++),": inittimes(10)"),10,,v);
	err += inittimes(strcat(str(tnum++),": inittimes(100,70)"),100,70,v);
	err += inittimes(strcat(str(tnum++),": inittimes(1000,4)"),1000,4,v);

	/*
	 * report results
	 */
	if (v > 1) {
		if (err) {
			print "***", err, "error(s) found in testall";
		} else {
			print "no errors in testall";
		}
	}
	return tnum;
}