1/* 2 * intnum - implementation of tanhsinh- and Gauss-Legendre quadrature 3 * 4 * Copyright (C) 2013,2021 Christoph Zurnieden 5 * 6 * Calc is open software; you can redistribute it and/or modify it under 7 * the terms of the version 2.1 of the GNU Lesser General Public License 8 * as published by the Free Software Foundation. 9 * 10 * Calc is distributed in the hope that it will be useful, but WITHOUT 11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 12 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General 13 * Public License for more details. 14 * 15 * A copy of version 2.1 of the GNU Lesser General Public License is 16 * distributed with calc under the filename COPYING-LGPL. You should have 17 * received a copy with calc; if not, write to Free Software Foundation, Inc. 18 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 19 */ 20 21 22static resource_debug_level; 23resource_debug_level = config("resource_debug", 0); 24 25 26read -once infinities; 27 28static __CZ__tanhsinh_x; 29static __CZ__tanhsinh_w; 30static __CZ__tanhsinh_order; 31static __CZ__tanhsinh_prec; 32 33define quadtsdeletenodes() 34{ 35 free(__CZ__tanhsinh_x); 36 free(__CZ__tanhsinh_w); 37 free(__CZ__tanhsinh_order); 38 free(__CZ__tanhsinh_prec); 39} 40 41define quadtscomputenodes(order, expo, eps) 42{ 43 local t cht sht chp sum k PI places; 44 local h t0 x w; 45 if (__CZ__tanhsinh_order == order && __CZ__tanhsinh_prec == eps) 46 return 1; 47 __CZ__tanhsinh_order = order; 48 __CZ__tanhsinh_prec = eps; 49 __CZ__tanhsinh_x = list(); 50 __CZ__tanhsinh_w = list(); 51 /* The tanhsinh algorithm needs a slightly higher precision than G-L */ 52 eps = epsilon(eps * 1e-2); 53 places = highbit(1 + int (1 / epsilon())) +1; 54 PI = pi(); 55 sum = 0; 56 t0 = 2 ^ (-expo); 57 h = 2 * t0; 58 /* 59 * The author wanted to use the mpmath trick here which was 60 * advertised---and reasonably so!---to be faster. Didn't work out 61 * so well with calc. 62 * PI4 = PI/4; 63 * expt0 = bround(exp(t0),places); 64 * a = bround( PI4 * expt0,places); 65 * b = bround(PI4 / expt0,places); 66 * udelta = bround(exp(h),places); 67 * urdelta = bround(1/udelta,places); 68 */ 69 /* make use of x(-t) = -x(t), w(-t) = w(t) */ 70 for (k = 0; k < 20 * order + 1; k++) { 71 /* 72 * x = tanh(pi/2 * sinh(t)) 73 * w = pi/2 * cosh(t) / cosh(pi/2 * sinh(t))^2 74 */ 75 t = bround(t0 + k * h, places); 76 77 cht = bround(cosh(t), places); 78 sht = bround(sinh(t), places); 79 chp = bround(cosh(0.5 * PI * sht), places); 80 x = bround(tanh(0.5 * PI * sht), places); 81 w = bround((PI * h * cht) / (2 * chp ^ 2), places); 82 /* 83 * c = bround(exp(a-b),places); 84 * d = bround(1/c,places); 85 * co =bround( (c+d)/2,places); 86 * si =bround( (c-d)/2,places); 87 * x = bround(si / co,places); 88 * w = bround((a+b) / co^2,places); 89 */ 90 if (abs(x - 1) <= eps) 91 break; 92 93 append(__CZ__tanhsinh_x, x); 94 append(__CZ__tanhsinh_w, w); 95 /* 96 * a *= udelta; 97 * b *= urdelta; 98 */ 99 } 100 101 /* Normalize the weights to make them add up to 2 (two) */ 102 /* 103 * for(k=0;k < size(__CZ__tanhsinh_w);k++) 104 * sum = bround(sum + __CZ__tanhsinh_w[k],places); 105 * sum *= 2; 106 * for(k=0;k < size(__CZ__tanhsinh_w);k++) 107 * __CZ__tanhsinh_w[k] = bround(2.0 * __CZ__tanhsinh_w[k] / sum,places); 108 */ 109 110 epsilon(eps); 111 return 1; 112} 113 114define quadtscore(a, b, n) 115{ 116 local k c d order eps places sum ret x x1 x2 xm w w1 w2 m sizel; 117 118 eps = epsilon(epsilon() * 1e-2); 119 places = highbit(1 + int (1 / epsilon())) +1; 120 m = int (4 + max(0, ln(places / 30.0) / ln(2))) + 2; 121 if (!isnull(n)) { 122 order = n; 123 m = ilog(order / 3, 2) + 1; 124 } else 125 order = 3 * 2 ^ (m - 1); 126 127 quadtscomputenodes(order, m, epsilon()); 128 sizel = size(__CZ__tanhsinh_w); 129 130 if (isinfinite(a) || isinfinite(b)) { 131 /* 132 * x 133 * t = ------------ 134 * 2 135 * sqrt(1 - y ) 136 */ 137 if (isninf(a) && ispinf(b)) { 138 for (k = 0; k < sizel; k++) { 139 x1 = __CZ__tanhsinh_x[k]; 140 x2 = -__CZ__tanhsinh_x[k]; 141 w1 = __CZ__tanhsinh_w[k]; 142 143 x = bround(x1 * (1 - x1 ^ 2) ^ (-1 / 2), places); 144 xm = bround(x2 * (1 - x2 ^ 2) ^ (-1 / 2), places); 145 w = bround(w1 * (((1 - x1 ^ 2) ^ (-1 / 2)) / (1 - x1 ^ 2)), 146 places); 147 w2 = bround(w1 * (((1 - x2 ^ 2) ^ (-1 / 2)) / (1 - x2 ^ 2)), 148 places); 149 sum += bround(w * f(x), places); 150 sum += bround(w2 * f(xm), places); 151 } 152 } 153 /* 154 * 1 155 * t = - - + b + 1 156 * x 157 */ 158 else if (isninf(a) && !iscinf(b)) { 159 for (k = 0; k < sizel; k++) { 160 x1 = __CZ__tanhsinh_x[k]; 161 x2 = -__CZ__tanhsinh_x[k]; 162 w1 = __CZ__tanhsinh_w[k]; 163 164 x = bround((b + 1) - (2 / (x1 + 1)), places); 165 xm = bround((b + 1) - (2 / (x2 + 1)), places); 166 w = bround(w1 * (1 / 2 * (2 / (x1 + 1)) ^ 2), places); 167 w2 = bround(w1 * (1 / 2 * (2 / (x2 + 1)) ^ 2), places); 168 sum += bround(w * f(x), places); 169 sum += bround(w2 * f(xm), places); 170 } 171 } 172 /* 173 * 1 174 * t = - + a - 1 175 * x 176 */ 177 else if (!iscinf(a) && ispinf(b)) { 178 for (k = 0; k < sizel; k++) { 179 x1 = __CZ__tanhsinh_x[k]; 180 x2 = -__CZ__tanhsinh_x[k]; 181 w1 = __CZ__tanhsinh_w[k]; 182 x = bround((a - 1) + (2 / (x1 + 1)), places); 183 xm = bround((a - 1) + (2 / (x2 + 1)), places); 184 w = bround(w1 * (((1 / 2) * (2 / (x1 + 1)) ^ 2)), places); 185 w2 = bround(w1 * (((1 / 2) * (2 / (x2 + 1)) ^ 2)), places); 186 sum += bround(w * f(x), places); 187 sum += bround(w2 * f(xm), places); 188 } 189 } else if (isninf(a) || isninf(b)) { 190 /*TODO: swap(a,b) and negate(w)? Lookup! */ 191 return newerror("quadtscore: reverse limits?"); 192 } else { 193 return 194 newerror("quadtscore: complex infinity not yet implemented"); 195 } 196 ret = sum; 197 } else { 198 /* Avoid rounding errors */ 199 if (a == -1 && b == 1) { 200 c = 1; 201 d = 0; 202 } else { 203 c = (b - a) / 2; 204 d = (b + a) / 2; 205 } 206 sum = 0; 207 for (k = 0; k < sizel; k++) { 208 sum += 209 bround(__CZ__tanhsinh_w[k] * f(c * __CZ__tanhsinh_x[k] + d), 210 places); 211 sum += 212 bround(__CZ__tanhsinh_w[k] * f(c * -__CZ__tanhsinh_x[k] + d), 213 places); 214 } 215 ret = c * sum; 216 } 217 epsilon(eps); 218 return ret; 219} 220 221static __CZ__quadts_error; 222 223define quadts(a, b, points) 224{ 225 local k sp results epsbits nsect interval length segment slope C ; 226 local x1 x2 y1 y2 sum D1 D2 D3 D4; 227 if (param(0) < 2) 228 return newerror("quadts: not enough arguments"); 229 epsbits = highbit(1 + int (1 / epsilon())) +1; 230 if (param(0) < 3 || isnull(points)) { 231 /* return as given */ 232 return quadtscore(a, b); 233 } else { 234 if ((isinfinite(a) || isinfinite(b)) 235 && (!ismat(points) && !islist(points))) 236 return 237 newerror(strcat 238 ("quadts: segments of infinite length ", 239 "are not yet supported")); 240 if (ismat(points) || islist(points)) { 241 sp = size(points); 242 if (sp == 0) 243 return 244 newerror(strcat 245 ("quadts: variable 'points` must be a list or ", 246 "1d-matrix of a length > 0")); 247 /* check if all points are numbers */ 248 for (k = 0; k < sp; k++) { 249 if (!isnum(points[k])) 250 return 251 newerror(strcat 252 ("quadts: elements of 'points` must be", 253 " numbers only")); 254 } 255 /* We have n-1 intervals and a and b, hence n-1 + 2 results */ 256 results = mat[sp + 1]; 257 if (a != points[0]) { 258 results[0] = quadtscore(a, points[0]); 259 } else { 260 results[0] = 0; 261 } 262 if (sp == 1) { 263 if (b != points[0]) { 264 results[1] = quadtscore(points[0], b); 265 } else { 266 results[1] = 0; 267 } 268 } else { 269 for (k = 1; k < sp; k++) { 270 results[k] = quadtscore(points[k - 1], points[k]); 271 } 272 if (b != points[k - 1]) { 273 results[k] = quadtscore(points[k - 1], b); 274 } else { 275 results[k] = 0; 276 } 277 } 278 } else { 279 if (!isint(points) || points <= 0) 280 return newerror(strcat("quadts: variable 'points` must be a ", 281 "list or a positive integer")); 282 /* Taking "points" as the number of equally spaced intervals */ 283 results = mat[points + 1]; 284 /* It is easy if a,b lie on the real line */ 285 if (isreal(a) && isreal(b)) { 286 length = abs(a - b); 287 segment = length / points; 288 289 for (k = 1; k <= points; k++) { 290 results[k - 1] = 291 quadtscore(a + (k - 1) * segment, a + k * segment); 292 } 293 } else { 294 /* We have at least one complex limit but treat "points" still 295 * as the number of equally spaced intervals on a straight line 296 * connecting a and b. Computing the segments here is a bit 297 * more complicated but not much, it should have been taught in 298 * high school. 299 * Other contours by way of a list of points */ 300 slope = (im(b) - im(a)) / (re(b) - re(a)); 301 C = (im(a) + slope) * re(a); 302 length = abs(re(a) - re(b)); 303 segment = length / points; 304 305 /* y = mx+C where m is the slope, x is the real part and y the 306 * imaginary part */ 307 if(re(a)>re(b))swap(a,b); 308 for (k = re(a); k <= (re(b)); k+=segment) { 309 x1 = slope*(k) + C; 310 results[k] = quadtscore(k + x1 * 1i); 311 } 312 } /* else of isreal */ 313 } /* else of ismat|islist */ 314 } /* else of isnull(points) */ 315 /* With a bit of undeserved luck we have a result by now. */ 316 sp = size(results); 317 for (k = 0; k < sp; k++) { 318 sum += results[k]; 319 } 320 return sum; 321} 322 323static __CZ__gl_x; 324static __CZ__gl_w; 325static __CZ__gl_order; 326static __CZ__gl_prec; 327 328define quadglcomputenodes(N) 329{ 330 local places k l x w t1 t2 t3 t4 t5 r tmp; 331 332 if (__CZ__gl_order == N && __CZ__gl_prec == epsilon()) 333 return; 334 335 __CZ__gl_x = mat[N]; 336 __CZ__gl_w = mat[N]; 337 __CZ__gl_order = N; 338 __CZ__gl_prec = epsilon(); 339 340 places = highbit(1 + int (1 / epsilon())) +1; 341 342 /* 343 * Compute roots and weights (doing it inline seems to be fastest) 344 * Trick shamelessly stolen from D. Bailey et .al (program "arprec") 345 */ 346 for (k = 1; k <= N//2; k++) { 347 r = bround(cos(pi() * (k - .25) / (N + .5)), places); 348 while (1) { 349 t1 = 1, t2 = 0; 350 for (l = 1; l <= N; l++) { 351 t3 = t2; 352 t2 = t1; 353 t1 = bround(((2 * l - 1) * r * t2 - (l - 1) * t3) / l, places); 354 } 355 t4 = bround(N * (r * t1 - t2) / ((r ^ 2) - 1), places); 356 t5 = r; 357 tmp = t1 / t4; 358 r = r - tmp; 359 if (abs(tmp) <= epsilon()) 360 break; 361 } 362 x = r; 363 w = bround(2 / ((1 - r ^ 2) * t4 ^ 2), places); 364 365 __CZ__gl_x[k - 1] = x; 366 __CZ__gl_w[k - 1] = w; 367 __CZ__gl_x[N - k] = -__CZ__gl_x[k - 1]; 368 __CZ__gl_w[N - k] = __CZ__gl_w[k - 1]; 369 } 370 return; 371} 372 373define quadgldeletenodes() 374{ 375 free(__CZ__gl_x); 376 free(__CZ__gl_w); 377 free(__CZ__gl_order); 378 free(__CZ__gl_prec); 379} 380 381define quadglcore(a, b, n) 382{ 383 local k c d digs order eps places sum ret err x x1 w w1 m; 384 local phalf x2 px1 spx1 u b1 a1 half; 385 386 eps = epsilon(epsilon() * 1e-2); 387 places = highbit(1 + int (1 / epsilon())) +1; 388 if (!isnull(n)) 389 order = n; 390 else { 391 m = int (4 + max(0, ln(places / 30.0) / ln(2))) + 2; 392 order = 3 * 2 ^ (m - 1); 393 } 394 395 396 quadglcomputenodes(order, 1); 397 398 if (isinfinite(a) || isinfinite(b)) { 399 if (isninf(a) && ispinf(b)) { 400 for (k = 0; k < order; k++) { 401 x1 = __CZ__gl_x[k]; 402 w1 = __CZ__gl_w[k]; 403 404 x = bround(x1 * (1 - x1 ^ 2) ^ (-1 / 2), places); 405 w = bround(w1 * (((1 - x1 ^ 2) ^ (-1 / 2)) / (1 - x1 ^ 2)), 406 places); 407 sum += bround(w * f(x), places); 408 } 409 } else if (isninf(a) && !iscinf(b)) { 410 for (k = 0; k < order; k++) { 411 x1 = __CZ__gl_x[k]; 412 w1 = __CZ__gl_w[k]; 413 414 x = bround((b + 1) - (2 / (x1 + 1)), places); 415 w = bround(w1 * (1 / 2 * (2 / (x1 + 1)) ^ 2), places); 416 sum += bround(w * f(x), places); 417 } 418 } else if (!iscinf(a) && ispinf(b)) { 419 for (k = 0; k < order; k++) { 420 x1 = __CZ__gl_x[k]; 421 w1 = __CZ__gl_w[k]; 422 x = bround((a - 1) + (2 / (x1 + 1)), places); 423 w = bround(w1 * (((1 / 2) * (2 / (x1 + 1)) ^ 2)), places); 424 sum += bround(w * f(x), places); 425 } 426 } else if (isninf(a) || isninf(b)) { 427 /*TODO: swap(a,b) and negate(w)? Lookup! */ 428 return newerror("quadglcore: reverse limits?"); 429 } else 430 return 431 newerror("quadglcore: complex infinity not yet implemented"); 432 ret = sum; 433 } else { 434 /* Avoid rounding errors */ 435 if (a == -1 && b == 1) { 436 c = 1; 437 d = 0; 438 } else { 439 c = (b - a) / 2; 440 d = (b + a) / 2; 441 } 442 sum = 0; 443 for (k = 0; k < order; k++) { 444 sum += bround(__CZ__gl_w[k] * f(c * __CZ__gl_x[k] + d), places); 445 } 446 ret = c * sum; 447 } 448 epsilon(eps); 449 return ret; 450} 451 452define quadgl(a, b, points) 453{ 454 local k sp results epsbits nsect interval length segment slope C x1 y1 x2 455 y2; 456 local sum D1 D2 D3 D4; 457 if (param(0) < 2) 458 return newerror("quadgl: not enough arguments"); 459 epsbits = highbit(1 + int (1 / epsilon())) +1; 460 if (isnull(points)) { 461 /* return as given */ 462 return quadglcore(a, b); 463 } else { 464 /* But if we could half the time needed to execute a single operation 465 * we could do all of it in just twice that time. */ 466 if (isinfinite(a) || isinfinite(b) 467 && (!ismat(points) && !islist(points))) 468 return 469 newerror(strcat 470 ("quadgl: multiple segments of infinite length ", 471 "are not yet supported")); 472 if (ismat(points) || islist(points)) { 473 sp = size(points); 474 if (sp == 0) 475 return 476 newerror(strcat 477 ("quadgl: variable 'points` must be a list or ", 478 "1d-matrix of a length > 0")); 479 /* check if all points are numbers */ 480 for (k = 0; k < sp; k++) { 481 if (!isnum(points[k])) 482 return 483 newerror(strcat 484 ("quadgl: elements of 'points` must be ", 485 "numbers only")); 486 } 487 /* We have n-1 intervals and a and b, hence n-1 + 2 results */ 488 results = mat[sp + 1]; 489 if (a != points[0]) { 490 results[0] = quadglcore(a, points[0]); 491 } else { 492 results[0] = 0; 493 } 494 if (sp == 1) { 495 if (b != points[0]) { 496 results[1] = quadglcore(points[0], b); 497 } else { 498 results[1] = 0; 499 } 500 } else { 501 for (k = 1; k < sp; k++) { 502 results[k] = quadglcore(points[k - 1], points[k]); 503 } 504 if (b != points[k - 1]) { 505 results[k] = quadglcore(points[k - 1], b); 506 } else { 507 results[k] = 0; 508 } 509 } 510 } else { 511 if (!isint(points) || points <= 0) 512 return newerror(strcat("quadgl: variable 'points` must be a ", 513 "list or a positive integer")); 514 /* Taking "points" as the number of equally spaced intervals */ 515 results = mat[points + 1]; 516 /* It is easy if a,b lie on the real line */ 517 if (isreal(a) && isreal(b)) { 518 length = abs(a - b); 519 segment = length / points; 520 521 for (k = 1; k <= points; k++) { 522 results[k - 1] = 523 quadglcore(a + (k - 1) * segment, a + k * segment); 524 } 525 } else { 526 /* Other contours by way of a list of points */ 527 slope = (im(b) - im(a)) / (re(b) - re(a)); 528 C = (im(a) + slope) * re(a); 529 length = abs(re(a) - re(b)); 530 segment = length / points; 531 532 /* y = mx+C where m is the slope, x is the real part and y the 533 * imaginary part */ 534 if(re(a)>re(b))swap(a,b); 535 for (k = re(a); k <= (re(b)); k+=segment) { 536 x1 = slope*(k) + C; 537 results[k] = quadglcore(k + x1 * 1i); 538 } 539 } /* else of isreal */ 540 } /* else of ismat|islist */ 541 } /* else of isnull(points) */ 542 /* With a bit of undeserved luck we have a result by now. */ 543 sp = size(results); 544 for (k = 0; k < sp; k++) { 545 sum += results[k]; 546 } 547 return sum; 548} 549 550define quad(a, b, points = -1, method = "tanhsinh") 551{ 552 if (isnull(a) || isnull(b) || param(0) < 2) 553 return newerror("quad: both limits must be given"); 554 if (isstr(a)) { 555 if (strncmp(a, "cinf", 1) == 0) 556 return 557 newerror(strcat 558 ("quad: complex infinity not yet supported, use", 559 " 'pinf' or 'ninf' respectively")); 560 } 561 if (isstr(b)) { 562 if (strncmp(b, "cinf", 1) == 0) 563 return 564 newerror(strcat 565 ("quad: complex infinity not yet supported, use", 566 " 'pinf' or 'ninf' respectively")); 567 } 568 569 if (param(0) == 3) { 570 if (isstr(points)) 571 method = points; 572 } 573 574 if (strncmp(method, "tanhsinh", 1) == 0) { 575 if (!isstr(points)) { 576 if (points == -1) { 577 return quadts(a, b); 578 } else { 579 return quadts(a, b, points); 580 } 581 } else { 582 return quadts(a, b); 583 } 584 } 585 586 if (strncmp(method, "gausslegendre", 1) == 0) { 587 if (!isstr(points)) { 588 if (points == -1) { 589 return quadgl(a, b); 590 } else { 591 return quadgl(a, b, points); 592 } 593 } else { 594 return quadgl(a, b); 595 } 596 } 597} 598 599define makerange(start, end, steps) 600{ 601 local ret k l step C length slope x1 x2 y1 y2; 602 local segment; 603 steps = int (steps); 604 if (steps < 1) { 605 return newerror("makerange: number of steps must be > 0"); 606 } 607 if (!isnum(start) || !isnum(end)) { 608 return newerror("makerange: only numbers are supported yet"); 609 } 610 if (isreal(start) && isreal(end)) { 611 step = (end - start) / (steps); 612 print step; 613 ret = mat[steps + 1]; 614 for (k = 0; k <= steps; k++) { 615 ret[k] = k * step + start; 616 } 617 } else { 618 ret = mat[steps + 1]; 619 if (re(start) > re(end)) { 620 swap(start, end); 621 } 622 623 slope = (im(end) - im(start)) / (re(end) - re(start)); 624 C = im(start) - slope * re(start); 625 length = abs(re(start) - re(end)); 626 segment = length / (steps); 627 628 for (k = re(start), l = 0; k <= (re(end)); k += segment, l++) { 629 x1 = slope * (k) + C; 630 ret[l] = k + x1 * 1i; 631 } 632 633 } 634 return ret; 635} 636 637define makecircle(radius, center, points) 638{ 639 local ret k a b twopi centerx centery; 640 if (!isint(points) || points < 2) { 641 return 642 newerror("makecircle: number of points is not a positive integer"); 643 } 644 if (!isnum(center)) { 645 return newerror("makecircle: center does not lie on the complex plane"); 646 } 647 if (!isreal(radius) || radius <= 0) { 648 return newerror("makecircle: radius is not a real > 0"); 649 } 650 ret = mat[points]; 651 twopi = 2 * pi(); 652 centerx = re(center); 653 centery = im(center); 654 for (k = 0; k < points; k++) { 655 a = centerx + radius * cos(twopi * k / points); 656 b = centery + radius * sin(twopi * k / points); 657 ret[k] = a + b * 1i; 658 } 659 return ret; 660} 661 662define makeellipse(angle, a, b, center, points) 663{ 664 local ret k x y twopi centerx centery; 665 if (!isint(points) || points < 2) { 666 return 667 newerror("makeellipse: number of points is not a positive integer"); 668 } 669 if (!isnum(center)) { 670 return 671 newerror("makeellipse: center does not lie on the complex plane"); 672 } 673 if (!isreal(a) || a <= 0) { 674 return newerror("makecircle: a is not a real > 0"); 675 } 676 if (!isreal(b) || b <= 0) { 677 return newerror("makecircle: b is not a real > 0"); 678 } 679 if (!isreal(angle)) { 680 return newerror("makecircle: angle is not a real"); 681 } 682 ret = mat[points]; 683 twopi = 2 * pi(); 684 centerx = re(center); 685 centery = im(center); 686 for (k = 0; k < points; k++) { 687 x = centerx + a * cos(twopi * k / points) * cos(angle) 688 - b * sin(twopi * k / points) * sin(angle); 689 y = centerx + a * cos(twopi * k / points) * sin(angle) 690 + b * sin(twopi * k / points) * cos(angle); 691 ret[k] = x + y * 1i; 692 } 693 return ret; 694} 695 696define makepoints() 697{ 698 local ret k; 699 ret = mat[param(0)]; 700 for (k = 0; k < param(0); k++) { 701 if (!isnum(param(k + 1))) { 702 return 703 newerror(strcat 704 ("makepoints: parameter number \"", str(k + 1), 705 "\" is not a number")); 706 } 707 ret[k] = param(k + 1); 708 } 709 return ret; 710} 711 712 713config("resource_debug", resource_debug_level),; 714if (config("resource_debug") & 3) { 715 print "quadtsdeletenodes()"; 716 print "quadtscomputenodes(order, expo, eps)"; 717 print "quadtscore(a,b,n)"; 718 print "quadts(a,b,points)"; 719 print "quadglcomputenodes(N)"; 720 print "quadgldeletenodes()"; 721 print "quadglcore(a,b,n)"; 722 print "quadgl(a,b,points)"; 723 print "quad(a,b,points=-1,method=\"tanhsinh\")"; 724 print "makerange(start, end, steps)"; 725 print "makecircle(radius, center, points)"; 726 print "makeellipse(angle, a, b, center, points)"; 727 print "makepoints(a1,[...])"; 728} 729