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