calc/cal/test2700.cal

/*
 * test2700 - 2700 series of the regress.cal test suite
 *
 * Copyright (C) 1999  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:	1995/11/01 22:52:25
 * File existed as early as:	1995
 *
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */

/*
 * The following resource file gives a severe test of sqrt(x,y,z) for
 * all 128 values of z, randomly produced real and complex x, and randomly
 * produced nonzero values for y.  After loading it, testcsqrt(n) will
 * test n combinations of x and y;  testcsqrt(str,n,2) will print 1 2 3 ...
 * indicating work in process; testcsqrt(str,n,3) will give information about
 * errors detected and will print values of x and y used.
 * I've also defined a function iscomsq(x) which does for complex as well
 * as real x what issq(x) currently does for real x.
 */


defaultverbose = 1;

define mknonnegreal() {
	switch(rand(8)) {
		case 0: return rand(20);
		case 1: return rand(20,1000);
		case 2: return rand(1,10000)/rand(1,100);
		case 3: return scale(mkposreal(), rand(1,100));
		case 4: return scale(mkposreal(), -rand(1,100));
		case 5: return rand(1, 1000) + scale(mkfrac(),-rand(1,100));
		case 6: return mkposreal()^2;
		case 7: return mkposreal() * (1+scale(mkfrac(),-rand(1,100)));
	}
}

define mkposreal() {
	local x;

	x = mknonnegreal();
	while (x == 0)
		x = mknonnegreal();
	return x;
}

define mkreal_2700() = rand(2) ? mknonnegreal() : -mknonnegreal();

define mknonzeroreal() = rand(2) ? mkposreal() : -mkposreal();

/* Number > 0 and < 1, almost uniformly distributed */
define mkposfrac() {
	local x,y;

	x = rand(1,1000);
	do
		y = rand(1,1000);
	while (y == x);
	if (x > y)
		swap(x,y);
	return x/y;
}

/* Nonzero > -1 and < 1 */
define mkfrac() = rand(2) ? mkposfrac() : -mkposfrac();

define mksquarereal() = mknonnegreal()^2;

/*
 * We might be able to do better than the following.  For nonsquare
 * positive integer less than 1e6, could use:
 *		 x = rand(1, 1000);
 *		 return rand(x^2 + 1, (x + 1)^2);
 * Maybe could do:
 *		do
 *			x = mkreal_2700();
 *		while
 *			(issq(x));
 * This would of course not be satisfactory for testing issq().
 */

define mknonsquarereal() = 22 * mkposreal()^2/7;

define mkcomplex_2700() = mkreal_2700() + 1i * mkreal_2700();

define testcsqrt(str, n, verbose)
{
	local x, y, z, m, i, p, v;

	if (isnull(verbose))
		verbose = defaultverbose;
	if (verbose > 0) {
		print str:":",:;
	}
	m = 0;
	for (i = 1; i <= n; i++) {
		if (verbose > 1) print i,:;
		x = rand(3) ? mkreal_2700(): mkcomplex_2700();
		y = scale(mknonzeroreal(), -100);
		if (verbose > 2)
			printf("%d: x = %d, y = %d\n", i, x, y);

		for (z = 0; z < 128; z++) {
			v = sqrt(x,y,z);
			p = checksqrt(x,y,z,v);
			if (p) {
			if (verbose > 0)
				printf(
				 "*** Type %d failure for x = %r, "
				 "y = %r, z = %d\n",
				    p, x, y, z);
				m++;
			}
		}
	}
	if (verbose > 0) {
		if (m) {
			printf("*** %d error(s)\n", m);
		} else {
			printf("no errors\n");
		}
	}
	return m;
}


define checksqrt(x,y,z,v)	/* Returns >0 if an error is detected */
{
	local A, B, X, Y, t1, t2, eps, u, n, f, s;

	A = re(x);
	B = im(x);
	X = re(v);
	Y = im(v);

	/* checking signs of X and Y */

	if (B == 0 && A <= 0)		/* t1 = sgn(re(tvsqrt)) */
		t1 = 0;
	else
		t1 = (z & 64) ? -1 : 1;

	t2 = B ? sgn(B) : (A < 0);	/* t2 = sgn(im(tvsqrt)) */
	if (z & 64)
		t2 = -t2;

	if (t1 == 0 && X != 0)
		return 1;

	if (t2 == 0 && Y != 0) {
		printf("x = %d, Y = %d, t2 = %d\n", x, Y, t2);
		return 2;
	}

	if (X && sgn(X) != t1)
		return 3;

	if (Y && sgn(Y) != t2)
		return 4;

	if (z & 32 && iscomsq(x))
		return 5 * (x != v^2);

	eps = (z & 16) ? abs(y)/2 : abs(y);
	u = sgn(y);

	/* Checking X */

	n = X/y;
	if (!isint(n))
		return 6;

	if (t1) {
		f = checkavrem(A, B, abs(X), eps);

		if (z & 16 && f < 0)
			return 7;
		if (!(z & 16) && f <= 0)
			return 8;

		if (!(z & 16) || f == 0) {
			s = X ? t1 * sgn(A - X^2 + B^2/4/X^2) : t1;
			if (s && !checkrounding(s,n,t1,u,z))
			return 9;
		}
	}

	/* Checking Y */

	n = Y/y;
	if (!isint(n))
		return 10;

	if (t2) {
		f = checkavrem(-A, B, abs(Y), eps);

		if (z & 16 && f < 0)
			return 11;
		if (!(z & 16) && f <= 0)
			return 12;

		if (!(z & 16) || f == 0) {
			s = Y ? t2 * sgn(-A - Y^2 + B^2/4/Y^2) : t2;
			if (s && !checkrounding(s,n,t2,u,z))
				return 13;
		}
	}
	return 0;
}

/*
 * Check that the calculated absolute value X of the real part of
 * sqrt(A + Bi) is between (true value - eps) and (true value + eps).
 * Returns -1 if it is outside, 0 if on boundary, 1 if between.
 */

define checkavrem(A, B, X, eps)
{
	local f;

	f = sgn(A - (X + eps)^2 + B^2/4/(X + eps)^2);
	if (f > 0)
		return -1;		/* X < tv - eps */
	if (f == 0)
		return 0;		/* X = tv - eps */

	if (X > eps) {
		f = sgn(A - (X - eps)^2 + B^2/4/(X - eps)^2);

		if (f < 0)
			return -1;	/* X > tv + eps */
		if (f == 0)
			return 0;	/* X = tv + eps */
	}
	return 1;		/* tv - eps < X < tv + eps */
}


define checkrounding(s,n,t,u,z)
{
	local w;

	switch (z & 15) {
		case 0: w = (s == u); break;
		case 1: w = (s == -u); break;
		case 2: w = (s == t); break;
		case 3: w = (s == -t); break;
		case 4: w = (s > 0); break;
		case 5: w = (s < 0); break;
		case 6: w = (s == u/t); break;
		case 7: w = (s == -u/t); break;
		case 8: w = iseven(n); break;
		case 9: w = isodd(n); break;
		case 10: w = (u/t > 0) ? iseven(n) : isodd(n); break;
		case 11: w = (u/t > 0) ? isodd(n) : iseven(n); break;
		case 12: w = (u > 0) ? iseven(n) : isodd(n); break;
		case 13: w = (u > 0) ? isodd(n) : iseven(n); break;
		case 14: w = (t > 0) ? iseven(n) : isodd(n); break;
		case 15: w = (t > 0) ? isodd(n) : iseven(n); break;
	}
	return w;
}

define iscomsq(x)
{
	local c;

	if (isreal(x))
		return issq(abs(x));
	c = norm(x);
	if (!issq(c))
		return 0;
	return issq((re(x) + sqrt(c,1,32))/2);
}

/*
 * test2700 - perform all of the above tests a bunch of times
 */
define test2700(verbose, tnum)
{
	local n;	/* test parameter */
	local ep;	/* test parameter */
	local i;

	/* set test parameters */
	n = 32;		/* internal test loop count */
	if (isnull(verbose)) {
		verbose = defaultverbose;
	}
	if (isnull(tnum)) {
		tnum = 1;	/* initial test number */
	}

	/*
	 * test a lot of stuff
	 */
	srand(2700e2700);
	for (i=0; i < n; ++i) {
		err += testcsqrt(strcat(str(tnum++),": complex sqrt",str(i)),
				 1, verbose);
	}
	if (verbose > 1) {
		if (err) {
			print "***", err, "error(s) found in testall";
		} else {
			print "no errors in testall";
		}
	}
	return tnum;
}