1/*
2 * test2700 - 2700 series of the regress.cal test suite
3 *
4 * Copyright (C) 1999,2021  Ernest Bowen and Landon Curt Noll
5 *
6 * Primary author:  Ernest Bowen
7 *
8 * Calc is open software; you can redistribute it and/or modify it under
9 * the terms of the version 2.1 of the GNU Lesser General Public License
10 * as published by the Free Software Foundation.
11 *
12 * Calc is distributed in the hope that it will be useful, but WITHOUT
13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
15 * Public License for more details.
16 *
17 * A copy of version 2.1 of the GNU Lesser General Public License is
18 * distributed with calc under the filename COPYING-LGPL.  You should have
19 * received a copy with calc; if not, write to Free Software Foundation, Inc.
20 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
21 *
22 * Under source code control:	1995/11/01 22:52:25
23 * File existed as early as:	1995
24 *
25 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
26 */
27
28/*
29 * The following resource file gives a severe test of sqrt(x,y,z) for
30 * all 128 values of z, randomly produced real and complex x, and randomly
31 * produced nonzero values for y.  After loading it, testcsqrt(n) will
32 * test n combinations of x and y;  testcsqrt(str,n,2) will print 1 2 3 ...
33 * indicating work in process; testcsqrt(str,n,3) will give information about
34 * errors detected and will print values of x and y used.
35 * I've also defined a function iscomsq(x) which does for complex as well
36 * as real x what issq(x) currently does for real x.
37 */
38
39
40defaultverbose = 1;
41
42define mknonnegreal() {
43	switch(rand(8)) {
44		case 0: return rand(20);
45		case 1: return rand(20,1000);
46		case 2: return rand(1,10000)/rand(1,100);
47		case 3: return scale(mkposreal(), rand(1,100));
48		case 4: return scale(mkposreal(), -rand(1,100));
49		case 5: return rand(1, 1000) + scale(mkfrac(),-rand(1,100));
50		case 6: return mkposreal()^2;
51		case 7: return mkposreal() * (1+scale(mkfrac(),-rand(1,100)));
52	}
53}
54
55define mkposreal() {
56	local x;
57
58	x = mknonnegreal();
59	while (x == 0)
60		x = mknonnegreal();
61	return x;
62}
63
64define mkreal_2700() = rand(2) ? mknonnegreal() : -mknonnegreal();
65
66define mknonzeroreal() = rand(2) ? mkposreal() : -mkposreal();
67
68/* Number > 0 and < 1, almost uniformly distributed */
69define mkposfrac() {
70	local x,y;
71
72	x = rand(1,1000);
73	do
74		y = rand(1,1000);
75	while (y == x);
76	if (x > y)
77		swap(x,y);
78	return x/y;
79}
80
81/* Nonzero > -1 and < 1 */
82define mkfrac() = rand(2) ? mkposfrac() : -mkposfrac();
83
84define mksquarereal() = mknonnegreal()^2;
85
86/*
87 * We might be able to do better than the following.  For non-square
88 * positive integer less than 1e6, could use:
89 *		 x = rand(1, 1000);
90 *		 return rand(x^2 + 1, (x + 1)^2);
91 * Maybe could do:
92 *		do
93 *			x = mkreal_2700();
94 *		while
95 *			(issq(x));
96 * This would of course not be satisfactory for testing issq().
97 */
98
99define mknonsquarereal() = 22 * mkposreal()^2/7;
100
101define mkcomplex_2700() = mkreal_2700() + 1i * mkreal_2700();
102
103define testcsqrt(str, n, verbose)
104{
105	local x, y, z, m, i, p, v;
106
107	if (isnull(verbose))
108		verbose = defaultverbose;
109	if (verbose > 0) {
110		print str:":",:;
111	}
112	m = 0;
113	for (i = 1; i <= n; i++) {
114		if (verbose > 1) print i,:;
115		x = rand(3) ? mkreal_2700(): mkcomplex_2700();
116		y = scale(mknonzeroreal(), -100);
117		if (verbose > 2)
118			printf("%d: x = %d, y = %d\n", i, x, y);
119
120		for (z = 0; z < 128; z++) {
121			v = sqrt(x,y,z);
122			p = checksqrt(x,y,z,v);
123			if (p) {
124			if (verbose > 0)
125				printf(
126				 "*** Type %d failure for x = %r, "
127				 "y = %r, z = %d\n",
128				    p, x, y, z);
129				m++;
130			}
131		}
132	}
133	if (verbose > 0) {
134		if (m) {
135			printf("*** %d error(s)\n", m);
136		} else {
137			printf("no errors\n");
138		}
139	}
140	return m;
141}
142
143
144define checksqrt(x,y,z,v)	/* Returns >0 if an error is detected */
145{
146	local A, B, X, Y, t1, t2, eps, u, n, f, s;
147
148	A = re(x);
149	B = im(x);
150	X = re(v);
151	Y = im(v);
152
153	/* checking signs of X and Y */
154
155	if (B == 0 && A <= 0)		/* t1 = sgn(re(tvsqrt)) */
156		t1 = 0;
157	else
158		t1 = (z & 64) ? -1 : 1;
159
160	t2 = B ? sgn(B) : (A < 0);	/* t2 = sgn(im(tvsqrt)) */
161	if (z & 64)
162		t2 = -t2;
163
164	if (t1 == 0 && X != 0)
165		return 1;
166
167	if (t2 == 0 && Y != 0) {
168		printf("x = %d, Y = %d, t2 = %d\n", x, Y, t2);
169		return 2;
170	}
171
172	if (X && sgn(X) != t1)
173		return 3;
174
175	if (Y && sgn(Y) != t2)
176		return 4;
177
178	if (z & 32 && iscomsq(x))
179		return 5 * (x != v^2);
180
181	eps = (z & 16) ? abs(y)/2 : abs(y);
182	u = sgn(y);
183
184	/* Checking X */
185
186	n = X/y;
187	if (!isint(n))
188		return 6;
189
190	if (t1) {
191		f = checkavrem(A, B, abs(X), eps);
192
193		if (z & 16 && f < 0)
194			return 7;
195		if (!(z & 16) && f <= 0)
196			return 8;
197
198		if (!(z & 16) || f == 0) {
199			s = X ? t1 * sgn(A - X^2 + B^2/4/X^2) : t1;
200			if (s && !checkrounding(s,n,t1,u,z))
201			return 9;
202		}
203	}
204
205	/* Checking Y */
206
207	n = Y/y;
208	if (!isint(n))
209		return 10;
210
211	if (t2) {
212		f = checkavrem(-A, B, abs(Y), eps);
213
214		if (z & 16 && f < 0)
215			return 11;
216		if (!(z & 16) && f <= 0)
217			return 12;
218
219		if (!(z & 16) || f == 0) {
220			s = Y ? t2 * sgn(-A - Y^2 + B^2/4/Y^2) : t2;
221			if (s && !checkrounding(s,n,t2,u,z))
222				return 13;
223		}
224	}
225	return 0;
226}
227
228/*
229 * Check that the calculated absolute value X of the real part of
230 * sqrt(A + Bi) is between (true value - eps) and (true value + eps).
231 * Returns -1 if it is outside, 0 if on boundary, 1 if between.
232 */
233
234define checkavrem(A, B, X, eps)
235{
236	local f;
237
238	f = sgn(A - (X + eps)^2 + B^2/4/(X + eps)^2);
239	if (f > 0)
240		return -1;		/* X < tv - eps */
241	if (f == 0)
242		return 0;		/* X = tv - eps */
243
244	if (X > eps) {
245		f = sgn(A - (X - eps)^2 + B^2/4/(X - eps)^2);
246
247		if (f < 0)
248			return -1;	/* X > tv + eps */
249		if (f == 0)
250			return 0;	/* X = tv + eps */
251	}
252	return 1;		/* tv - eps < X < tv + eps */
253}
254
255
256define checkrounding(s,n,t,u,z)
257{
258	local w;
259
260	switch (z & 15) {
261		case 0: w = (s == u); break;
262		case 1: w = (s == -u); break;
263		case 2: w = (s == t); break;
264		case 3: w = (s == -t); break;
265		case 4: w = (s > 0); break;
266		case 5: w = (s < 0); break;
267		case 6: w = (s == u/t); break;
268		case 7: w = (s == -u/t); break;
269		case 8: w = iseven(n); break;
270		case 9: w = isodd(n); break;
271		case 10: w = (u/t > 0) ? iseven(n) : isodd(n); break;
272		case 11: w = (u/t > 0) ? isodd(n) : iseven(n); break;
273		case 12: w = (u > 0) ? iseven(n) : isodd(n); break;
274		case 13: w = (u > 0) ? isodd(n) : iseven(n); break;
275		case 14: w = (t > 0) ? iseven(n) : isodd(n); break;
276		case 15: w = (t > 0) ? isodd(n) : iseven(n); break;
277	}
278	return w;
279}
280
281define iscomsq(x)
282{
283	local c;
284
285	if (isreal(x))
286		return issq(abs(x));
287	c = norm(x);
288	if (!issq(c))
289		return 0;
290	return issq((re(x) + sqrt(c,1,32))/2);
291}
292
293/*
294 * test2700 - perform all of the above tests a bunch of times
295 */
296define test2700(verbose, tnum)
297{
298	local n;	/* test parameter */
299	local ep;	/* test parameter */
300	local i;
301
302	/* set test parameters */
303	n = 32;		/* internal test loop count */
304	if (isnull(verbose)) {
305		verbose = defaultverbose;
306	}
307	if (isnull(tnum)) {
308		tnum = 1;	/* initial test number */
309	}
310
311	/*
312	 * test a lot of stuff
313	 */
314	srand(2700e2700);
315	for (i=0; i < n; ++i) {
316		err += testcsqrt(strcat(str(tnum++),": complex sqrt",str(i)),
317				 1, verbose);
318	}
319	if (verbose > 1) {
320		if (err) {
321			print "***", err, "error(s) found in testall";
322		} else {
323			print "no errors in testall";
324		}
325	}
326	return tnum;
327}
328