sult(x,y,z,a) /* tests correctness of a = appr(x,y,z) */ { local r, n, s, v; if (y == 0) return (a != x); r = x - a; n = a/y; if (!isint(n)) return 1; if (abs(r) >= abs(y)) return 1; if (r == 0) return 0; if (z & 16) { if (abs(r) > abs(y)/2) return 1; if (abs(r) < abs(y)/2) return 0; z &= 15; } s = sgn(r); switch (z) { case 0: v = (s == sgn(y)); break; case 1: v = (s == -sgn(y)); break; case 2: v = (s == sgn(x)); break; case 3: v = (s == -sgn(x)); break; case 4: v = (s > 0); break; case 5: v = (s < 0); break; case 6: v = (s == sgn(x/y)); break; case 7: v = (s == -sgn(x/y)); break; case 8: v = iseven(n); break; case 9: v = isodd(n); break; case 10: v = (x/y > 0) ? iseven(n) : isodd(n); break; case 11: v = (x/y > 0) ? isodd(n) : iseven(n); break; case 12: v = (y > 0) ? iseven(n) : isodd(n); break; case 13: v = (y > 0) ? isodd(n) : iseven(n); break; case 14: v = (x > 0) ? iseven(n) : isodd(n); break; case 15: v = (x > 0) ? isodd(n) : iseven(n); breÒÒak; } return !v; } /* * test2600 - perform all of the above tests a bunch of times */ define test2600(verbose, tnum) { local n; /* test parameter */ local ep; /* test parameter */ local i; /* set test parameters */ n = 5; /* internal test loop count */ if (isnull(verbose)) { verbose = defaultverbose; } if (isnull(tnum)) { tnum = 1; /* initial test number */ } /* * test a lot of stuff */ srand(2600e2600); ep = 1e-250; err += testismult(strcat(str(tnum++), ": mult"), n*20, verbose); err += testappr(strcat(str(tnum++), ": appr"), n*40, verbose); err += testexp(strcat(str(tnum++),": exp"), n, ep, verbose); err += testln(strcat(str(tnum++),": ln"), n, ep, verbose); err += testpower(strcat(str(tnum++),": power"), n, rand(2,10), ep, verbose); err += testgcd(strcat(str(tnum++),": gcd"), n, ep, verbose); for (i=0; i < 32; ++i) { config("sqrt", i); err += testsqrt(strcat(str(tnum++),": sqrt",str(i)), n*10, ep, verbose); } err += testpower2(strcat(str(tnum++),": power"), n*4, ep, verbose); if (verbose > 1) { if (err) { print "***", err, "error(s) found in test2600"; } else { print "no errors in test2600"; } } return tnum; } æbindingsç randrun.calè set8700.calé mersenne.calê chrem.calë dotest.calì intfile.calípix.calî test2300.calï test3400.calðalg_config.calñ test4100.calò test1700.calóchi.calôrandmprime.calõbeer.calöpoly.cal÷¸ hello.cal# bindings - default key bindings for calc line editing functions # # Copyright (C) 1999 David I. Bell # # 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. # # @(#) $Revision: 30.1 $ # @(#) $Id: bindings,v 30.1 2007/03/16 11:09:54 chongo Exp $ # @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/bindings,v $ # # Under source code control: 1993/05/02 20:09:19 # File existed as early as: 1993 # # Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ # NOTE: This facility is ignored if calc was compiled with GNU-readline. # In that case, the standard readline mechanisms (see readline(3)) # are used in place of those found below. map base-map default insert-char ^@ set-mark ^A start-of-line ^B backward-char ^D delete-char ^E end-of-line ^F forward-char ^H backward-kill-char ^J new-line ^K kill-line ^L refresh-line ^M new-line ^N forward-history ^O save-line ^P backward-history ^R reverse-search ^T swap-chars ^U flush-input ^V quote-char ^W kill-region ^Y yank ^? backward-kill-char ^[ ignore-char esc-map map esc-map default ignore-char base-map G start-of-line H backward-history P forward-history K backward-char M forward-char O end-of-line S delete-char g goto-line s backward-word t forward-word d forward-kill-word u uppercase-word l lowercase-word h list-history ^[ flush-input [ arrow-key /* * randrun - perform a run test on rand() * * Copyright (C) 1999 David I. Bell * * 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. * * @(#) $Revision: 30.1 $ * @(#) $Id: randrun.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $ * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/randrun.cal,v $ * * Under source code control: 1995/02/12 20:00:06 * File existed as early as: 1995 * * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ */ /* * If X(j) < X(j+1) < ... X(j+k) >= X(j+k+1), then we have a run of 'k'. * We ignore the run breaker, X(j+k+1), and start with X(j+k+2) when * considering a new run in order to make our runs chi independent. * * See Knuth's "Art of Computer Programming - 2nd edition", * Volume 2 ("Seminumerical Algorithms"), Section 3.3.2. * "G. Run test", pp. 65-68, * "problem #14", pp. 74, 536. * * We use the suggestion in problem #14 to allow an application of the * chi-square test and to make estimating the run length probs easy. */ define randrun(run_cnt) { local i; /* index */ local max_run; /* longest run */ local long_run_cnt; /* number of runs longer than MAX_RUN */ local run; /* current run length */ local tally_sum; /* sum of all tally values */ local last; /* last random number */ local current; /* current random number */ local MAX_RUN = 9; /* max run we will keep track of */ local mat tally[1:MAX_RUN]; /* tally of length of a rise run of 'x' */ local mat prob[1:MAX_RUN]; /* prob[x] = probability of 'x' length run */ /* * parse args */ if (param(0) == 0) { run_cnt = 65536; } /* * run setup */ max_run = 0; /* no runs yet */ long_run_cnt = 0; /* no long runs set */ current = rand(); /* our first number */ run = 1; /* * compute the run length probabilities * * A run length of 'r' occurs with a probability of: * * 1/r! - 1/(r+1)! */ for (i=1; i <= MAX_RUN; ++i) { prob[i] = 1.0/fact(i) - 1.0/fact(i+1); } /* * look at a number of random number trials */ for (i=0; i < run_cnt; ++i) { /* get our current number */ last = current; current = rand(); /* look for a run break */ if (current < last) { /* record the stats */ if (run > max_run) { max_run = run; } if (run > MAX_RUN) { ++long_run_cnt; } else { ++tally[run]; } /* start a new run */ current = rand(); run = 1; /* note the continuing run */ } else { ++run; } } /* determine the number of runs found */ tally_sum = matsum(tally) + long_run_cnt; /* * print the stats */ printf("rand run test used %d values to produce %d runs\n", run_cnt, tally_sum); for (i=1; i <= MAX_RUN; ++i) { printf("length=%d\tprob=%9.7f\texpect=%d \tcount=%d \terr=%9.7f\n", i, prob[i], round(tally_sum*prob[i]), tally[i], (tally[i] - round(tally_sum*prob[i]))/tally_sum); } printf("length>%d\t\t\t\t\tcount=%d\n", MAX_RUN, long_run_cnt); printf("max length=%d\n", max_run); } if (config("resource_debug") & 3) { print "randrun([run_length]) defined"; } /* * set8700 - environment for dotest line tests for the 8700 set of regress.cal * * Copyright (C) 2006 Ernest Bowen and Landon Curt Noll * * 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. * * @(#) $Revision: 30.1 $ * @(#) $Id: set8700.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $ * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/set8700.cal,v $ * * Under source code control: 2006/05/20 14:10:11 * File existed as early as: 2006 * * chongo /\oo/\ http://www.isthe.com/chongo/ * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ */ /* * setup global variables for dotest() to use with set8700.set */ global set8700_A; global set8700_B; global set8700_M; global set8700_M1; global set8700_M2; global set8700_L; global set8700_L1; global set8700_L2; global set8700_O; global set8700_P; global set8700_P1; global set8700_P2; global set8700_Q; global set8700_R; global set8700_S; global set8700_X; global set8700_Y; global set8700_x; global set8700_y; define set8700_getA1() = set8700_A; define set8700_getA2() { return set8700_A; } define set8700_getvar() {local a = 42; protect(a,256); return a;} define set8700_f(set8700_x) = set8700_x^2; define set8700_g(set8700_x) { if (isodd(set8700_x)) protect(set8700_x, 256); return set8700_x; } obj set8700_point { set8700_x, set8700_y, set8700_z } global mat set8700_c[] = { 1, 2+3i, -5+4i, 5i+6, -7i }; global mat set8700_e[] = { 0, 1, 0, 0, 2, -3/2, 2, -1/2, -3, 0.5, -1.0, 0.5, 1.0, 0.0, 0.0, 0.0 }; /* * mersenne - perform a primality test of 2^p-1, for prime p>1 * * Copyright (C) 1999 David I. Bell and Landon Curt Noll * * Primary author: David I. Bell * * 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. * * @(#) $Revision: 30.1 $ * @(#) $Id: mersenne.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $ * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/mersenne.cal,v $ * * Under source code control: 1991/05/22 21:56:36 * File existed as early as: 1991 * * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ */ /* * NOTE: See lucas.cal for a more general routine. */ define mersenne(p) { local u, i, p_mask; /* firewall */ if (! isint(p)) quit "p is not an integer"; /* two is a special case */ if (p == 2) return 1; /* if p is not prime, then 2^p-1 is not prime */ if (! ptest(p,1)) return 0; /* lltest: u(i+1) = u(i)^2 - 2 mod 2^p-1 */ u = 4; for (i = 2; i < p; ++i) { u = hnrmod(u^2 - 2, 1, p, -1); } /* 2^p-1 is prime iff u(p) = 0 mod 2^p-1 */ return (u == 0); } /* * chrem - chinese remainder theorem/problem solver * * 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. * * @(#) $Revision: 30.1 $ * @(#) $Id: chrem.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $ * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/chrem.cal,v $ * * Under source code control: 1992/09/26 01:00:47 * File existed as early as: 1992 * * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ */ /* * When possible, chrem finds solutions for x of a set of congruences * of the form: * * x = r1 (mod m1) * x = r2 (mod m2) * ... * * where the residues r1, r2, ... and the moduli m1, m2, ... are * given integers. The Chinese remainder theorem states that if * m1, m2, ... are relatively prime in pairs, the above congruences * have a unique solution modulo m1 * m2 * ... If m1, m2, ... * are not relatively prime in pairs, it is possible that no solution * exists. If solutions exist, the general solution is expressible as: * * x = r (mod m) * * where m = lcm(m1,m2,...), and if m > 0, 0 <= r < m. This * solution may be interpreted as: * * x = r + k * m [[NOTE 1]] * * where k is an arbitrary integer. * *** * * usage: * * chrem(r1,m1 [,r2,m2, ...]) * * r1, r2, ... remainder integers or null values * m1, m2, ... moduli integers * * chrem(r_list, [m_list]) * * r_list list (r1,r2, ...) * m_list list (m1,m2, ...) * * If m_list is omitted, then 'defaultmlist' is used. * This default list is a global value that may be changed * by the user. Initially it is the first 8 primes. * * If a remainder is null(), then the corresponding congruence is * ignored. This is useful when working with a fixed list of moduli. * * If there are more remainders than moduli, then the later moduli are * ignored. * * The moduli may be any integers, not necessarily relatively prime in * pairs (as required for the Chinese remainder theorem). Any moduli * may be zero; x = r (mod 0) has the meaning of x = r. * * returns: * * If args were integer pairs: * * r ('r' is defined above, see [[NOTE 1]]) * * If 1 or 2 list args were given: * * (r, m) ('r' and 'm' are defined above, see [[NOTE 1]]) * * NOTE: In all cases, null() is returned if there is no solution. * *** * * This function may be used to solve the following historical problems: * * Sun-Tsu, 1st century A.D. * * To find a number for which the reminders after division by 3, 5, 7 * are 2, 3, 2, respectively: * * chrem(2,3,3,5,2,7) ---> 23 * * Fibonacci, 13th century A.D. * * To find a number divisible by 7 which leaves remainder 1 when * divided by 2, 3, 4, 5, or 6: * * * chrem(list(0,1,1,1,1,1),list(7,2,3,4,5,6)) ---> (301,420) * * i.e., any value that is 301 mod 420. */ static defaultmlist = list(2,3,5,7,11,13,17,19); /* The first eight primes */ define chrem() { local argc; /* number of args given */ local rlist; /* reminder list - ri */ local mlist; /* modulus list - mi */ local list_args; /* true => args given are lists, not r1,m1, ... */ local m,z,r,y,d,t,x,u,i; /* * parse args */ argc = param(0); if (argc == 0) { quit "usage: chrem(r1,m1 [,r2,m2 ...]) or chrem(r_list, m_list)"; } list_args = islist(param(1)); if (list_args) { rlist = param(1); mlist = (argc == 1) ? defaultmlist : param(2); if (size(rlist) > size(mlist)) { quit "too many residues"; } } else { if (argc % 2 == 1) { quit "odd number integers given"; } rlist = list(); mlist = list(); for (i=1; i <= argc; i+=2) { push(rlist, param(i)); push(mlist, param(i+1)); } } /* * solve the problem found in rlist & mlist */ m = 1; z = 0; while (size(rlist)) { r=pop(rlist); y=abs(pop(mlist)); if (r==null()) continue; if (m) { if (y) { d = t = z - r; m = lcm(x=m, y); while (d % y) { u = x; x %= y; swap(x,y); if (y==0) return; z += (t *= -u/y); } } else { if ((r % m) != (z % m)) return; else { m = 0; z = r; } } } else if (((y) && (r % y != z % y)) || (r != z)) return; } if (m) { z %= m; if (z < 0) z += m; } /* * return information as required */ if (list_args) { return list(z,m); } else { return z; } } if (config("resource_debug") & 3) { print "chrem(r1,m1 [,r2,m2 ...]) defined"; print "chrem(rlist [,mlist]) defined"; } /* * dotest - test truth statements found in line tests of dotest_testline file * * This file was created by Ernest Bowen * and modified by Landon Curt Noll. * * This dotest_code has been placed in the public domain. Please do not * copyright this dotest_code. * * ERNEST BOWEN AND LANDON CURT NOLL DISCLAIMS ALL WARRANTIES WITH REGARD TO * THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MER- * CHANTABILITY AND FITNESS. IN NO EVENT SHALL LANDON CURT * NOLL BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL * DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF * USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. * * @(#) $Revision: 30.2 $ * @(#) $Id: dotest.cal,v 30.2 2007/03/16 11:09:54 chongo Exp $ * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/dotest.cal,v $ * * This file is not covered under version 2.1 of the GNU LGPL. * * Under source dotest_code control: 2006/03/08 05:54:09 * File existed as early as: 2006 */ /* * dotest - perform tests from dotest_testline file * * given: * dotest_file filename containing single test lines * dotest_code regress.cal test number to use (def: 0) * dotest_maxcond max error conditions allowed (def: <0 ==> 2^31-1) * * returns: * number of line test failures * * NOTE: All variables used by the dotest() function start with "dotest_". * The dotest_file and dotest_read should not use any variable * that starts with "dotest_". */ define dotest(dotest_file, dotest_code = 0, dotest_maxcond = -1) { local dotest_f_file; /* open file containing test lines */ local dotest_testline; /* test line */ local dotest_testeval; /* eval value from dotest_testline test line */ local dotest_tmperrcnt; /* temp error count after line test */ local dotest_errcnt; /* total number of errors */ local dotest_failcnt; /* number of line tests failed */ local dotest_testnum; /* number of test lines evaluated */ local dotest_linenum; /* test line number */ local dotest_old_errmax; /* value of errmax() prior to calling */ local dotest_old_errcount; /* value of errcount() prior to calling */ /* * preserve calling stats */ dotest_old_errmax = errmax(); dotest_old_errcount = errcount(0); /* * initialize test accounting */ dotest_errcnt = errcount(); dotest_failcnt = 0; dotest_testnum = 0; dotest_linenum = 0; /* * setup error accounting for dotest */ if (dotest_maxcond >= 0 && dotest_maxcond < 2147483647) { errmax(dotest_maxcond + dotest_old_errcount + 1),; } else { errmax(2147483647),; } /* * open the test line file */ printf("%d-: opening line file: %d", dotest_code, dotest_file); dotest_f_file = fpathopen(dotest_file, "r"); if (!isfile(dotest_f_file)) { printf("**** Unable to file or open file \"%s\"\n", dotest_file); quit; } printf('%d: testing "%s"\n', dotest_code, dotest_file); /* * perform dotest_testline test on each line of the file */ for (;;) { /* get the next test line */ dotest_testline = fgets(dotest_f_file); ++dotest_linenum; if (iserror(dotest_testline)) { quit "**** Error while reading file"; } else if (isnull(dotest_testline)) { /* EOF - end of test file */ break; } /* skip empty lines */ if (dotest_testline == "\n") { continue; } /* evaluate the test line */ dotest_testeval = eval(dotest_testline); /* ignore white space or comment lines */ if (isnull(dotest_testeval)) { continue; } /* look for test line parse errors */ if (iserror(dotest_testeval)) { printf("**** evaluation error: "); ++dotest_failcnt; /* look for test line dotest_failcnt */ } else if (dotest_testeval != 1) { printf("**** did not return 1: "); ++dotest_failcnt; } /* show the test line we just performed */ printf("%d-%d: %s", dotest_code, dotest_linenum, dotest_testline); /* error accounting */ dotest_tmperrcnt = errcount() - dotest_errcnt; if (dotest_tmperrcnt > 0) { /* report any other errors */ if (dotest_tmperrcnt > 1) { printf("%d-%d: NOTE: %d error conditions(s): %s\n", dotest_code, dotest_linenum, dotest_tmperrcnt); } /* report the calc error string */ printf("%d-%d: NOTE: last error string: %s\n", dotest_code, dotest_linenum, strerror()); /* new error count level */ dotest_errcnt = errcount(); if (dotest_maxcond >= 0 && dotest_old_errcount-dotest_errcnt > dotest_maxcond) { printf("%d-%d: total error conditions: %d > %d\n", dotest_code, dotest_linenum, dotest_maxcond, dotest_old_errcount-dotest_errcnt); } } } /* * test the close of the line file */ printf("%d-: detected %d error condition(s), many of which may be OK\n", dotest_code, dotest_old_errcount-dotest_errcnt); printf("%d-: closing line file: %d\n", dotest_code, dotest_file); fclose(dotest_f_file); /* * test line file accounting */ if (dotest_failcnt > 0) { printf("**** %d-: %d test failure(s) in %d line(s)\n", dotest_code, dotest_failcnt, dotest_linenum); } else { printf("%d-: no failure(s) in %d line(s)\n", dotest_code, dotest_linenum); } /* * preppare to return to the caller environment * * We increase the caller's error count by the number * of line tests that failed, not the number of internal * errors that were noted. */ errmax(dotest_old_errmax),; errcount(dotest_old_errcount + dotest_failcnt),; /* * All Done!!! -- Jessica Noll, Age 2 */ return dotest_failcnt; } /* * intfile - integer to file and file to integer conversion * * Copyright (C) 2001 Landon Curt Noll * * 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. * * @(#) $Revision: 30.1 $ * @(#) $Id: intfile.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $ * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/intfile.cal,v $ * * Under source code control: 2001/03/31 08:13:11 * File existed as early as: 2001 * * chongo /\oo/\ http://www.isthe.com/chongo/ * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ */ /* * NOTE: Because leading HALF values are trimmed from integer, a file * that begins with lots of 0 bits (in the case of big endian) * or that ends with lots of 0 bits (in the case of little endian) * will be changed when the subsequent integer is written back. */ /* * file2be - convert a file into an big endian integer * * given: * filename filename to read * * returns: * integer read from its contents on big endian order */ define file2be(filename) { local fd; /* open file */ local ret; /* integer to return */ local c; /* character read from the file */ local i; /* * open the file for reading */ fd = fopen(filename, "rb"); if (!isfile(fd)) quit "file2be: cannot open file for reading"; /* * read the contents of the file * * The first octets become the most significant bits of the integer. */ ret = 0; while (! isnull(c = fgetc(fd))) { ret <<= 8; ret += ord(c); } /* * cleanup and return the integer */ fclose(fd); return ret; } /* * file2le - convert a file into an little endian integer * * given: * filename filename to read * * returns: * integer read from its contents on little endian order */ define file2le(filename) { local fd; /* open file */ local ret; /* integer to return */ local c; /* character read from the file */ local shft; /* bit shift for the c value */ local i; /* * open the file for reading */ fd = fopen(filename, "rb"); if (!isfile(fd)) quit "file2le: cannot open file for reading"; /* * read the contents of the file into a string * * The first octets become are the least significant bits of the integer. */ ret = 0; shft = 0; while (! isnull(c = fgetc(fd))) { ret |= (ord(c) << shft); shft += 8; } /* * cleanup and return the integer */ fclose(fd); return ret; } /* * be2file - convert a big endian integer into a file * * given: * v integer to write to the file * filename filename to write * * returns: * The number of octets written to the file. * * NOTE: The absolute value of the integer is written to the file. */ define be2file(v, filename) { local fd; /* open file */ local octlen; /* length of v in octets */ local i; /* * firewall */ if (!isint(v)) { quit "be2file: 1st arg not an integer"; } v = abs(v); /* * open the file for writing */ fd = fopen(filename, "wb"); if (!isfile(fd)) quit "be2file: cannot open file for writing"; /* * write the octets to the file * * The most significant bits of the integer become the first file octets. */ octlen = int((highbit(v)+8) / 8); for (i=octlen-1; i >= 0; --i) { fputc(fd, char(v >> (i*8))); } /* * cleanup */ fclose(fd); return octlen; } /* * le2file - convert a little endian integer into a file * * given: * v integer to write to the file * filename filename to write * * returns: * The number of octets written to the file. * * NOTE: The absolute value of the integer is written to the file. */ define le2file(v, filename) { local fd; /* open file */ local cnt; /* octets written */ /* * firewall */ if (!isint(v)) { quit "be2file: 1st arg not an integer"; } v = abs(v); /* * open the file for writing */ fd = fopen(filename, "wb"); if (!isfile(fd)) quit "le2file: cannot open file for writing"; /* * Write the octets to the file. * * The least significant bits of the integer become the first file octets. */ cnt = 0; while (v > 0) { fputc(fd, char(v)); v >>= 8; ++cnt; } /* * cleanup */ fclose(fd); return cnt; } /* * pix - iterative method of finding the number of primes less than x * * Copyright (C) 1999 Landon Curt Noll * * 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. * * @(#) $Revision: 30.1 $ * @(#) $Id: pix.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $ * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/pix.cal,v $ * * Under source code control: 1996/07/09 03:14:14 * File existed as early as: 1996 * * chongo /\oo/\ http://www.isthe.com/chongo/ * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ */ /* * Here is an iterative method of finding the number of primes less than * or equal to a given number. This method is from "Computer Recreations" * June 1996 issue of Scientific American. * * NOTE: For reasonable values of x, the builtin function pix(x) is * much faster. This code is provided because the method * is interesting. */ define pi_of_x(x) { local An; /* A(n) */ local An1; /* A(n-1) */ local An2; /* A(n-2) */ local An3; /* A(n-3) */ local primes; /* number of primes found */ local n; /* loop counter */ /* * setup */ An1 = 2; An2 = 0; An3 = 3; primes = 1; /* * main A(n+1)=A(n-1)+A(n-2) sequence loop */ for (n = 3; n < x; ++n) { An = An2 + An3; An3 = An2; An2 = An1; An1 = An; if (An % n == 0) ++primes; } return primes; } /* * test2300 - 2300 series of the regress.cal test suite * * Copyright (C) 1999 Landon Curt Noll * * 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. * * @(#) $Revision: 30.1 $ * @(#) $Id: test2300.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $ * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test2300.cal,v $ * * Under source code control: 1995/07/09 06:12:13 * File existed as early as: 1995 * * chongo /\oo/\ http://www.isthe.com/chongo/ * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ */ obj matrix {m} /* * matrix_inc - increment the matrix inside the object */ define matrix_inc(a) { local i; /* increment each matrix member */ for (i= 0; i < size(a.m); i++) ++a.m[[i]]; return a; } /* * matrix_dec - decrement the matrix inside the object */ define matrix_dec(a) { local i; /* decrement each matrix member */ for (i= 0; i < size(a.m); i++) --a.m[[i]]; return a; } /* * mkmat - load the matrix inside the object */ define mkmat() { local s, M, i, v; /* firewall */ s = param(0); if (s == 0) quit "Need at least one argument"; /* create the matrix */ mat M[s]; /* load the matrix with the args */ for (i = 0; i < s; i++) M[i] = param(i + 1); /* create the object with the matrix */ obj matrix v; v.m = M; return v; } /* * ckmat - check if the matrix inside an object has a set of given values */ define ckmat() { local s, a, i; /* firewall */ s = param(0); if (s < 2) quit "Need at least two arguments"; /* get the object to test */ a = param(1); /* verify the matrix in the object is the right size */ if (size(a.m) != s-1) { return 0; } /* check each matrix element with the args passed */ for (i = 2; i <= s; i++) { if (a.m[i-2] != param(i)) { /* args do not match */ return 0; } } /* args match the matrix in the object */ return 1; } /* * test3400 - 3400 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. * * @(#) $Revision: 30.1 $ * @(#) $Id: test3400.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $ * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test3400.cal,v $ * * Under source code control: 1995/12/02 05:20:11 * File existed as early as: 1995 * * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ */ /* * tests of performance of some trigonometric functions * * test3401 tests abs(acot(cot(x)) - x) <= eps for x = k * eps < pi * test3402 tests abs(tan(x/2) - csc(x) + cot(x)) <= eps * test3403 tests abs(tan(x) - cot(x) + 2 * cot(2 * x)) <= eps * test3404 tests abs(cot(x/2) - csc(x) - cot(x)) <= eps * test3405 tests atan(tan(x)) == x for x = k * eps, abs(x) <= pi/2 * test3406 tests abs(sec(x) - sec(x + 2 * N * pi)) <= eps * * To run say, test1 n times give instruction test1(n, eps); eps * defaults to epsilon(). * * Here pi1k is pi to 1000 decimal places; x is a random real number * except when x is described as k * eps, in which case k is a random * integer such that x is in the specified range. * * In the last test N is a large random integer, but it is assumed * that eps is large compared with N * 1e-1000. * * I am surprised that test3406 seems to give no errors - I had expected * that the two sides might differ by eps. [[test changed to test eps error]] */ defaultverbose = 1; /* default verbose value */ global pi1k = pi(1e-1000); define test3401(str, n, eps, verbose) { local i, m, x, y, N; if (isnull(verbose)) verbose = defaultverbose; if (verbose > 0) { print str:":",:; } if (isnull(n)) n = 250; if (isnull(eps)) eps = epsilon(); m = 0; N = pi(eps)/eps; for (i = 0; i < n; i++) { x = rand(1, N) * eps; y = cot(x, eps); if (verbose > 1) printf("%r\n", x); if (abs(acot(y, eps) - x) > eps) { if (verbose > 1) { printf("*** Failure for x = %r\n", x); } m++; } } if (verbose > 0) { if (m) { printf("*** %d error(s)\n", m); } else { printf("no errors\n"); } } return m; } define test3402(str, n, eps, verbose) { local i, m, x, y, N; if (isnull(verbose)) verbose = defaultverbose; if (verbose > 0) { print str:":",:; } if (isnull(n)) n = 250; if (isnull(eps)) eps = epsilon(); eps = abs(eps); m = 0; N = 1e10; for (i = 0; i < n; i++) { x = rand(-N, N)/rand(1, N); y = tan(x/2, eps) - csc(x,eps) + cot(x,eps); if (verbose > 1) printf("%r\n", x); if (abs(y) > eps) { if (verbose > 1) { printf("*** Failure for x = %r\n", x); } m++; } } if (verbose > 0) { if (m) { printf("*** %d error(s)\n", m); } else { printf("no errors\n"); } } return m; } define test3403(str, n, eps, verbose) { local i, m, x, y, N; if (isnull(verbose)) verbose = defaultverbose; if (verbose > 0) { print str:":",:; } if (isnull(n)) n = 250; if (isnull(eps)) eps = epsilon(); eps = abs(eps); m = 0; N = 1e10; for (i = 0; i < n; i++) { x = rand(-N, N)/rand(1, N); y = tan(x, eps) - cot(x,eps) + 2 * cot(2 * x,eps); if (verbose > 1) printf("%r\n", x); if (abs(y) > eps) { m++; if (verbose > 1) { printf("*** Failure for x = %r\n", x); } } } if (verbose > 0) { if (m) { printf("*** %d error(s)\n", m); } else { printf("no errors\n"); } } return m; } define test3404(str, n, eps, verbose) { local i, m, x, y, N; if (isnull(verbose)) verbose = defaultverbose; if (verbose > 0) { print str:":",:; } if (isnull(n)) n = 250; if (isnull(eps)) eps = epsilon(); eps = abs(eps); m = 0; N = 1e10; for (i = 0; i < n; i++) { x = rand(-N, N)/rand(1, N); y = cot(x/2, eps) - csc(x,eps) - cot(x,eps); if (verbose > 1) printf("%r\n", x); if (abs(y) > eps) { m++; if (verbose > 1) { printf("*** Failure for x = %r\n", x); } } } if (verbose > 0) { if (m) { printf("*** %d error(s)\n", m); } else { printf("no errors\n"); } } return m; } define test3405(str, n, eps, verbose) { local i, m, x, y, N; if (isnull(verbose)) verbose = defaultverbose; if (verbose > 0) { print str:":",:; } if (isnull(n)) n = 250; if (isnull(eps)) eps = epsilon(); m = 0; N = pi(eps)/eps; N = quo(N, 2, 0); for (i = 0; i < n; i++) { x = rand(-N, N) * eps; y = tan(x, eps); if (verbose > 1) printf("%r\n", x); if (atan(y, eps) != x) { m++; if (verbose > 1) { printf("*** Failure for x = %r\n", x); } } } if (verbose > 0) { if (m) { printf("*** %d error(s)\n", m); } else { printf("no errors\n"); } } return m; } define test3406(str, n, eps, verbose) { local i, m, x, y, z, N; if (isnull(verbose)) verbose = defaultverbose; if (verbose > 0) { print str:":",:; } if (isnull(n)) n = 250; if (isnull(eps)) eps = epsilon(); m = 0; for (i = 0; i < n; i++) { x = rand(-1e10, 1e10)/rand(1, 1e10); N = rand(-1e10, 1e10); y = sec(x, eps); z = sec(x + 2 * N * pi1k, eps); if (verbose > 1) printf("%r, %d\n", x, N); if (abs(y-z) > eps) { m++; if (verbose > 1) { printf("*** Failure for x = %r\n", x); } } } if (verbose > 0) { if (m) { printf("*** %d error(s)\n", m); } else { printf("no errors\n"); } } return m; } /* * test3400 - perform all of the above tests */ define test3400(verbose, tnum) { local n; /* test parameter */ local eps; /* test parameter */ local i; /* * set test parameters */ if (isnull(verbose)) { verbose = defaultverbose; } n = 250; eps = epsilon(); srand(3400e3400); /* * test a lot of stuff */ err += test3401(strcat(str(tnum++), \ ": acot(cot(x))"), n, eps, verbose); err += test3402(strcat(str(tnum++), \ ": tan(x/2)-csc(x)+cot(x)"), n, eps, verbose); err += test3403(strcat(str(tnum++), \ ": tan(x)-cot(x)+2*cot(2*x)"), n, eps, verbose); err += test3404(strcat(str(tnum++), \ ": cot(x/2)-csc(x)-cot(x)"), n, eps, verbose); err += test3405(strcat(str(tnum++), \ ": atan(tan(x))"), n, eps, verbose); err += test3406(strcat(str(tnum++), \ ": sec(x)-sec(x+2*N*pi)"), n, eps, verbose); /* * test results */ if (verbose > 1) { if (err) { print "***", err, "error(s) found in test3400"; } else { print "no errors in test3400"; } } return tnum; } /* * alg_config - help determine optimal values for algorithm levels * * Copyright (C) 2006 Landon Curt Noll * * 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. * * @(#) $Revision: 30.1 $ * @(#) $Id: alg_config.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $ * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/alg_config.cal,v $ * * Under source code control: 2006/06/07 14:10:11 * File existed as early as: 2006 * * chongo /\oo/\ http://www.isthe.com/chongo/ * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ */ global test_time; /* try for this many seconds in loop test */ /* * mul_loop - measure the CPU time to perform a set of multiply loops * * given: * repeat number of multiply loops to perform * x array of 5 values, each the same length in BASEB-bit words * * NOTE: When their lengths are 1 BASEB-bit word, then a * dummy loop of simple constants are used. Thus the * length == 1 is an approximation of loop overhead. * * returns: * approximate runtime to perform repeat the multiply loops * * NOTE: This is an internal support function that is normally * not called directly from the command line. Call the * function best_mul2() instead. */ define mul_loop(repeat, x) { local start; /* start of execution */ local end; /* end of execution */ local answer; /* multiplicand */ local len; /* length of each element */ local baseb_bytes; /* bytes in a BASEB-bit word */ local i; /* firewall */ if (!isint(repeat) || repeat < 0) { quit "mul_loop: 1st arg: repeat must be an integer > 0"; } if (size(*x) != 5) { quit "mul_loop: 2nd arg matrix does not have 5 elements"; } if (matdim(*x) != 1) { quit "mul_loop: 2nd arg matrix is not 1 dimensional"; } if (matmin(*x, 1) != 0) { quit "mul_loop: 2nd arg matrix index range does not start with 0"; } if (matmax(*x, 1) != 4) { quit "mul_loop: 2nd arg matrix index range does not end with 4"; } baseb_bytes = config("baseb") / 8; len = sizeof((*x)[0]) / baseb_bytes; for (i=1; i < 4; ++i) { if ((sizeof((*x)[i]) / baseb_bytes) != len) { quit "mul_loop: 2nd arg matrix elements are not of equal BASEB-bit word length"; } } /* multiply pairwise, all sets of a given length */ start = usertime(); for (i=0; i < repeat; ++i) { if (len == 1) { /* we use len == 1 to test this tester loop overhead */ answer = 0 * 0; answer = 0 * 0; answer = 0 * 0; answer = 0 * 0; /**/ answer = 0 * 0; answer = 0 * 0; answer = 0 * 0; answer = 0 * 0; /**/ answer = 0 * 0; answer = 0 * 0; answer = 0 * 0; answer = 0 * 0; /**/ answer = 0 * 0; answer = 0 * 0; answer = 0 * 0; answer = 0 * 0; /**/ answer = 0 * 0; answer = 0 * 0; answer = 0 * 0; answer = 0 * 0; } else { answer = (*x)[0] * (*x)[1]; answer = (*x)[0] * (*x)[2]; answer = (*x)[0] * (*x)[3]; answer = (*x)[0] * (*x)[4]; /**/ answer = (*x)[1] * (*x)[0]; answer = (*x)[1] * (*x)[2]; answer = (*x)[1] * (*x)[3]; answer = (*x)[1] * (*x)[4]; /**/ answer = (*x)[2] * (*x)[0]; answer = (*x)[2] * (*x)[1]; answer = (*x)[2] * (*x)[3]; answer = (*x)[2] * (*x)[4]; /**/ answer = (*x)[3] * (*x)[0]; answer = (*x)[3] * (*x)[1]; answer = (*x)[3] * (*x)[2]; answer = (*x)[3] * (*x)[4]; /**/ answer = (*x)[4] * (*x)[0]; answer = (*x)[4] * (*x)[1]; answer = (*x)[4] * (*x)[2]; answer = (*x)[4] * (*x)[3]; } } /* * return duration */ end = usertime(); return end-start; } /* * mul_ratio - ratio of rates of 1st and 2nd multiply algorithms * * given: * len length in BASEB-bit words to multiply * * return: * ratio of (1st / 2nd) algorithm rate * * When want to determine a rate to a precision of 1 part in 1000. * Most systems today return CPU time to at least 10 msec precision. * So to get rates to that precision, we need to time loops to at * least 1000 times as long as the precision (10 msec * 1000) * which usually requires timing of loops that last 10 seconds or more. * * NOTE: This is an internal support function that is normally * not called directly from the command line. Call the * function best_mul2() instead. */ define mul_ratio(len) { local mat x[5]; /* array of values for mul_loop to multiply */ local mat one[5]; /* array if single BASEB-bit values */ local baseb; /* calc word size in bits */ local orig_cfg; /* caller configuration */ local loops; /* number of multiply loops to time */ local tlen; /* time to perform some number of loops */ local tover; /* est of time for loop overhead */ local alg1_rate; /* loop rate of 1st algorithm */ local alg2_rate; /* loop rate of 2nd algorithm */ local i; /* * firewall */ if (!isint(len) || len < 2) { quit "mul_ratio: 1st arg: len is not an integer > 1"; } /* * remember the caller's config state */ orig_cfg = config("all"); config("mul2", 0),; config("sq2", 0),; config("pow2", 0),; config("redc2", 0),; config("tilde", 0),; /* * initialize x, the values we will multiply * * We want these tests to be repeatable as possible, so we will seed * the PRNG in a deterministic way. */ baseb = config("baseb"); srand(sha1(sha1(baseb, config("version")))); for (i=0; i < 5; ++i) { /* force the values to be a full len words long */ x[i] = ((1<<(((len-1) * baseb) + baseb-1)) | randbit(((len-1) * baseb) + baseb-2)); /* single BASEB-bit values */ one[i] = 1; } /* * determine the number of loops needed to test 1st alg */ config("mul2", 2^31-1),; loops = 1/2; do { loops *= 2; tlen = mul_loop(loops, &x); if (config("user_debug") > 3) { printf("\t alg1 loops %d took %.3f sec\n", loops, tlen); } } while (tlen < 1.0); /* * determine the 1st algorithm rate */ loops = max(1, ceil(loops * test_time / tlen)); if (loops < 8) { if (config("user_debug") > 1) { printf(" we must expand loop test time to more than %d secs\n", ceil(test_time * (8 / loops))); } loops = 8; } if (config("user_debug") > 3) { printf("\t will try alg1 %d loops\n", loops); } tlen = mul_loop(loops, &x); if (config("user_debug") > 3) { printf("\t alg1 time = %.3f secs\n", tlen); } tover = mul_loop(loops, &one); if (config("user_debug") > 3) { printf("\t alg1 overhead look %.3f secs\n", tover); } if (tlen <= tover) { quit "mul_ratio: overhead >= loop time"; } alg1_rate = loops / (tlen - tover); if (config("user_debug") > 2) { printf("\tmultiply alg1 rate = %.3f loopsets/sec\n", alg1_rate); } if (alg1_rate <= 0.0) { quit "mul_ratio: alg1 rate was <= 0.0"; } /* * determine the number of loops needed to test 1st alg */ config("mul2", 2),; loops = 1/2; do { loops *= 2; tlen = mul_loop(loops, &x); if (config("user_debug") > 3) { printf("\t alg2 loops %d took %.3f sec\n", loops, tlen); } } while (tlen < 1.0); /* * determine the 2nd algorithm rate */ loops = max(1, ceil(loops * test_time / tlen)); if (loops < 8) { if (config("user_debug") > 1) { printf(" we must expand loop test time to more than %d secs\n", ceil(test_time * (8 / loops))); } loops = 8; } tlen = mul_loop(loops, &x); if (config("user_debug") > 3) { printf("\t alg2 time = %.3f secs\n", tlen); } tover = mul_loop(loops, &one); if (config("user_debug") > 3) { printf("\t alg2 overhead look %.3f secs\n", tover); } if (tlen <= tover) { quit "mul_ratio: overhead >= loop time"; } alg2_rate = loops / (tlen - tover); if (config("user_debug") > 2) { printf("\tmultiply alg2 rate = %.3f loopsets/sec\n", alg2_rate); } if (alg2_rate <= 0.0) { quit "mul_ratio: alg2 rate was <= 0.0"; } /* * restore old config */ config("all", orig_cfg),; /* * return alg1 / alg2 rate ratio */ return (alg1_rate / alg2_rate); } /* * best_mul2 - determine the best config("mul2") parameter * * NOTE: Due to precision problems with CPU measurements, it is not * unusual for the output of this function to vary slightly * from run to run. * * NOTE: This function is designed to take a long time to run. * We recommend setting: * * config("user_debug", 2) * * so that yon can watch the progress of this function. */ define best_mul2() { local ratio; /* previously calculated alg1/alg2 ratio */ local low; /* low loop value tested */ local high; /* low loop value tested */ local expand; /* how fast to expand the length */ /* * setup */ test_time = 15.0; if (config("user_debug") > 0) { printf("will start with loop test time of %d secs\n", test_time); } /* * firewall - must have a >1 ratio for the initial length */ high = 16; if (config("user_debug") > 0) { printf("testing multiply alg1/alg2 ratio for len = %d\n", high); } ratio = mul_ratio(high); if (config("user_debug") > 1) { printf(" multiply alg1/alg2 ratio = %.3f\n", ratio); } if (ratio <= 1.0) { quit "best_mul2: tests imply config(\"mul2\") should be < 4"; } /* * expand lengths until the ratio flips */ do { /* * determine the paramters of the next ratio test * * We will multiplicatively expand our test level until * the ratio drops below 1.0. */ expand = ((ratio >= 3.5) ? 16 : 2^round(ratio)); low = high; high *= expand; if (config("user_debug") > 1) { printf(" expand the next test range by a factor of %d\n", expand); } /* * determine the alg1/alg2 test ratio for this new length */ if (high >= 2^31) { quit "best_mul2: tests imply config(\"mul2\") should be >= 2^31"; } if (config("user_debug") > 0) { printf("testing multiply alg1/alg2 ratio for len = %d\n", high); } ratio = mul_ratio(high); if (config("user_debug") > 1) { printf(" multiply alg1/alg2 ratio = %.3f\n", ratio); } } while (ratio >= 1.0); if (config("user_debug") > 0) { printf("alg1/alg2 ratio now < 1.0, starting binary search between %d and %d\n", low, high); } /* * binary search between low and high, for where ratio is just under 1.0 */ while (low+1 < high) { /* try the mid-point */ if (config("user_debug") > 0) { printf("testing multiply alg1/alg2 ratio for len = %d\n", int((low+high)/2)); } ratio = mul_ratio(int((low+high)/2)); if (config("user_debug") > 1) { printf(" multiply alg1/alg2 ratio = %.3f\n", ratio); } /* bump lower range up if we went over */ if (ratio >= 1.0) { if (config("user_debug") > 2) { printf("\tmove low from %d up to %d\n", low, int((low+high)/2)); } low = int((low+high)/2); /* drop higher range down if we went under */ } else { if (config("user_debug") > 2) { printf("\tmove high from %d down to %d\n", high, int((low+high)/2)); } high = int((low+high)/2); } } /* * return on the suggested config("mul2") value */ if (config("user_debug") > 0) { printf("best value of config(\"mul2\") is %d\n", (ratio >= 1.0) ? high : low); } return ((ratio >= 1.0) ? high : low); } /* * sq_loop - measure the CPU time to perform a set of square loops * * given: * repeat number :Ò;Ò<Ò=Ò>Ò?Ò@ÒAÒBÒCÒDÒEÒFÒGÒHÒIÒJÒKÒLÒMÒNÒOÒPÒof square loops to perform * x array of 5 values, each the same length in BASEB-bit words * * NOTE: When their lengths are 1 BASEB-bit word, then a * dummy loop of simple constants are used. Thus the * length == 1 is an approximation of loop overhead. * returns: * approximate runtime to perform a square loop * * NOTE: This is an internal support function that is normally * not called directly from the command line. Call the * function best_sq2() instead. */ define sq_loop(repeat, x) { local start; /* start of execution */ local end; /* end of execution */ local answer; /* squared value */ local len; /* length of each element */ local baseb_bytes; /* bytes in a BASEB-bit word */ local i; /* firewall */ if (!isint(repeat) || repeat < 0) { quit "sq_loop: 1st arg: repeat must be an integer > 0"; } if (size(*x) != 5) { quit "sq_loop: 2nd arg matrix does not have 5 elements"; } if (matdim(*x) != 1) { quit "sq_loop: 2nd arg matrix is not 1 dimensional"; } if (matmin(*x, 1) != 0) { quit "sq_loop: 2nd arg matrix index range does not start with 0"; } if (matmax(*x, 1) != 4) { quit "sq_loop: 2nd arg matrix index range does not end with 4"; } baseb_bytes = config("baseb") / 8; len = sizeof((*x)[0]) / baseb_bytes; for (i=1; i < 4; ++i) { if ((sizeof((*x)[i]) / baseb_bytes) != len) { quit "sq_loop: 2nd arg matrix elements are not of equal BASEB-bit word length"; } } /* square pairwise, all sets of a given length */ start = usertime(); for (i=0; i < repeat; ++i) { if (len == 1) { /* we use len == 1 to test this tester loop overhead */ answer = 0^2; answer = 0^2; answer = 0^2; answer = 0^2; answer = 0^2; /**/ answer = 0^2; answer = 0^2; answer = 0^2; answer = 0^2; answer = 0^2; /**/ answer = 0^2; answer = 0^2; answer = 0^2; answer = 0^2; answer = 0^2; /**/ answer = 0^2; answer = 0^2; answer = 0^2; answer = 0^2; answer = 0^2; } else { /* one square loop */ answer = (*x)[0]^2; answer = (*x)[1]^2; answer = (*x)[2]^2; answer = (*x)[3]^2; answer = (*x)[4]^2; /**/ answer = (*x)[0]^2; answer = (*x)[1]^2; answer = (*x)[2]^2; answer = (*x)[3]^2; answer = (*x)[4]^2; /**/ answer = (*x)[0]^2; answer = (*x)[1]^2; answer = (*x)[2]^2; answer = (*x)[3]^2; answer = (*x)[4]^2; /**/ answer = (*x)[0]^2; answer = (*x)[1]^2; answer = (*x)[2]^2; answer = (*x)[3]^2; answer = (*x)[4]^2; } } /* * return duration */ end = usertime(); return end-start; } /* * sq_ratio - ratio of rates of 1st and 2nd square algorithms * * given: * len length in BASEB-bit words to square * * return: * ratio of (1st / 2nd) algorithm rates * * When want to determine a rate to a precision of 1 part in 1000. * Most systems today return CPU time to at least 10 msec precision. * So to get rates to that precision, we need to time loops to at * least 1000 times as long as the precision (10 msec * 1000) * which usually requires timing of loops that last 10 seconds or more. * * NOTE: This is an internal support function that is normally * not called directly from the command line. Call the * function best_sq2() instead. */ define sq_ratio(len) { local mat x[5]; /* array of values for sq_loop to square */ local mat one[5]; /* array if single BASEB-bit values */ local baseb; /* calc word size in bits */ local orig_cfg; /* caller configuration */ local loops; /* number of square loops to time */ local tlen; /* time to perform some number of loops */ local tover; /* est of time for loop overhead */ local alg1_rate; /* loop rate of 1st algorithm */ local alg2_rate; /* loop rate of 2nd algorithm */ local i; /* * firewall */ if (!isint(len) || len < 2) { quit "sq_ratio: 1st arg: len is not an integer > 1"; } /* * remember the caller's config state */ orig_cfg = config("all"); config("mul2", 0),; config("sq2", 0),; config("pow2", 0),; config("redc2", 0),; config("tilde", 0),; /* * initialize x, the values we will square * * We want these tests to be repeatable as possible, so we will seed * the PRNG in a deterministic way. */ baseb = config("baseb"); srand(sha1(sha1(baseb, config("version")))); for (i=0; i < 5; ++i) { /* force the values to be a full len words long */ x[i] = ((1<<(((len-1) * baseb) + baseb-1)) | randbit(((len-1) * baseb) + baseb-2)); /* single BASEB-bit values */ one[i] = 1; } /* * determine the number of loops needed to test 1st alg */ config("sq2", 2^31-1),; loops = 1/2; do { loops *= 2; tlen = sq_loop(loops, &x); if (config("user_debug") > 3) { printf("\t alg1 loops %d took %.3f sec\n", loops, tlen); } } while (tlen < 1.0); /* * determine the 1st algorithm rate */ loops = max(1, ceil(loops * test_time / tlen)); if (loops < 8) { if (config("user_debug") > 1) { printf(" we must expand loop test time to more than %d secs\n", ceil(test_time * (8 / loops))); } loops = 8; } tlen = sq_loop(loops, &x); if (config("user_debug") > 3) { printf("\t alg1 time = %.3f secs\n", tlen); } tover = sq_loop(loops, &one); if (config("user_debug") > 3) { printf("\t alg1 overhead look %.3f secs\n", tover); } if (tlen <= tover) { quit "sq_ratio: overhead >= loop time"; } alg1_rate = loops / (tlen - tover); if (config("user_debug") > 2) { printf("\tsquare alg1 rate = %.3f loopsets/sec\n", alg1_rate); } if (alg1_rate <= 0.0) { quit "sq_ratio: alg1 rate was <= 0.0"; } /* * determine the number of loops needed to test 1st alg */ config("sq2", 2),; loops = 1/2; do { loops *= 2; tlen = sq_loop(loops, &x); if (config("user_debug") > 3) { printf("\t alg2 loops %d took %.3f sec\n", loops, tlen); }