Release calc version 2.10.2t30

lib/Makefile

@@ -0,0 +1,116 @@
+#
+# lib - makefile for calc library scripts
+#
+# Copyright (c) 1996 David I. Bell and Landon Curt Noll
+# Permission is granted to use, distribute, or modify this source,
+# provided that this copyright notice remains intact.
+#
+# Arbitrary precision calculator.
+#
+# calculator by David I. Bell
+# makefile by Landon Curt Noll
+
+# required vars
+#
+SHELL = /bin/sh
+MAKE_FILE = Makefile
+
+# Normally, the upper level makefile will set these values.  We provide
+# a default here just in case you want to build from this directory.
+#
+# where to install things
+TOPDIR= /usr/local/lib
+#TOPDIR= /usr/lib
+#TOPDIR= /usr/libdata
+
+LIBDIR= ${TOPDIR}/calc
+
+# Makefile debug
+#
+# Q=@	do not echo internal makefile actions (quiet mode)
+# Q=	echo internal makefile actions (debug / verbose mode)
+#
+#Q=
+Q=@
+
+# The calc files to install
+#
+CALC_FILES= README bigprime.cal deg.cal ellip.cal lucas.cal lucas_chk.cal \
+	lucas_tbl.cal mersenne.cal mod.cal pell.cal pi.cal \
+	pollard.cal poly.cal psqrt.cal quat.cal regress.cal solve.cal \
+	sumsq.cal surd.cal unitfrac.cal varargs.cal chrem.cal mfactor.cal \
+	bindings altbind randmprime.cal test1700.cal randrun.cal \
+	randbitrun.cal cryrand.cal bernoulli.cal test2300.cal test2600.cal \
+	test2700.cal test3100.cal test3300.cal test3400.cal prompt.cal \
+	test3500.cal seedrandom.cal test4000.cal test4100.cal test4600.cal
+
+# These files are found (but not built) in the distribution
+#
+DISTLIST= ${CALC_FILES} ${MAKE_FILE}
+
+SHELL= /bin/sh
+
+all: ${CALC_FILES} ${MAKE_FILE} .all
+
+# used by the upper level Makefile to determine of we have done all
+#
+.all:
+	rm -f .all
+	touch .all
+
+##
+#
+# File list generation.  You can ignore this section.
+#
+#
+# We will form the names of source files as if they were in a
+# sub-directory called calc/lib.
+#
+# NOTE: Due to bogus shells found on one common system we must have
+#       an non-emoty else clause for every if condition.  *sigh*
+#
+##
+
+distlist: ${DISTLIST}
+	${Q}for i in ${DISTLIST}; do \
+		echo calc/lib/$$i; \
+	done
+
+# The bsdi distribution has generated files as well as distributed files.
+#
+bsdilist: ${DISTLIST}
+	${Q}for i in ${DISTLIST}; do \
+		echo calc/lib/$$i; \
+	done
+
+clean:
+
+clobber:
+	rm -f .all
+
+install: all
+	-${Q}if [ ! -d ${TOPDIR} ]; then \
+		echo mkdir ${TOPDIR}; \
+		mkdir ${TOPDIR}; \
+	else \
+	    true; \
+	fi
+	-${Q}if [ ! -d ${LIBDIR} ]; then \
+		echo mkdir ${LIBDIR}; \
+		mkdir ${LIBDIR}; \
+	else \
+	    true; \
+	fi
+	${Q}for i in ${CALC_FILES}; do \
+		echo rm -f ${LIBDIR}/$$i; \
+		rm -f ${LIBDIR}/$$i; \
+		echo cp $$i ${LIBDIR}; \
+		cp $$i ${LIBDIR}; \
+		echo chmod 0444 ${LIBDIR}/$$i; \
+		chmod 0444 ${LIBDIR}/$$i; \
+	done
+	${Q}echo remove files that are obsolete
+	-rm -f nextprime.cal nextprim.cal
+	-rm -f test1000.cal test2000.cal ${LIBDIR}/test2000.cal
+	-rm -f ${LIBDIR}/nextprime.cal ${LIBDIR}/nextprim.cal
+	-rm -f ${LIBDIR}/test1000.cal

lib/README

@@ -0,0 +1,489 @@
+# Copyright (c) 1996 David I. Bell and Landon Curt Noll
+# Permission is granted to use, distribute, or modify this source,
+# provided that this copyright notice remains intact.
+
+The following calc library files are provided because they serve as
+examples of how use the calc language, and because the authors thought
+them to be useful!
+
+If you write something that you think is useful, please send it to:
+
+    dbell@auug.org.au
+    chongo@toad.com                 {uunet,pyramid,sun}!hoptoad!chongo
+
+By convention, a lib file only defines and/or initializes functions,
+objects and variables.  (The regression test is an exception.)  Also by
+convention, the a usage message regarding each important object and
+function is printed at the time of the read.
+
+If a lib file needs to load another lib file, it should use the -once
+version of read:
+
+    /* pull in needed library files */
+    read -once "surd"
+    read -once "lucas"
+
+This will cause the needed library files to be read once.  If these
+files have already been read, the read -once will act as a noop.
+
+By convention, the global variable  lib_debug  is used to control
+the verbosity of debug information printed by lib files.  By default,
+the lib_debug has a value of 0.  If lib_debug < 0, then no debug
+messages are printed.  If lib_debug >= 0, then only usage message
+regarding each important object are printed at the time of the read.
+If lib_debug == 0, then only such usage messages are printed; no
+other debug information is printed.
+
+To conform to the above convention, your lib files should end with
+lines of the form:
+
+	global lib_debug;
+	if (lib_debug >= 0) {
+	    print "funcA(side_a, side_b, side_c) defined";
+	    print "funcB(size, mass) defined";
+	}
+
+
+=-=
+
+
+bernoulli.cal
+
+    B(n)
+
+    Calculate the nth Bernoulli number.
+
+
+bigprime.cal
+
+    bigprime(a, m, p)
+
+    A prime test, base a, on p*2^x+1 for even x>m.
+
+
+chrem.cal
+
+    chrem(r1,m1 [,r2,m2, ...])
+    chrem(rlist, mlist)
+
+    Chinese remainder theorem/problem solver.
+
+
+cryrand.cal
+
+    obj cryobj
+    cryrand(len)
+    scryrand([seed, [len1, len2]])
+    scryrand(seed, ip, iq, ir)
+    random([a, [b]])
+    srandom(seed)
+    randstate([cryobj | 0])
+
+    cryptographically strong pseudo-romandom number generator
+
+
+deg.cal
+
+    dms(deg, min, sec)
+    dms_add(a, b)
+    dms_neg(a)
+    dms_sub(a, b)
+    dms_mul(a, b)
+    dms_print(a)
+
+    Calculate in degrees, minutes, and seconds.
+
+
+ellip.cal
+
+    factor(iN, ia, B, force)
+
+    Attempt to factor using the elliptic functions: y^2 = x^3 + a*x + b.
+
+
+lucas.cal
+
+    lucas(h, n)
+
+    Perform a primality test of h*2^n-1, with 1<=h<2*n.
+
+
+lucas_chk.cal
+
+    lucas_chk(high_n)
+
+    Test all primes of the form h*2^n-1, with 1<=h<200 and n <= high_n.
+    Requires lucas.cal to be loaded.  The highest useful high_n is 1000.
+
+    Used by regress.cal during the 2100 test set.
+
+
+lucas_tbl.cal
+
+    Lucasian criteria for primality tables.
+
+
+mersenne.cal
+
+    mersenne(p)
+
+    Perform a primality test of 2^p-1, for prime p>1.
+
+
+mfactor.cal
+
+    mfactor(n [, start_k [, rept_loop])
+
+    Return the lowest factor of 2^n-1, for n > 0.  Starts looking for factors
+    at 2*start_k*n+1.  By default, start_k == 1.
+
+    Be default, mfactor() does not report the search progress.  When
+    rept_loop > 0, then a report is given every 4*rept_loop loops.
+
+
+mod.cal
+
+    mod(a)
+    mod_print(a)
+    mod_one()
+    mod_cmp(a, b)
+    mod_rel(a, b)
+    mod_add(a, b)
+    mod_sub(a, b)
+    mod_neg(a)
+    mod_mul(a, b)
+    mod_square(a)
+    mod_inc(a)
+    mod_dec(a)
+    mod_inv(a)
+    mod_div(a, b)
+    mod_pow(a, b)
+
+    Routines to handle numbers modulo a specified number.
+
+
+pell.cal
+
+    pellx(D)
+    pell(D)
+
+    Solve Pell's equation; Returns the solution X to: X^2 - D * Y^2 = 1.
+    Type the solution to pells equation for a particular D.
+
+
+pi.cal
+
+    qpi(epsilon)
+
+    Calculate pi within the specified epsilon using the quartic convergence
+    iteration.
+
+
+pollard.cal
+
+    factor(N, N, ai, af)
+
+    Factor using Pollard's p-1 method.
+
+
+poly.cal
+
+    Calculate with polynomials of one variable.  There are many functions.
+    Read the documentation in the library file.
+
+
+prompt.cal
+
+    adder()
+    showvalues(str)
+
+    Demonstration of some uses of prompt() and eval().
+
+
+psqrt.cal
+
+    psqrt(u, p)
+
+    Calculate square roots modulo a prime
+
+
+quat.cal
+
+    quat(a, b, c, d)
+    quat_print(a)
+    quat_norm(a)
+    quat_abs(a, e)
+    quat_conj(a)
+    quat_add(a, b)
+    quat_sub(a, b)
+    quat_inc(a)
+    quat_dec(a)
+    quat_neg(a)
+    quat_mul(a, b)
+    quat_div(a, b)
+    quat_inv(a)
+    quat_scale(a, b)
+    quat_shift(a, b)
+
+    Calculate using quaternions of the form: a + bi + cj + dk.  In these
+    functions, quaternians are manipulated in the form: s + v, where
+    s is a scalar and v is a vector of size 3.
+
+
+randbitrun.cal
+
+    randbitrun([run_cnt])
+
+    Using randbit(1) to generate a sequence of random bits, determine if
+    the number and kength of identical bits runs match what is expected.
+    By default, run_cnt is to test the next 65536 random values.
+
+
+randmprime.cal
+
+    randmprime(bits, seed [,dbg])
+
+    Find a prime of the form h*2^n-1 >= 2^bits for some given x.  The initial
+    search points for 'h' and 'n' are selected by a cryptographic pseudo-random
+    number generator.  The optional argument, dbg, if set to 1, 2 or 3
+    turn on various debugging print statements.
+
+
+randrun.cal
+
+    randrun([run_cnt])
+
+    Perform the "G. Run test" (pp. 65-68) as found in Knuth's "Art of
+    Computer Programming - 2nd edition", Volume 2, Section 3.3.2 on
+    the builtin rand() function.  This function will generate run_cnt
+    64 bit values.  By default, run_cnt is to test the next 65536
+    random values.
+
+
+regress.cal
+
+    Test the correct execution of the calculator by reading this library file.
+    Errors are reported with '****' mssages, or worse.  :-)
+
+
+seedrandom.cal
+
+    seedrandom(seed1, seed2, bitsize [,trials])
+
+    Given:
+	seed1 - a large random value (at least 10^20 and perhaps < 10^93)
+	seed2 - a large random value (at least 10^20 and perhaps < 10^93)
+ 	size - min Blum modulus as a power of 2 (at least 100, perhaps > 1024)
+	trials - number of ptest() trials (default 25) (optional arg)
+
+    Returns:
+	the previous random state
+
+    Seed the cryptographically strong Blum generator.  This functions allows
+    one to use the raw srandom() without the burden of finding appropriate
+    Blum primes for the modulus.
+
+
+solve.cal
+
+    solve(low, high, epsilon)
+
+    Solve the equation f(x) = 0 to within the desired error value for x.
+    The function 'f' must be defined outside of this routine, and the low
+    and high values are guesses which must produce values with opposite signs.
+
+
+sumsq.cal
+
+    ss(p)
+
+    Determine the unique two positive integers whose squares sum to the
+    specified prime.  This is always possible for all primes of the form
+    4N+1, and always impossible for primes of the form 4N-1.
+
+
+surd.cal
+
+    surd(a, b)
+    surd_print(a)
+    surd_conj(a)
+    surd_norm(a)
+    surd_value(a, xepsilon)
+    surd_add(a, b)
+    surd_sub(a, b)
+    surd_inc(a)
+    surd_dec(a)
+    surd_neg(a)
+    surd_mul(a, b)
+    surd_square(a)
+    surd_scale(a, b)
+    surd_shift(a, b)
+    surd_div(a, b)
+    surd_inv(a)
+    surd_sgn(a)
+    surd_cmp(a, b)
+    surd_rel(a, b)
+
+    Calculate using quadratic surds of the form: a + b * sqrt(D).
+
+
+test1700.cal
+
+    value
+
+    This script is used by regress.cal to test the read and use keywords.
+
+
+test2600.cal
+
+    global defaultverbose
+    global err
+    testismult(str, n, verbose)
+    testsqrt(str, n, eps, verbose)
+    testexp(str, n, eps, verbose)
+    testln(str, n, eps, verbose)
+    testpower(str, n, b, eps, verbose)
+    testgcd(str, n, verbose)
+    cpow(x, n, eps)
+    cexp(x, eps)
+    cln(x, eps)
+    mkreal()
+    mkcomplex()
+    mkbigreal()
+    mksmallreal()
+    testappr(str, n, verbose)
+    checkappr(x, y, z, verbose)
+    checkresult(x, y, z, a)
+    test2600(verbose, tnum)
+
+    This script is used by regress.cal to test some of builtin functions
+    in terms of accuracy and roundoff.
+
+
+test2700.cal
+
+    global defaultverbose
+    mknonnegreal()
+    mkposreal()
+    mkreal_2700()
+    mknonzeroreal()
+    mkposfrac()
+    mkfrac()
+    mksquarereal()
+    mknonsquarereal()
+    mkcomplex_2700()
+    testcsqrt(str, n, verbose)
+    checksqrt(x, y, z, v)
+    checkavrem(A, B, X, eps)
+    checkrounding(s, n, t, u, z)
+    iscomsq(x)
+    test2700(verbose, tnum)
+
+    This script is used by regress.cal to test sqrt() for real and complex
+    values.
+
+
+test3100.cal
+
+    obj res
+    global md
+    res_test(a)
+    res_sub(a, b)
+    res_mul(a, b)
+    res_neg(a)
+    res_inv(a)
+    res(x)
+
+    This script is used by regress.cal to test determinants of a matrix
+
+
+test3300.cal
+
+    global defaultverbose
+    global err
+    testi(str, n, N, verbose)
+    testr(str, n, N, verbose)
+    test3300(verbose, tnum)
+
+    This script is used by regress.cal to provide for more determinant tests.
+
+
+test3400.cal
+
+    global defaultverbose
+    global err
+    test1(str, n, eps, verbose)
+    test2(str, n, eps, verbose)
+    test3(str, n, eps, verbose)
+    test4(str, n, eps, verbose)
+    test5(str, n, eps, verbose)
+    test6(str, n, eps, verbose)
+    test3400(verbose, tnum)
+
+    This script is used by regress.cal to test trig functions.
+    containing objects.
+
+test4000.cal
+
+    global defaultverbose
+    global err
+    global BASEB
+    global BASE
+    global COUNT
+    global SKIP
+    global RESIDUE
+    global MODULUS
+    global K1
+    global H1
+    global K2
+    global H2
+    global K3
+    global H3
+    plen(N) defined
+    rlen(N) defined
+    clen(N) defined
+    ptimes(str, N, n, count, skip, verbose) defined
+    ctimes(str, N, n, count, skip, verbose) defined
+    crtimes(str, a, b, n, count, skip, verbose) defined
+    ntimes(str, N, n, count, skip, residue, mod, verbose) defined
+    testnextcand(str, N, n, cnt, skip, res, mod, verbose) defined
+    testnext1(x, y, count, skip, residue, modulus) defined
+    testprevcand(str, N, n, cnt, skip, res, mod, verbose) defined
+    testprev1(x, y, count, skip, residue, modulus) defined
+    test4000(verbose, tnum) defined
+
+    This script is used by regress.cal to test ptest, nextcand and
+    prevcand buildins.
+
+test4100.cal
+
+    global defaultverbose
+    global err
+    global K1
+    global K2
+    global BASEB
+    global BASE
+    rlen_4100(N) defined
+    olen(N) defined
+    test1(x, y, m, k, z1, z2) defined
+    testall(str, n, N, M, verbose) defined
+    times(str, N, n, verbose) defined
+    powtimes(str, N1, N2, n, verbose) defined
+    inittimes(str, N, n, verbose) defined
+    test4100(verbose, tnum) defined
+
+    This script is used by regress.cal to test REDC operations.
+
+unitfrac.cal
+
+    unitfrac(x)
+
+    Represent a fraction as sum of distinct unit fractions.
+
+
+varargs.cal
+
+    sc(a, b, ...)
+
+    Example program to use 'varargs'.  Program to sum the cubes of all
+    the specified numbers.

lib/altbind

@@ -0,0 +1,45 @@
+# Alternate key bindings for calc line editing functions
+
+map	base-map
+default	insert-char
+^@	set-mark
+^A	start-of-line
+^B	backward-char
+^D	quit
+^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
+^?	delete-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

lib/bernoulli.cal

@@ -0,0 +1,67 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Calculate the Nth Bernoulli number B(n).
+ * This uses the following symbolic formula to calculate B(n):
+ *
+ *	(b+1)^(n+1) - b^(n+1) = 0
+ *
+ * where b is a dummy value, and each power b^i gets replaced by B(i).
+ * For example, for n = 3:
+ *	(b+1)^4 - b^4 = 0
+ *	b^4 + 4*b^3 + 6*b^2 + 4*b + 1 - b^4 = 0
+ *	4*b^3 + 6*b^2 + 4*b + 1 = 0
+ *	4*B(3) + 6*B(2) + 4*B(1) + 1 = 0
+ *	B(3) = -(6*B(2) + 4*B(1) + 1) / 4
+ *
+ * The combinatorial factors in the expansion of the above formula are
+ * calculated interatively, and we use the fact that B(2i+1) = 0 if i > 0.
+ * Since all previous B(n)'s are needed to calculate a particular B(n), all
+ * values obtained are saved in an array for ease in repeated calculations.
+ */
+static Bnmax;
+static mat Bn[1001];
+
+
+define B(n)
+{
+	local	nn, np1, i, sum, mulval, divval, combval;
+
+	if (!isint(n) || (n < 0))
+		quit "Non-negative integer required for Bernoulli";
+
+	if (n == 0)
+		return 1;
+	if (n == 1)
+		return -1/2;
+	if (isodd(n))
+		return 0;
+	if (n > 1000)
+		quit "Very large Bernoulli";
+
+	if (n <= Bnmax)
+		return Bn[n];
+
+	for (nn = Bnmax + 2; nn <= n; nn+=2) {
+		np1 = nn + 1;
+		mulval = np1;
+		divval = 1;
+		combval = 1;
+		sum = 1 - np1 / 2;
+		for (i = 2; i < np1; i+=2) {
+			combval = combval * mulval-- / divval++;
+			combval = combval * mulval-- / divval++;
+			sum += combval * Bn[i];
+		}
+		Bn[nn] = -sum / np1;
+	}
+	Bnmax = n;
+	return Bn[n];
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "B(n) defined";
+}

lib/bigprime.cal

@@ -0,0 +1,32 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * A prime test, base a, on p*2^x+1 for even x>m.
+ */
+
+define bigprime(a, m, p)
+{
+	local n1, n;
+
+	n1 = 2^m * p;
+	for (;;) {
+		m++;
+		n1 += n1;
+		n = n1 + 1;
+		if (isodd(m))
+			continue;
+		print m;
+		if (pmod(a, n1 / 2, n) != n1)
+			continue;
+		if (pmod(a, n1 / p, n) == 1)
+			continue;
+		print "	" : n;
+	}
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "bigprime(a, m, p) defined";
+}

lib/bindings

@@ -0,0 +1,45 @@
+# Default key bindings for calc line editing functions
+
+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

lib/chrem.cal

@@ -0,0 +1,181 @@
+/*
+ * chrem - Chinese remainder theorem/problem solver
+ *
+ * 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.
+ *
+ * Written by: Ernest W Bowen <ernie@neumann.une.edu.au>
+ * Interface by: Landon Curt Noll <chongo@toad.com>
+ */
+
+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;
+    }
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "chrem(r1,m1 [,r2,m2 ...]) defined";
+    print "chrem(rlist [,mlist]) defined";
+}

lib/deg.cal

@@ -0,0 +1,124 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Calculate in degrees, minutes, and seconds.
+ */
+
+obj dms {deg, min, sec};
+
+define dms(deg, min, sec)
+{
+	local ans;
+
+	if (isnull(sec))
+		sec = 0;
+	if (isnull(min))
+		min = 0;
+	obj dms ans;
+	ans.deg = deg;
+	ans.min = min;
+	ans.sec = sec;
+	fixdms(&ans);
+	return ans;
+}
+
+
+define dms_add(a, b)
+{
+	local obj dms	ans;
+
+	ans.deg = 0;
+	ans.min = 0;
+	ans.sec = 0;
+	if (istype(a, ans)) {
+		ans.deg += a.deg;
+		ans.min += a.min;
+		ans.sec += a.sec;
+	} else
+		ans.deg += a;
+	if (istype(b, ans)) {
+		ans.deg += b.deg;
+		ans.min += b.min;
+		ans.sec += b.sec;
+	} else
+		ans.deg += b;
+	fixdms(&ans);
+	return ans;
+}
+
+
+define dms_neg(a)
+{
+	local obj dms	ans;
+
+	ans.deg = -ans.deg;
+	ans.min = -ans.min;
+	ans.sec = -ans.sec;
+	return ans;
+}
+
+
+define dms_sub(a, b)
+{
+	return a - b;
+}
+
+
+define dms_mul(a, b)
+{
+	local obj dms	ans;
+
+	if (istype(a, ans) && istype(b, ans))
+		quit "Cannot multiply degrees together";
+	if (istype(a, ans)) {
+		ans.deg = a.deg * b;
+		ans.min = a.min * b;
+		ans.sec = a.sec * b;
+	} else {
+		ans.deg = b.deg * a;
+		ans.min = b.min * a;
+		ans.sec = b.sec * a;
+	}
+	fixdms(&ans);
+	return ans;
+}
+
+
+define dms_print(a)
+{
+	print a.deg : 'd' : a.min : 'm' : a.sec : 's' :;
+}
+
+
+define dms_abs(a)
+{
+	return a.deg + a.min / 60 + a.sec / 3600;
+}
+
+
+define fixdms(a)
+{
+	a.min += frac(a.deg) * 60;
+	a.deg = int(a.deg);
+	a.sec += frac(a.min) * 60;
+	a.min = int(a.min);
+	a.min += a.sec // 60;
+	a.sec %= 60;
+	a.deg += a.min // 60;
+	a.min %= 60;
+	a.deg %= 360;
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "obj dms {deg, min, sec} defined";
+    print "dms(deg, min, sec) defined";
+    print "dms_add(a, b) defined";
+    print "dms_neg(a) defined";
+    print "dms_sub(a, b) defined";
+    print "dms_mul(a, b) defined";
+    print "dms_print(a) defined";
+    print "dms_abs(a) defined";
+}

lib/ellip.cal

@@ -0,0 +1,172 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Attempt to factor numbers using elliptic functions.
+ * 	y^2 = x^3 + a*x + b   (mod N).
+ *
+ * Many points (x,y) (mod N) are found that solve the above equation,
+ * starting from a trivial solution and 'multiplying' that point together
+ * to generate high powers of the point, looking for such a point whose
+ * order contains a common factor with N.  The order of the group of points
+ * varies almost randomly within a certain interval for each choice of a
+ * and b, and thus each choice provides an independent opportunity to
+ * factor N.  To generate a trivial solution, a is chosen and then b is
+ * selected so that (1,1) is a solution.  The multiplication is done using
+ * the basic fact that the equation is a cubic, and so if a line hits the
+ * curve in two rational points, then the third intersection point must
+ * also be rational.  Thus by drawing lines between known rational points
+ * the number of rational solutions can be made very large.  When modular
+ * arithmetic is used, solving for the third point requires the taking of a
+ * modular inverse (instead of division), and if this fails, then the GCD
+ * of the failing value and N provides a factor of N.  This description is
+ * only an approximation, read "A Course in Number Theory and Cryptography"
+ * by Neal Koblitz for a good explanation.
+ *
+ * factor(iN, ia, B, force)
+ *	iN is the number to be factored.
+ *	ia is the initial value of a in the equation, and each successive
+ *	value of a is an independent attempt at factoring (default 1).
+ *	B is the limit of the primes that make up the high power that the
+ *	point is raised to for each factoring attempt (default 100).
+ *	force is a flag to attempt to factor numbers even if they are
+ *	thought to already be prime (default FALSE).
+ *
+ * Making B larger makes the power the point being raised to contain more
+ * prime factors, thus increasing the chance that the order of the point
+ * will be made up of those factors.  The higher B is then, the greater
+ * the chance that any individual attempt will find a factor.  However,
+ * a higher B also slows down the number of independent functions being
+ * examined.  The order of the point for any particular function might
+ * contain a large prime and so won't succeed even for a really large B,
+ * whereas the next function might have an order which is quickly found.
+ * So you want to trade off the depth of a particular search with the
+ * number of searches made.  For example, for factoring 30 digits, I make
+ * B be about 1000 (probably still too small).
+ *
+ * If you have lots of machines available, then you can run parallel
+ * factoring attempts for the same number by giving different starting
+ * values of ia for each machine (e.g. 1000, 2000, 3000).
+ *
+ * The output as the function is running is (occasionally) the value of a
+ * when a new function is started, the prime that is being included in the
+ * high power being calculated, and the current point which is the result
+ * of the powers so far.
+ *
+ * If a factor is found, it is returned and is also saved in the global
+ * variable f.  The number being factored is also saved in the global
+ * variable N.
+ */
+
+obj point {x, y};
+global	N;		/* number to factor */
+global	a;		/* first coefficient */
+global	b;		/* second coefficient */
+global	f;		/* found factor */
+
+
+define factor(iN, ia, B, force)
+{
+	local	C, x, p;
+
+	if (!force && ptest(iN, 50))
+		return 1;
+	if (isnull(B))
+		B = 100;
+	if (isnull(ia))
+		ia = 1;
+	obj point x;
+	a = ia;
+	b = -ia;
+	N = iN;
+	C = isqrt(N);
+	C = 2 * C + 2 * isqrt(C) + 1;
+	f = 0;
+	while (f == 0) {
+		print "A =", a;
+		x.x = 1;
+		x.y = 1;
+		print 2, x;
+		x = x ^ (2 ^ (highbit(C) + 1));
+		for (p = 3; ((p < B) && (f == 0)); p += 2) {
+			if (!ptest(p, 1))
+				continue;
+			print p, x;
+			x = x ^ (p ^ ((highbit(C) // highbit(p)) + 1));
+		}
+		a++;
+		b--;
+	}
+	return f;
+}
+
+
+define point_print(p)
+{
+	print "(" : p.x : "," : p.y : ")" :;
+}
+
+
+define point_mul(p1, p2)
+{
+	local	r, m;
+
+	if (p2 == 1)
+		return p1;
+	if (p1 == p2)
+		return point_square(&p1);
+	obj point r;
+	m = (minv(p2.x - p1.x, N) * (p2.y - p1.y)) % N;
+	if (m == 0) {
+		if (f == 0)
+			f = gcd(p2.x - p1.x, N);
+		r.x = 1;
+		r.y = 1;
+		return r;
+	}
+	r.x = (m^2 - p1.x - p2.x) % N;
+	r.y = ((m * (p1.x - r.x)) - p1.y) % N;
+	return r;
+}
+
+
+define point_square(p)
+{
+	local	r, m;
+
+	obj point r;
+	m = ((3 * p.x^2 + a) * minv(p.y << 1, N)) % N;
+	if (m == 0) {
+		if (f == 0)
+			f = gcd(p.y << 1, N);
+		r.x = 1;
+		r.y = 1;
+		return r;
+	}
+	r.x = (m^2 - p.x - p.x) % N;
+	r.y = ((m * (p.x - r.x)) - p.y) % N;
+	return r;
+}
+
+
+define point_pow(p, pow)
+{
+	local bit, r, t;
+
+	r = 1;
+	if (isodd(pow))
+		r = p;
+	t = p;
+	for (bit = 2; ((bit <= pow) && (f == 0)); bit <<= 1) {
+		t = point_square(&t);
+		if (bit & pow)
+			r = point_mul(&t, &r);
+	}
+	return r;
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "factor(N, I, B, force) defined";
+}

lib/lucas_chk.cal

@@ -0,0 +1,381 @@
+/*
+ * Copyright (c) 1995 Landon Curt Noll
+ *
+ * Permission to use, copy, modify, and distribute this software and
+ * its documentation for any purpose and without fee is hereby granted,
+ * provided that the above copyright, this permission notice and text
+ * this comment, and the disclaimer below appear in all of the following:
+ *
+ *	supporting documentation
+ *	source copies
+ *	source works derived from this source
+ *	binaries derived from this source or from derived source
+ *
+ * LANDON CURT NOLL DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+ * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY 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.
+ *
+ * chongo was here	/\../\		chongo@toad.com
+ */
+/*
+ * primes of the form h*2^n-1 for 1<=h<200 and 1<=n<1000
+ *
+ * For all 0 <= i < prime_cnt, h_p[i]*2^n_p[i]-1 is prime.
+ *
+ * These values were taken from:
+ *
+ *	"Prime numbers and Computer Methods for Factorization", by Hans Riesel,
+ *	Birkhauser, 1985, pp 384-387.
+ *
+ * This routine assumes that the file "lucas.cal" has been loaded.
+ *
+ * NOTE: There are several errors in Riesel's table that have been corrected
+ *	 in this file:
+ *
+ *		193*2^87-1 is prime
+ *		193*2^97-1 is NOT prime
+ *		199*2^211-1 is prime
+ *		199*2^221-1 is NOT prime
+ */
+
+static prime_cnt = 1145;	/* number of primes in the list */
+
+/* h = prime parameters */
+static mat h_p[prime_cnt] = {
+	1, 1, 1, 1, 1, 1, 1, 1, 1, 1,			/* element 0 */
+	1, 1, 1, 1, 3, 3, 3, 3, 3, 3,
+	3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
+	3, 3, 3, 3, 3, 3, 3, 3, 3, 5,
+	5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
+	5, 5, 5, 5, 5, 5, 7, 7, 7, 7,
+	7, 7, 7, 7, 9, 9, 9, 9, 9, 9,
+	9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
+	9, 9, 9, 11, 11, 11, 11, 11, 11, 11,
+	11, 11, 11, 13, 13, 13, 13, 13, 13, 15,
+	15, 15, 15, 15, 15, 15, 15, 15, 15, 15,		/* 100 */
+	15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
+	15, 15, 17, 17, 17, 17, 17, 17, 17, 17,
+	17, 17, 17, 17, 17, 17, 17, 17, 17, 17,
+	17, 17, 19, 19, 19, 19, 19, 19, 19, 19,
+	19, 19, 19, 19, 19, 19, 19, 19, 19, 19,
+	19, 19, 21, 21, 21, 21, 21, 21, 21, 21,
+	21, 21, 21, 21, 21, 21, 21, 21, 23, 23,
+	23, 23, 23, 23, 23, 23, 23, 25, 25, 25,
+	25, 25, 25, 25, 25, 25, 25, 25, 25, 25,
+	25, 25, 25, 27, 27, 27, 27, 27, 27, 27,		/* 200 */
+	27, 27, 27, 27, 27, 27, 27, 27, 27, 27,
+	27, 27, 27, 27, 27, 27, 27, 29, 29, 29,
+	29, 29, 31, 31, 31, 31, 31, 31, 31, 31,
+	31, 31, 31, 31, 31, 31, 31, 31, 31, 31,
+	33, 33, 33, 33, 33, 33, 33, 33, 33, 33,
+	33, 33, 33, 33, 33, 33, 33, 33, 33, 33,
+	33, 33, 33, 33, 35, 35, 35, 35, 35, 35,
+	35, 35, 35, 35, 35, 35, 35, 35, 35, 35,
+	35, 37, 39, 39, 39, 39, 39, 39, 39, 39,
+	39, 41, 41, 41, 41, 41, 41, 41, 41, 41,		/* 300 */
+	41, 41, 41, 41, 43, 43, 43, 43, 43, 45,
+	45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
+	45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
+	45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
+	45, 45, 45, 45, 45, 47, 47, 47, 47, 49,
+	49, 49, 49, 49, 49, 49, 49, 49, 49, 49,
+	49, 49, 49, 49, 49, 49, 51, 51, 51, 51,
+	51, 51, 51, 51, 51, 51, 51, 51, 51, 51,
+	51, 53, 53, 53, 53, 53, 53, 53, 53, 53,
+	53, 55, 55, 55, 55, 55, 55, 55, 55, 55,		/* 400 */
+	55, 55, 55, 55, 55, 55, 55, 55, 55, 55,
+	57, 57, 57, 57, 57, 57, 57, 57, 57, 57,
+	57, 57, 57, 57, 57, 57, 57, 57, 59, 59,
+	59, 59, 59, 59, 61, 61, 61, 61, 61, 61,
+	61, 61, 61, 61, 61, 61, 61, 61, 61, 61,
+	61, 63, 63, 63, 63, 63, 63, 63, 63, 63,
+	63, 63, 63, 63, 63, 63, 63, 63, 63, 63,
+	63, 63, 63, 63, 65, 65, 65, 65, 65, 65,
+	65, 65, 65, 65, 65, 65, 65, 65, 65, 65,
+	65, 65, 67, 67, 67, 67, 67, 67, 67, 67,		/* 500 */
+	69, 69, 69, 69, 69, 69, 69, 69, 69, 69,
+	69, 69, 69, 69, 69, 69, 69, 69, 69, 69,
+	69, 69, 69, 69, 69, 69, 69, 69, 69, 69,
+	69, 69, 71, 71, 71, 73, 73, 73, 73, 73,
+	73, 75, 75, 75, 75, 75, 75, 75, 75, 75,
+	75, 75, 75, 75, 75, 75, 75, 75, 75, 75,
+	75, 75, 75, 75, 75, 75, 75, 77, 77, 77,
+	77, 77, 77, 77, 77, 77, 77, 77, 77, 79,
+	79, 79, 79, 79, 79, 79, 79, 79, 79, 79,
+	81, 81, 81, 81, 81, 81, 81, 81, 81, 81,		/* 600 */
+	81, 81, 81, 83, 83, 83, 83, 83, 83, 83,
+	83, 83, 83, 83, 83, 83, 83, 83, 83, 83,
+	83, 83, 83, 83, 83, 85, 85, 85, 85, 85,
+	85, 85, 85, 85, 87, 87, 87, 87, 87, 87,
+	87, 87, 87, 87, 87, 87, 87, 87, 87, 87,
+	87, 87, 87, 87, 87, 87, 89, 89, 89, 89,
+	89, 89, 89, 89, 89, 91, 91, 91, 91, 91,
+	91, 91, 91, 91, 91, 91, 91, 91, 91, 91,
+	91, 91, 91, 91, 91, 91, 91, 93, 93, 93,
+	93, 93, 93, 93, 93, 93, 93, 93, 93, 93,		/* 700 */
+	93, 93, 93, 93, 93, 95, 95, 95, 95, 95,
+	95, 95, 95, 95, 95, 97, 97, 97, 97, 97,
+	99, 99, 99, 99, 99, 99, 99, 99, 99, 99,
+	99, 99, 99, 99, 99, 99, 101, 101, 101, 101,
+	103, 103, 103, 103, 103, 103, 103, 103, 103, 103,
+	103, 103, 103, 105, 105, 105, 105, 105, 105, 105,
+	105, 105, 105, 105, 105, 105, 105, 105, 105, 105,
+	105, 105, 107, 107, 107, 107, 107, 107, 107, 107,
+	107, 107, 107, 107, 107, 107, 109, 109, 109, 109,
+	109, 113, 113, 113, 113, 113, 113, 113, 113, 113,	/* 800 */
+	113, 115, 115, 115, 115, 115, 115, 115, 115, 115,
+	115, 115, 115, 115, 115, 115, 115, 119, 119, 119,
+	119, 119, 119, 119, 119, 121, 121, 121, 121, 121,
+	121, 121, 121, 121, 121, 121, 121, 125, 125, 125,
+	125, 125, 125, 127, 127, 131, 131, 131, 131, 131,
+	131, 131, 131, 131, 131, 133, 133, 133, 133, 133,
+	133, 133, 133, 133, 133, 133, 133, 133, 137, 137,
+	137, 137, 139, 139, 139, 139, 139, 139, 139, 139,
+	139, 139, 139, 139, 139, 139, 139, 139, 139, 139,
+	139, 139, 139, 139, 139, 139, 139, 139, 139, 143,	/* 900 */
+	143, 143, 143, 143, 143, 143, 143, 143, 143, 143,
+	143, 143, 143, 143, 143, 143, 143, 143, 143, 143,
+	143, 143, 143, 145, 145, 145, 145, 145, 145, 145,
+	145, 145, 145, 145, 149, 149, 149, 149, 149, 149,
+	149, 151, 151, 151, 155, 155, 155, 155, 155, 155,
+	155, 155, 155, 155, 155, 155, 157, 157, 157, 157,
+	157, 157, 157, 157, 157, 161, 161, 161, 161, 161,
+	161, 161, 161, 161, 161, 161, 161, 161, 161, 161,
+	163, 163, 163, 163, 167, 167, 167, 167, 167, 167,
+	167, 167, 167, 167, 167, 167, 169, 169, 169, 169,	/* 1000 */
+	169, 169, 169, 169, 169, 169, 169, 169, 173, 173,
+	173, 173, 173, 173, 173, 173, 173, 173, 173, 173,
+	173, 173, 173, 173, 175, 175, 175, 175, 175, 175,
+	175, 175, 175, 175, 175, 175, 175, 175, 175, 175,
+	179, 179, 179, 181, 181, 181, 181, 181, 181, 181,
+	181, 181, 181, 181, 181, 181, 181, 181, 181, 181,
+	181, 181, 181, 181, 181, 181, 181, 181, 185, 185,
+	185, 185, 185, 185, 185, 185, 185, 185, 187, 187,
+	187, 187, 187, 191, 193, 193, 193, 193, 193, 193,
+	193, 197, 197, 197, 197, 197, 197, 197, 197, 197,	/* 1100 */
+	197, 197, 197, 197, 197, 197, 197, 197, 197, 199,
+	199, 199, 199, 199, 199, 199, 199, 199, 199, 199,
+	199, 199, 199, 199, 199, 199, 199, 199, 199, 199,
+	199, 199, 199, 199, 199
+};
+
+
+/* n (exponent) prime parameters */
+static mat n_p[prime_cnt] = {
+	2, 3, 5, 7, 13, 17, 19, 31, 61, 89,		/* element 0 */
+	107, 127, 521, 607, 1, 2, 3, 4, 6, 7,
+	11, 18, 34, 38, 43, 55, 64, 76, 94, 103,
+	143, 206, 216, 306, 324, 391, 458, 470, 827, 2,
+	4, 8, 10, 12, 14, 18, 32, 48, 54, 72,
+	148, 184, 248, 270, 274, 420, 1, 5, 9, 17,
+	21, 29, 45, 177, 1, 3, 7, 13, 15, 21,
+	43, 63, 99, 109, 159, 211, 309, 343, 415, 469,
+	781, 871, 939, 2, 26, 50, 54, 126, 134, 246,
+	354, 362, 950, 3, 7, 23, 287, 291, 795, 1,
+	2, 4, 5, 10, 14, 17, 31, 41, 73, 80,		/* 100 */
+	82, 116, 125, 145, 157, 172, 202, 224, 266, 289,
+	293, 463, 2, 4, 6, 16, 20, 36, 54, 60,
+	96, 124, 150, 252, 356, 460, 612, 654, 664, 698,
+	702, 972, 1, 3, 5, 21, 41, 49, 89, 133,
+	141, 165, 189, 293, 305, 395, 651, 665, 771, 801,
+	923, 953, 1, 2, 3, 7, 10, 13, 18, 27,
+	37, 51, 74, 157, 271, 458, 530, 891, 4, 6,
+	12, 46, 72, 244, 264, 544, 888, 3, 9, 11,
+	17, 23, 35, 39, 75, 105, 107, 155, 215, 335,
+	635, 651, 687, 1, 2, 4, 5, 8, 10, 14,		/* 200 */
+	28, 37, 38, 70, 121, 122, 160, 170, 253, 329,
+	362, 454, 485, 500, 574, 892, 962, 4, 16, 76,
+	148, 184, 1, 5, 7, 11, 13, 23, 33, 35,
+	37, 47, 115, 205, 235, 271, 409, 739, 837, 887,
+	2, 3, 6, 8, 10, 22, 35, 42, 43, 46,
+	56, 91, 102, 106, 142, 190, 208, 266, 330, 360,
+	382, 462, 503, 815, 2, 6, 10, 20, 44, 114,
+	146, 156, 174, 260, 306, 380, 654, 686, 702, 814,
+	906, 1, 3, 24, 105, 153, 188, 605, 795, 813,
+	839, 2, 10, 14, 18, 50, 114, 122, 294, 362,	/* 300 */
+	554, 582, 638, 758, 7, 31, 67, 251, 767, 1,
+	2, 3, 4, 5, 6, 8, 9, 14, 15, 16,
+	22, 28, 29, 36, 37, 54, 59, 85, 93, 117,
+	119, 161, 189, 193, 256, 308, 322, 327, 411, 466,
+	577, 591, 902, 928, 946, 4, 14, 70, 78, 1,
+	5, 7, 9, 13, 15, 29, 33, 39, 55, 81,
+	95, 205, 279, 581, 807, 813, 1, 9, 10, 19,
+	22, 57, 69, 97, 141, 169, 171, 195, 238, 735,
+	885, 2, 6, 8, 42, 50, 62, 362, 488, 642,
+	846, 1, 3, 5, 7, 15, 33, 41, 57, 69,		/* 400 */
+	75, 77, 131, 133, 153, 247, 305, 351, 409, 471,
+	1, 2, 4, 5, 8, 10, 20, 22, 25, 26,
+	32, 44, 62, 77, 158, 317, 500, 713, 12, 16,
+	72, 160, 256, 916, 3, 5, 9, 13, 17, 19,
+	25, 39, 63, 67, 75, 119, 147, 225, 419, 715,
+	895, 2, 3, 8, 11, 14, 16, 28, 32, 39,
+	66, 68, 91, 98, 116, 126, 164, 191, 298, 323,
+	443, 714, 758, 759, 4, 6, 12, 22, 28, 52,
+	78, 94, 124, 162, 174, 192, 204, 304, 376, 808,
+	930, 972, 5, 9, 21, 45, 65, 77, 273, 677,	/* 500 */
+	1, 4, 5, 7, 9, 11, 13, 17, 19, 23,
+	29, 37, 49, 61, 79, 99, 121, 133, 141, 164,
+	173, 181, 185, 193, 233, 299, 313, 351, 377, 540,
+	569, 909, 2, 14, 410, 7, 11, 19, 71, 79,
+	131, 1, 3, 5, 6, 18, 19, 20, 22, 28,
+	29, 39, 43, 49, 75, 85, 92, 111, 126, 136,
+	159, 162, 237, 349, 381, 767, 969, 2, 4, 14,
+	26, 58, 60, 64, 100, 122, 212, 566, 638, 1,
+	3, 7, 15, 43, 57, 61, 75, 145, 217, 247,
+	3, 5, 11, 17, 21, 27, 81, 101, 107, 327,	/* 600 */
+	383, 387, 941, 2, 4, 8, 10, 14, 18, 22,
+	24, 26, 28, 36, 42, 58, 64, 78, 158, 198,
+	206, 424, 550, 676, 904, 5, 11, 71, 113, 115,
+	355, 473, 563, 883, 1, 2, 8, 9, 10, 12,
+	22, 29, 32, 50, 57, 69, 81, 122, 138, 200,
+	296, 514, 656, 682, 778, 881, 4, 8, 12, 24,
+	48, 52, 64, 84, 96, 1, 3, 9, 13, 15,
+	17, 19, 23, 47, 57, 67, 73, 77, 81, 83,
+	191, 301, 321, 435, 867, 869, 917, 3, 4, 7,
+	10, 15, 18, 19, 24, 27, 39, 60, 84, 111,	/* 700 */
+	171, 192, 222, 639, 954, 2, 6, 26, 32, 66,
+	128, 170, 288, 320, 470, 1, 9, 45, 177, 585,
+	1, 4, 5, 7, 8, 11, 19, 25, 28, 35,
+	65, 79, 212, 271, 361, 461, 10, 18, 54, 70,
+	3, 7, 11, 19, 63, 75, 95, 127, 155, 163,
+	171, 283, 563, 2, 3, 5, 6, 8, 9, 25,
+	32, 65, 113, 119, 155, 177, 299, 335, 426, 462,
+	617, 896, 10, 12, 18, 24, 28, 40, 90, 132,
+	214, 238, 322, 532, 858, 940, 9, 149, 177, 419,
+	617, 8, 14, 74, 80, 274, 334, 590, 608, 614,	/* 800 */
+	650, 1, 3, 11, 13, 19, 21, 31, 49, 59,
+	69, 73, 115, 129, 397, 623, 769, 12, 16, 52,
+	160, 192, 216, 376, 436, 1, 3, 21, 27, 37,
+	43, 91, 117, 141, 163, 373, 421, 2, 4, 44,
+	182, 496, 904, 25, 113, 2, 14, 34, 38, 42,
+	78, 90, 178, 778, 974, 3, 11, 15, 19, 31,
+	59, 75, 103, 163, 235, 375, 615, 767, 2, 18,
+	38, 62, 1, 5, 7, 9, 15, 19, 21, 35,
+	37, 39, 41, 49, 69, 111, 115, 141, 159, 181,
+	201, 217, 487, 567, 677, 765, 811, 841, 917, 2,	/* 900 */
+	4, 6, 8, 12, 18, 26, 32, 34, 36, 42,
+	60, 78, 82, 84, 88, 154, 174, 208, 256, 366,
+	448, 478, 746, 5, 13, 15, 31, 77, 151, 181,
+	245, 445, 447, 883, 4, 16, 48, 60, 240, 256,
+	304, 5, 221, 641, 2, 8, 14, 16, 44, 46,
+	82, 172, 196, 254, 556, 806, 1, 5, 33, 121,
+	125, 305, 445, 473, 513, 2, 6, 18, 22, 34,
+	54, 98, 122, 146, 222, 306, 422, 654, 682, 862,
+	3, 31, 63, 303, 4, 6, 8, 10, 16, 32,
+	38, 42, 52, 456, 576, 668, 1, 5, 11, 17,	/* 1000 */
+	67, 137, 157, 203, 209, 227, 263, 917, 2, 4,
+	6, 16, 32, 50, 76, 80, 96, 104, 162, 212,
+	230, 260, 480, 612, 1, 3, 9, 21, 23, 41,
+	47, 57, 69, 83, 193, 249, 291, 421, 433, 997,
+	8, 68, 108, 3, 5, 7, 9, 11, 17, 23,
+	31, 35, 43, 47, 83, 85, 99, 101, 195, 267,
+	281, 363, 391, 519, 623, 653, 673, 701, 2, 6,
+	10, 18, 26, 40, 46, 78, 230, 542, 1, 17,
+	21, 53, 253, 226, 3, 15, 27, 63, 87, 135,
+	543, 2, 16, 20, 22, 40, 82, 112, 178, 230,	/* 1100 */
+	302, 304, 328, 374, 442, 472, 500, 580, 694, 1,
+	5, 7, 15, 19, 23, 25, 27, 43, 65, 99,
+	125, 141, 165, 201, 211, 331, 369, 389, 445, 461,
+	463, 467, 513, 583, 835
+};
+
+
+/* obtain our required libs */
+read -once "lucas.cal";
+
+
+/*
+ * lucas_chk - check the lucas function on known primes
+ *
+ * This function tests entries in the above h_p, n_p table
+ * when n_p is below a given limit.
+ *
+ * input:
+ *	high_n	skip tests on n_p[i] > high_n
+ *	[quiet]	if given and != 0, then do not print individual test results
+ *
+ * returns:
+ *	1	all is ok
+ *	0	something went wrong
+ */
+define
+lucas_chk(high_n, quiet)
+{
+	local i;	/* index */
+	local result;	/* 0 => non-prime, 1 => prime, -1 => bad test */
+	local error;	/* number of errors and bad tests found */
+
+	/*
+	 * firewall
+	 */
+	if (!isint(high_n)) {
+		ldebug("test_lucas", "high_n is non-int");
+		quit "FATAL: bad args: high_n must be an integer";
+	}
+	if (param(0) == 1) {
+		quiet = 0;
+	}
+
+	/*
+	 * scan thru the above prime table
+	 */
+	error = 0;
+	for (i=0; i < prime_cnt; ++i) {
+
+		/* skip primes where h>=2^n */
+		if (highbit(h_p[i]) >= n_p[i]) {
+			if (lib_debug > 0) {
+				print "h>=2^n skip:", h_p[i]:"*2^":n_p[i]:"-1";
+			}
+			continue;
+		}
+
+		/* test the prime if it is small enough */
+		if (n_p[i] <= high_n) {
+
+			/* test the table value */
+			result = lucas(h_p[i], n_p[i]);
+
+			/* report the test */
+			if (result == 0) {
+				print "ERROR, bad primality test of",\
+				    h_p[i]:"*2^":n_p[i]:"-1";
+				++error;
+			} else if (result == 1) {
+				if (quiet == 0) {
+					print h_p[i]:"*2^":n_p[i]:"-1 is prime";
+				}
+			} else if (result == -1) {
+				print "ERROR, failed to compute v(1) for",\
+				    h_p[i]:"*2^":n_p[i]:"-1";
+				++error;
+			} else {
+				print "ERROR, bogus return value:", result;
+				++error;
+			}
+		}
+	}
+
+	/* return the full status */
+	if (error == 0) {
+		if (quiet == 0) {
+			print "lucas_chk(":high_n:") passed";
+		}
+		return 1;
+	} else if (error == 1) {
+		print "lucas_chk(":high_n:") failed", error, "test";
+		return 0;
+	} else {
+		print "lucas_chk(":high_n:") failed", error, "tests";
+		return 0;
+	}
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "lucas_chk(high_n) defined";
+}

lib/lucas_tbl.cal

@@ -0,0 +1,158 @@
+/*
+ * Copyright (c) 1995 Landon Curt Noll
+ *
+ * Permission to use, copy, modify, and distribute this software and
+ * its documentation for any purpose and without fee is hereby granted,
+ * provided that the above copyright, this permission notice and text
+ * this comment, and the disclaimer below appear in all of the following:
+ *
+ *	supporting documentation
+ *	source copies
+ *	source works derived from this source
+ *	binaries derived from this source or from derived source
+ *
+ * LANDON CURT NOLL DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+ * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY 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.
+ *
+ * chongo was here	/\../\		chongo@toad.com
+ */
+/*
+ * Lucasian criteria for primality
+ *
+ * The following table is taken from:
+ *
+ *	"Lucasian Criteria for the Primality of N=h*2^n-1", by Hans Riesel,
+ *	Mathematics of Computation, Vol 23 #108, p 872.
+ *
+ * The index of the *_val[] arrays correspond to the v(1) values found
+ * in the table.  That is, for v(1) == x:
+ *
+ *	D == d_val[x]
+ *	a == a_val[x]
+ *	b == b_val[x]
+ *	r == r_val[x]		(r == abs(a^2 - b^2*D))
+ *
+ *
+ * Note that when *_val[i] is not a number, the related v(1) value
+ * is not found in Table 1.
+ */
+
+trymax = 100;
+mat d_val[trymax+1];
+mat a_val[trymax+1];
+mat b_val[trymax+1];
+mat r_val[trymax+1];
+/* v1= 0	    INVALID */
+/* v1= 1	    INVALID */
+/* v1= 2	    INVALID */
+d_val[ 3]=   5;  a_val[ 3]= 1;  b_val[ 3]=1;  r_val[ 3]=4;
+d_val[ 4]=   3;  a_val[ 4]= 1;  b_val[ 4]=1;  r_val[ 4]=2;
+d_val[ 5]=  21;  a_val[ 5]= 3;  b_val[ 5]=1;  r_val[ 5]=12;
+d_val[ 6]=   2;  a_val[ 6]= 1;  b_val[ 6]=1;  r_val[ 6]=1;
+/* v1= 7	    INVALID */
+d_val[ 8]=  15;  a_val[ 8]= 3;  b_val[ 8]=1;  r_val[ 8]=6;
+d_val[ 9]=  77;  a_val[ 9]= 7;  b_val[ 9]=1;  r_val[ 9]=28;
+d_val[10]=   6;  a_val[10]= 2;  b_val[10]=1;  r_val[10]=2;
+d_val[11]=  13;  a_val[11]= 3;  b_val[11]=1;  r_val[11]=4;
+d_val[12]=  35;  a_val[12]= 5;  b_val[12]=1;  r_val[12]=10;
+d_val[13]= 165;  a_val[13]=11;  b_val[13]=1;  r_val[13]=44;
+/* v1=14	    INVALID */
+d_val[15]= 221;  a_val[15]=13;  b_val[15]=1;  r_val[15]=52;
+d_val[16]=   7;  a_val[16]= 3;  b_val[16]=1;  r_val[16]=2;
+d_val[17]= 285;  a_val[17]=15;  b_val[17]=1;  r_val[17]=60;
+/* v1=18	    INVALID */
+d_val[19]= 357;  a_val[19]=17;  b_val[19]=1;  r_val[19]=68;
+d_val[20]=  11;  a_val[20]= 3;  b_val[20]=1;  r_val[20]=2;
+d_val[21]= 437;  a_val[21]=19;  b_val[21]=1;  r_val[21]=76;
+d_val[22]=  30;  a_val[22]= 5;  b_val[22]=1;  r_val[22]=5;
+/* v1=23	    INVALID */
+d_val[24]= 143;  a_val[24]=11;  b_val[24]=1;  r_val[24]=22;
+d_val[25]=  69;  a_val[25]= 9;  b_val[25]=1;  r_val[25]=12;
+d_val[26]=  42;  a_val[26]= 6;  b_val[26]=1;  r_val[26]=6;
+d_val[27]=  29;  a_val[27]= 5;  b_val[27]=1;  r_val[27]=4;
+d_val[28]= 195;  a_val[28]=13;  b_val[28]=1;  r_val[28]=26;
+d_val[29]=  93;  a_val[29]= 9;  b_val[29]=1;  r_val[29]=12;
+d_val[30]=  14;  a_val[30]= 4;  b_val[30]=1;  r_val[30]=2;
+d_val[31]= 957;  a_val[31]=29;  b_val[31]=1;  r_val[31]=116;
+d_val[32]= 255;  a_val[32]=15;  b_val[32]=1;  r_val[32]=30;
+d_val[33]=1085;  a_val[33]=31;  b_val[33]=1;  r_val[33]=124;
+/* v1=34	    INVALID */
+d_val[35]=1221;  a_val[35]=33;  b_val[35]=1;  r_val[35]=132;
+d_val[36]= 323;  a_val[36]=17;  b_val[36]=1;  r_val[36]=34;
+d_val[37]=1365;  a_val[37]=35;  b_val[37]=1;  r_val[37]=140;
+d_val[38]=  10;  a_val[38]= 3;  b_val[38]=1;  r_val[38]=1;
+d_val[39]=1517;  a_val[39]=37;  b_val[39]=1;  r_val[39]=148;
+d_val[40]= 399;  a_val[40]=19;  b_val[40]=1;  r_val[40]=38;
+d_val[41]=1677;  a_val[41]=39;  b_val[41]=1;  r_val[41]=156;
+d_val[42]= 110;  a_val[42]=10;  b_val[42]=1;  r_val[42]=10;
+d_val[43]= 205;  a_val[43]=15;  b_val[43]=1;  r_val[43]=20;
+d_val[44]= 483;  a_val[44]=21;  b_val[44]=1;  r_val[44]=42;
+d_val[45]=2021;  a_val[45]=43;  b_val[45]=1;  r_val[45]=172;
+d_val[46]=  33;  a_val[46]= 6;  b_val[46]=1;  r_val[46]=3;
+/* v1=47	    INVALID */
+d_val[48]=  23;  a_val[48]= 5;  b_val[48]=1;  r_val[48]=2;
+d_val[49]=2397;  a_val[49]=47;  b_val[49]=1;  r_val[49]=188;
+d_val[50]=  39;  a_val[50]= 6;  b_val[50]=1;  r_val[50]=3;
+d_val[51]=  53;  a_val[51]= 7;  b_val[51]=1;  r_val[51]=4;
+/* v1=52	    INVALID */
+d_val[53]=2805;  a_val[53]=51;  b_val[53]=1;  r_val[53]=204;
+d_val[54]= 182;  a_val[54]=13;  b_val[54]=1;  r_val[54]=13;
+d_val[55]=3021;  a_val[55]=53;  b_val[55]=1;  r_val[55]=212;
+d_val[56]=  87;  a_val[56]= 9;  b_val[56]=1;  r_val[56]=6;
+d_val[57]=3245;  a_val[57]=55;  b_val[57]=1;  r_val[57]=220;
+d_val[58]= 210;  a_val[58]=14;  b_val[58]=1;  r_val[58]=14;
+d_val[59]=3477;  a_val[59]=57;  b_val[59]=1;  r_val[59]=228;
+d_val[60]= 899;  a_val[60]=29;  b_val[60]=1;  r_val[60]=58;
+d_val[61]= 413;  a_val[61]=21;  b_val[61]=1;  r_val[61]=28;
+/* v1=62	    INVALID */
+d_val[63]=3965;  a_val[63]=61;  b_val[63]=1;  r_val[63]=244;
+d_val[64]=1023;  a_val[64]=31;  b_val[64]=1;  r_val[64]=62;
+d_val[65]= 469;  a_val[65]=21;  b_val[65]=1;  r_val[65]=28;
+d_val[66]=  17;  a_val[66]= 4;  b_val[66]=1;  r_val[66]=1;
+d_val[67]=4485;  a_val[67]=65;  b_val[67]=1;  r_val[67]=260;
+d_val[68]=1155;  a_val[68]=33;  b_val[68]=1;  r_val[68]=66;
+d_val[69]=4757;  a_val[69]=67;  b_val[69]=1;  r_val[69]=268;
+d_val[70]=  34;  a_val[70]= 6;  b_val[70]=1;  r_val[70]=2;
+d_val[71]=5037;  a_val[71]=69;  b_val[71]=1;  r_val[71]=276;
+d_val[72]=1295;  a_val[72]=35;  b_val[72]=1;  r_val[72]=70;
+d_val[73]= 213;  a_val[73]=15;  b_val[73]=1;  r_val[73]=12;
+d_val[74]=  38;  a_val[74]= 6;  b_val[74]=1;  r_val[74]=2;
+d_val[75]=5621;  a_val[75]=73;  b_val[75]=1;  r_val[75]=292;
+d_val[76]=1443;  a_val[76]=37;  b_val[76]=1;  r_val[76]=74;
+d_val[77]= 237;  a_val[77]=15;  b_val[77]=1;  r_val[77]=12;
+d_val[78]=  95;  a_val[78]=10;  b_val[78]=1;  r_val[78]=5;
+/* v1=79	    INVALID */
+d_val[80]=1599;  a_val[80]=39;  b_val[80]=1;  r_val[80]=78;
+d_val[81]=6557;  a_val[81]=79;  b_val[81]=1;  r_val[81]=316;
+d_val[82]= 105;  a_val[82]=10;  b_val[82]=1;  r_val[82]=5;
+d_val[83]=  85;  a_val[83]= 9;  b_val[83]=1;  r_val[83]=4;
+d_val[84]=1763;  a_val[84]=41;  b_val[84]=1;  r_val[84]=82;
+d_val[85]=7221;  a_val[85]=83;  b_val[85]=1;  r_val[85]=332;
+d_val[86]= 462;  a_val[86]=21;  b_val[86]=1;  r_val[86]=21;
+d_val[87]=7565;  a_val[87]=85;  b_val[87]=1;  r_val[87]=340;
+d_val[88]= 215;  a_val[88]=15;  b_val[88]=1;  r_val[88]=10;
+d_val[89]=7917;  a_val[89]=87;  b_val[89]=1;  r_val[89]=348;
+d_val[90]= 506;  a_val[90]=22;  b_val[90]=1;  r_val[90]=22;
+d_val[91]=8277;  a_val[91]=89;  b_val[91]=1;  r_val[91]=356;
+d_val[92]= 235;  a_val[92]=15;  b_val[92]=1;  r_val[92]=10;
+d_val[93]=8645;  a_val[93]=91;  b_val[93]=1;  r_val[93]=364;
+d_val[94]= 138;  a_val[94]=12;  b_val[94]=1;  r_val[94]=6;
+d_val[95]=9021;  a_val[95]=93;  b_val[95]=1;  r_val[95]=372;
+d_val[96]=  47;  a_val[96]= 7;  b_val[96]=1;  r_val[96]=2;
+d_val[97]=1045;  a_val[97]=33;  b_val[97]=1;  r_val[97]=44;
+/* v1=98	    INVALID */
+d_val[99]=9797;  a_val[99]=97;  b_val[99]=1;  r_val[99]=388;
+d_val[100]=  51; a_val[100]= 7; b_val[100]=1; r_val[100]=2;
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "d_val[100] defined";
+    print "a_val[100] defined";
+    print "b_val[100] defined";
+    print "r_val[100] defined";
+}

lib/mersenne.cal

@@ -0,0 +1,44 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Perform a primality test of 2^p-1, for prime p>1.
+ */
+
+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,10))
+		return 0;
+
+	/* calculate 2^p-1 for later mods */
+	p_mask = 2^p - 1;
+
+	/* lltest: u(i+1) = u(i)^2 - 2 mod 2^p-1 */
+	u = 4;
+	for (i = 2; i < p; ++i) {
+		u = u^2 - 2;
+		u = u&p_mask + u>>p;
+		if (u > p_mask)
+			u = u&p_mask + 1;
+	}
+
+	/* 2^p-1 is prime iff u(p) = 0 mod 2^p-1 */
+	return (u == p_mask);
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "mersenne(p) defined";
+}

lib/mfactor.cal

@@ -0,0 +1,157 @@
+/*
+ * Copyright (c) 1996 Landon Curt Noll
+ *
+ * Permission to use, copy, modify, and distribute this software and
+ * its documentation for any purpose and without fee is hereby granted,
+ * provided that the above copyright, this permission notice and text
+ * this comment, and the disclaimer below appear in all of the following:
+ *
+ *	supporting documentation
+ *	source copies
+ *	source works derived from this source
+ *	binaries derived from this source or from derived source
+ *
+ * LANDON CURT NOLL DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+ * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY 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.
+ *
+ * chongo was here	/\../\		chongo@toad.com
+ */
+
+
+/*
+ * mfactor - find a factor of a Mersenne Number
+ *
+ * Mersenne numbers are numbers of the form:
+ *
+ *	2^n-1	for integer n > 0
+ *
+ * We know that factors of a Mersenne number are of the form:
+ *
+ *	2*k*n+1   and   +/- 1 mod 8
+ *
+ * given:
+ *	n		attempt to factor M(n) = 2^n-1
+ *	start_k		the value k in 2*k*n+1 to start the search
+ *	rept_loop	loop cycle reporting, 0 => none
+ *
+ * returns:
+ *	factor of M(n)
+ */
+define mfactor(n, start_k, rept_loop)
+{
+	local q;	/* test factor 2*k*n+1 */
+	local k;	/* k in 2*k*n+1 */
+	local step2;	/* 2*n */
+	local step6;	/* 6*n */
+	local mod8;	/* q mod 8 */
+	local loop;	/* report loop count */
+
+	/*
+	 * firewall
+	 */
+	if (!isint(n) || n <= 0) {
+		quit "n must be an integer > 0";
+	}
+	if (isnull(start_k)) {
+		start_k = 1;
+	} else if (!isint(start_k) || start_k <= 0) {
+		quit "start_k must be an integer > 0";
+	}
+	if (!isint(rept_loop)) {
+		rept_loop = 0;
+	}
+
+	/*
+	 * setup
+	 */
+	step2 = 2*n;
+	step6 = 6*n;
+	k = start_k - 1;
+	q = 2*k*n+1;
+	/* step2 to the first factor candidate */
+	do {
+		q += step2;
+		mod8 = mod(q,8);
+		++k;
+	} while (mod8 != 1 && mod8 != 7);
+
+	/*
+	 * At this point we are at either at the first or second
+	 * of two consequtive factor candidates depending on if
+	 * the next to k values are 1 and 7 mod 8.
+	 *
+	 * The loops below assume that we will test, bump k by 1
+	 * (move to the 2nd consequtive factor candidate), test and
+	 * bump k by 3 (move to the first of the next consequtive
+	 * factor candidate pair).
+	 *
+	 * In order to prepair, we need to move to the first of
+	 * a consequtive factor candidate pair.  If we happen to
+	 * be on a the 2nd of a pair, we will test it outside
+	 * of the loop and bump to the first of the next pair.
+	 */
+	mod8 = mod(q+step2,8);
+	if (mod8 != 1 && mod8 != 7) {
+		/*
+		 * q is the 2nd of a consequtive factor candidate pair
+		 * so we test q now and bump k by 3.
+		 */
+		if (pmod(2,n,q) == 1) {
+			/* q was a factor afterall, no need to do more! */
+			return q;
+		}
+		q += step6;
+		k += 3;
+	}
+
+	/*
+	 * look for a factor
+	 */
+	loop = k;
+	while (pmod(2,n,q) != 1) {
+
+		/*
+		 * determine if we need to report
+		 */
+		if (rept_loop > 0) {
+			if (rept_loop <= ++loop) {
+				/* report this loop */
+				printf("at 2*%d*%d+1, cpu: %f\n",
+					k, n, runtime());
+				fflush(files(1));
+				loop = 0;
+			}
+			k += 4;
+		}
+
+		/* 
+		 * 1st of a consequtive factor candidate pair is not
+		 * a factor, try the 2nd of that pair
+		 */
+		q += step2;
+		if (pmod(2,n,q) == 1) {
+			break; 	/* factor found */
+		}
+
+		/*
+		 * 2nd of a consequtive factor candidate pair is not
+		 * a factor, try the next pair
+		 */
+		q += step6;
+	}
+
+	/*
+	 * return the factor found
+	 */
+	return q;
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "mfactor(n [, start_k [, rept_loop]])"
+}

lib/mod.cal

@@ -0,0 +1,211 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Routines to handle numbers modulo a specified number.
+ *	a (mod N)
+ */
+
+obj mod {a};			/* definition of the object */
+
+global mod_value = 100;		/* modulus value (value of N) */
+
+
+define mod(a)
+{
+	local obj mod	x;
+
+	if (!isreal(a) || !isint(a))
+		quit "Bad argument for mod function";
+	x.a = a % mod_value;
+	return x;
+}
+
+
+define mod_print(a)
+{
+	if (digits(mod_value) <= 20)
+		print a.a, "(mod", mod_value : ")" :;
+	else
+		print a.a, "(mod N)" :;
+}
+
+
+define mod_one()
+{
+	return mod(1);
+}
+
+
+define mod_cmp(a, b)
+{
+	if (isnum(a))
+		return (a % mod_value) != b.a;
+	if (isnum(b))
+		return (b % mod_value) != a.a;
+	return a.a != b.a;
+}
+
+
+define mod_rel(a, b)
+{
+	if (isnum(a))
+		a = mod(a);
+	if (isnum(b))
+		b = mod(b);
+	if (a.a < b.a)
+		return -1;
+	return a.a != b.a;
+}
+
+
+define mod_add(a, b)
+{
+	local obj mod	x;
+
+	if (isnum(b)) {
+		if (!isint(b))
+			quit "Adding non-integer";
+		x.a = (a.a + b) % mod_value;
+		return x;
+	}
+	if (isnum(a)) {
+		if (!isint(a))
+			quit "Adding non-integer";
+		x.a = (a + b.a) % mod_value;
+		return x;
+	}
+	x.a = (a.a + b.a) % mod_value;
+	return x;
+}
+
+
+define mod_sub(a, b)
+{
+	return a + (-b);
+}
+
+
+define mod_neg(a)
+{
+	local obj mod	x;
+
+	x.a = mod_value - a.a;
+	return x;
+}
+
+
+define mod_mul(a, b)
+{
+	local obj mod	x;
+
+	if (isnum(b)) {
+		if (!isint(b))
+			quit "Multiplying by non-integer";
+		x.a = (a.a * b) % mod_value;
+		return x;
+	}
+	if (isnum(a)) {
+		if (!isint(a))
+			quit "Multiplying by non-integer";
+		x.a = (a * b.a) % mod_value;
+		return x;
+	}
+	x.a = (a.a * b.a) % mod_value;
+	return x;
+}
+
+
+define mod_square(a)
+{
+	local obj mod	x;
+
+	x.a = a.a^2 % mod_value;
+	return x;
+}
+
+
+define mod_inc(a)
+{
+	local x;
+
+	x = a;
+	if (++x.a == mod_value)
+		x.a = 0;
+	return x;
+}
+
+
+define mod_dec(a)
+{
+	local x;
+
+	x = a;
+	if (--x.a < 0)
+		x.a = mod_value - 1;
+	return x;
+}
+
+
+define mod_inv(a)
+{
+	local obj mod	x;
+
+	x.a = minv(a.a, mod_value);
+	return x;
+}
+
+
+define mod_div(a, b)
+{
+	local c, x, y;
+
+	obj mod x, y;
+	if (isnum(a))
+		a = mod(a);
+	if (isnum(b))
+		b = mod(b);
+	c = gcd(a.a, b.a);
+	x.a = a.a / c;
+	y.a = b.a / c;
+	return x * inverse(y);
+}
+
+
+define mod_pow(a, b)
+{
+	local x, y, z;
+
+	obj mod x;
+	y = a;
+	z = b;
+	if (b < 0) {
+		y = inverse(a);
+		z = -b;
+	}
+	x.a = pmod(y.a, z, mod_value);
+	return x;
+}
+
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "obj mod {a} defined";
+    print "mod(a) defined";
+    print "mod_print(a) defined";
+    print "mod_one(a) defined";
+    print "mod_cmp(a, b) defined";
+    print "mod_rel(a, b) defined";
+    print "mod_add(a, b) defined";
+    print "mod_sub(a, b) defined";
+    print "mod_mod(a, b) defined";
+    print "mod_square(a) defined";
+    print "mod_inc(a) defined";
+    print "mod_dec(a) defined";
+    print "mod_inv(a) defined";
+    print "mod_div(a, b) defined";
+    print "mod_pow(a, b) defined";
+    print "mod_value defined";
+    print "set mod_value as needed";
+}

lib/pell.cal

@@ -0,0 +1,74 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Solve Pell's equation; Returns the solution X to: X^2 - D * Y^2 = 1.
+ * Type the solution to pells equation for a particular D.
+ */
+
+define pell(D)
+{
+	local X, Y;
+
+	X = pellx(D);
+	if (isnull(X)) {
+		print "D=":D:" is square";
+		return;
+	}
+	Y = isqrt((X^2 - 1) / D);
+	print X : "^2 - " : D : "*" : Y : "^2 = " : X^2 - D*Y^2;
+}
+
+
+/*
+ * Function to solve Pell's equation
+ * Returns the solution X to:
+ *	X^2 - D * Y^2 = 1
+ */
+define pellx(D)
+{
+	local R, Rp, U, Up, V, Vp, A, T, Q1, Q2, n;
+	local mat ans[2,2];
+	local mat tmp[2,2];
+
+	R = isqrt(D);
+	Vp = D - R^2;
+	if (Vp == 0)
+		return;
+	Rp = R + R;
+	U = Rp;
+	Up = U;
+	V = 1;
+	A = 0;
+	n = 0;
+	ans[0,0] = 1;
+	ans[1,1] = 1;
+	tmp[0,1] = 1;
+	tmp[1,0] = 1;
+	do {
+		T = V;
+		V = A * (Up - U) + Vp;
+		Vp = T;
+		A = U // V;
+		Up = U;
+		U = Rp - U % V;
+		tmp[0,0] = A;
+		ans *= tmp;
+		n++;
+	} while (A != Rp);
+	Q2 = ans[[1]];
+	Q1 = isqrt(Q2^2 * D + 1);
+	if (isodd(n)) {
+		T = Q1^2 + D * Q2^2;
+		Q2 = Q1 * Q2 * 2;
+		Q1 = T;
+	}
+	return Q1;
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "pell(D) defined";
+    print "pellx(D) defined";
+}

lib/pi.cal

@@ -0,0 +1,54 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Calculate pi within the specified epsilon using the quartic convergence
+ * iteration.
+ */
+
+define qpi(epsilon)
+{
+	local niter, yn, ym, tm, an, am, t, tn, sqrt2, epsilon2, count, digits;
+	local bits, bits2;
+
+	if (isnull(epsilon))
+		epsilon = epsilon();
+	digits = digits(1/epsilon);
+	if	(digits <=  8) { niter = 1; epsilon =	1e-8; }
+	else if (digits <= 40) { niter = 2; epsilon =  1e-40; }
+	else if (digits <= 170) { niter = 3; epsilon = 1e-170; }
+	else if (digits <= 693) { niter = 4; epsilon = 1e-693; }
+	else {
+		niter = 4;
+		t = 693;
+		while (t < digits) {
+			++niter;
+			t *= 4;
+		}
+	}
+	epsilon2 = epsilon/(digits/10 + 1);
+	digits = digits(1/epsilon2);
+	sqrt2 = sqrt(2, epsilon2);
+	bits = abs(ilog2(epsilon)) + 1;
+	bits2 = abs(ilog2(epsilon2)) + 1;
+	yn = sqrt2 - 1;
+	an = 6 - 4 * sqrt2;
+	tn = 2;
+	for (count = 0; count < niter; count++) {
+		ym = yn;
+		am = an;
+		tn *= 4;
+		t = sqrt(sqrt(1-ym^4, epsilon2), epsilon2);
+		yn = (1-t)/(1+t);
+		an = (1+yn)^4*am-tn*yn*(1+yn+yn^2);
+		yn = bround(yn, bits2);
+		an = bround(an, bits2);
+	}
+	return (bround(1/an, bits));
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "qpi(epsilon) defined";
+}

lib/pollard.cal

@@ -0,0 +1,35 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Factor using Pollard's p-1 method.
+ */
+
+define factor(N, B, ai, af)
+{
+	local	a, k, i, d;
+
+	if (isnull(B))
+		B = 1000;
+	if (isnull(ai))
+		ai = 2;
+	if (isnull(af))
+		af = ai + 20;
+	k = lcmfact(B);
+	d = lfactor(N, B);
+	if (d > 1)
+		return d;
+	for (a = ai; a <= af; a++) {
+		i = pmod(a, k, N);
+		d = gcd(i - 1, N);
+		if ((d > 1) && (d != N))
+			return d;
+	}
+	return 1;
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "factor(N, B, ai, af) defined";
+}

lib/poly.cal

@@ -0,0 +1,728 @@
+/*
+ * A collection of functions designed for calculations involving
+ * 	polynomials in one variable (by Ernest W. Bowen).
+ *
+ * On starting the program the independent variable has identifier x
+ *	and name "x", i.e. the user can refer to it as x, the
+ *	computer displays it as "x".  The name of the independent
+ *	variable is stored as varname, so, for example, varname = "alpha"
+ *	will change its name to "alpha".  At any time, the independent
+ *	variable has only one name.  For some purposes, a name like
+ *	"sin(t)" or "(a + b)" or "\lambda" might be useful;
+ *	names like "*" or "-27" are legal but might give expressions
+ *	that are difficult to intepret.
+ *
+ * Polynomial expressions may be constructed from numbers and the
+ *	independent variable and other polynomials by the algebraic
+ *	operations +, -, *, ^, and if the result is a polynomial /.
+ *	The operations // and % are defined to have the quotient and
+ *	remainder meanings as usually defined for polynomials.
+ *
+ * When polynomials are assigned to idenfifiers, it is convenient to
+ *	think of the polynomials as values.  For example, p = (x - 1)^2
+ *	assigns to p a polynomial value in the same way as q = (7 - 1)^2
+ *	would assign to q a number value.  As with number expressions
+ *	involving operations, the expression used to define the
+ *	polynomial is usually lost; in the above example, the normal
+ *	computer display for p will be  x^2 - 2x + 1.  Different
+ *	identifiers may of course have the same polynomial value.
+ *
+ * The polynomial we think of as a_0 + a_1 * x + ... + a_n * x^n,
+ *	for number coefficients a_0, a_1, ... a_n may also be
+ *	constructed as pol(a_0, a_1, ..., a_n).  Note that here the
+ *	coefficients are to be in ascending power order.  The independent
+ *	variable is pol(0,1), so to use t, say, as an identifier for
+ *	this, one may assign  t = pol(0,1).  To simultaneously specify
+ *	an identifier and a name for the independent variable, there is
+ *	the instruction var, used as in identifier = var(name).  For
+ *	example, to use "t" in the way "x" is initially, one may give
+ *	the instruction  t = var("t").
+ *
+ * There are four parameters pmode, order, iod and ims for controlling
+ *	the format in which polynomials are displayed.
+ *	The parameter pmode may have values "alg" or "list": the
+ *	former gives a display as an algebraic formula, while
+ *	the latter only lists the coefficients.  Whether the terms or
+ *	coefficients are in ascending or descending power order is
+ *	controlled by order being "up" or "down".  If the
+ *	parameter iod (for integer-only display), the polynomial
+ *	is expressed in terms of a polynomial whose coefficients are
+ *	integers with gcd = 1, the leading coefficient having positive
+ *	real part, with where necessary a leading multiplying integer,
+ *	a Gaussian integer multiplier if the coefficients are complex
+ *	with a common complex factor, and a trailing divisor integer.
+ *	If a non-zero value is assigned to the parameter ims,
+ *	multiplication signs will be inserted where appropriate;
+ *	this may be useful if the expression is to be copied to a
+ *	program or a string to be used with eval.
+ *
+ * For evaluation of polynomials the standard function is ev(p, t).
+ *	If p is a polynomial and t anything for which the relevant
+ *	operations can be performed, this returns the value of p
+ *	at t.  The function ev(p, t) also accepts lists or matrices
+ *	as possible values for p; each element of p is then evaluated
+ *	at t.  For other p, t is ignored and the value of p is returned.
+ *	If an identifier, a, say, is used for the polynomial, list or
+ *	matrix p, the definition
+ *			define a(t) = ev(a, t);
+ *	permits a(t) to be used for the value of a at t as if the
+ *	polynomial, list or matrix were a function.  For example,
+ *	if a = 1 + x^2, a(2) will return the value 5, just as if
+ *			define a(t) = 1 + t^2;
+ *	had been used.   However, when the polynomial definition is
+ *	used, changing the polynomial a will change a(t) to the value
+ *	of the new polynomial at t.  For example,
+ *	after
+ *		L = list(x, x^2, x^3, x^4);
+		define a(t) = ev(a, t);
+ *	the loop
+ *		for (i = 0; i < 4; i++)
+ *			print ev(L[[i]], 5);
+ *	may be replaced by
+ *		for (i = 0; i < 4; i++) {
+ *			a = L[[i]];
+ *			print a(5);
+ *		}
+ *
+ * Matrices with polynomial elements may be added, subtracted and
+ *	multiplied as long as the usual rules for compatibility are
+ *	observed.  Also, matrices may be multiplied by polynomials,
+ *	i.e. if p is a 	polynomial and A a matrix whose elements
+ *	may be numbers or polynomials, p * A returns the matrix of
+ *	the same shape as A with each element multiplied by p.
+ *	Square matrices may also be 'substituted for the variable' in
+ *	polynomials, e.g. if A is an m x m matrix, and
+ *	p = x^2 + 3 * x + 2, ev(p, A) returns the same as
+ *	A^2 + 3 * A + 2 * I, where I is the unit m x m matrix.
+ *
+ * On starting this program, three demonstration polynomials a, b, c
+ *	have been defined.  The functions a(t), b(t), c(t) corresponding
+ *	to a, b, c, and x(t) corresponding to x, have also been
+ *	defined, so the usual function notation can be used for
+ *	evaluations of a, b, c and x.  For x, as long as x identifies
+ *	the independent variable, x(t) should return the value of t,
+ *	i.e. it acts as an identity function.
+ *
+ * Functions defined include:
+ *
+ *	monic(a) returns the monic multiple of a, i.e., if a != 0,
+ *		the multiple of	a with leading coefficient 1
+ *	conj(a) returns the complex conjugate of a
+ *	ispmult(a,b) returns 1 or 0 according as a is or is not
+ *		a polynomial multiple of b
+ *	pgcd(a,b) returns the monic gcd of a and b
+ *	pfgcd(a,b) returns a list of three polynomials (g, u, v)
+ *		where g = pgcd(a,b) and g = u * a + v * b.
+ *	plcm(a,b) returns the monic lcm of a and b
+ *
+ *	interp(X,Y,t) returns the value at t of the polynomial given
+ *		by Newtonian divided difference interpolation, where
+ *		X is a list of x-values, Y a list of corresponding
+ *		y-values.  If t is omitted, the interpolating
+ *		polynomial is returned.  A y-value may be replaced by
+ *		list (y, y_1, y_2, ...), where y_1, y_2, ... are
+ *		the reduced derivatives at the corresponding x;
+ *		i.e. y_r is the r-th derivative divided by fact(r).
+ *	mdet(A) returns the determinant of the square matrix A,
+ *		computed by an algorithm that does not require
+ *		inverses;  the built-in det function usually fails
+ *		for matrices with polynomial elements.
+ *	D(a,n) returns the n-th derivative of a; if n is omitted,
+ *		the first derivative is returned.
+ *
+ * A first-time user can see what the initially defined polynomials
+ *	a, b and c are, and experiment with the algebraic operations
+ *	and other functions that have been defined by giving
+ *	instructions like:
+ *			a
+ *			b
+ *			c
+ *			(x^2 + 1) * a
+ *			a^27
+ *			a * b
+ *			a % b
+ *			a // b
+ *			a(1 + x)
+ *			a(b)
+ *			conj(c)
+ *			g = pgcd(a, b)
+ *			g
+ *			a / g
+ *			D(a)
+ *			mat A[2,2] = {1 + x, x^2, 3, 4*x}
+ *			mdet(A)
+ *			D(A)
+ *			A^2
+ *			define A(t) = ev(A, t)
+ *			A(2)
+ *			A(1 + x)
+ *			define L(t) = ev(L, t)
+ *			L = list(x, x^2, x^3, x^4)
+ *			L(5)
+ *			a(L)
+ *			interp(list(0,1,2,3), list(2,3,5,7))
+ *			interp(list(0,1,2), list(0,list(1,0),2))
+ *
+ * One check on some of the functions is provided by the Cayley-Hamilton
+ *	theorem:  if A is any m x m matrix and I the m x m unit matrix,
+ *	and x is pol(0,1),
+ *			ev(mdet(x * I - A), A)
+ *	should return the zero m x m matrix.
+ */
+
+obj poly {p};
+
+define pol() {
+	local u,i,s;
+	obj poly u;
+	s = list();
+	for (i=1; i<= param(0); i++) append (s,param(i));
+	i=size(s) -1;
+	while (i>=0 && s[[i]]==0) {i--; remove(s)}
+	u.p = s;
+	return u;
+}
+
+define ispoly(a) {
+	local y;
+	obj poly y;
+	return istype(a,y);
+}
+
+define findlist(a) {
+	if (ispoly(a)) return a.p;
+	if (a) return list(a);
+	return list();
+}
+
+pmode = "alg";	/* The other acceptable pmode is "list" */
+ims = 0;	/* To be non-zero if multiplication signs to be inserted */
+iod = 0;	/* To be non-zero for integer-only display */
+order = "down"	/* Determines order in which coefficients displayed */
+
+define poly_print(a) {
+	local f, g, t;
+	if (size(a.p) == 0) {
+		print 0:;
+		return;
+	}
+	if (iod) {
+		g = gcdcoeffs(a);
+		t = a.p[[size(a.p) - 1]] / g;
+		if (re(t) < 0) { t = -t; g = -g;}
+		if (g != 1) {
+			if (!isreal(t)) {
+				if (im(t) > re(t)) g *= 1i;
+				else if (im(t) <= -re(t)) g *= -1i;
+			}
+			if (isreal(g)) f = g;
+			else f = gcd(re(g), im(g));
+			if (num(f) != 1) {
+				print num(f):;
+				if (ims) print"*":;
+			}
+			if (!isreal(g)) {
+				printf("(%d)", g/f);
+				if (ims) print"*":;
+			}
+			if (pmode == "alg") print"(":;
+			polyprint(1/g * a);
+			if (pmode == "alg") print")":;
+			if (den(f) > 1) print "/":den(f):;
+			return;
+		}
+	}
+	polyprint(a);
+}
+
+define polyprint(a) {
+	local s,n,i,c;
+	s = a.p;
+	n=size(s) - 1;
+	if (pmode=="alg") {
+		if (order == "up") {
+			i = 0;
+			while (!s[[i]]) i++;
+			pterm (s[[i]], i);
+			for (i++ ; i <= n; i++) {
+				c = s[[i]];
+				if (c) {
+					if (isreal(c)) {
+						if (c > 0) print" + ":;
+						else {
+							print" - ":;
+							c = -c;
+						}
+					}
+					else print " + ":;
+					pterm(c,i);
+				}
+			}
+			return;
+		}
+		if (order == "down") {
+			pterm(s[[n]],n);
+			for (i=n-1; i>=0; i--) {
+				c = s[[i]];
+				if (c) {
+					if (isreal(c)) {
+						if (c > 0) print" + ":;
+						else {
+							print" - ":;
+							c = -c;
+						}
+					}
+					else print " + ":;
+					pterm(c,i);
+				}
+			}
+			return;
+		}
+		quit "order to be up or down";
+	}
+	if (pmode=="list") {
+		plist(s);
+		return;
+	}
+	print pmode,:"is unknown mode";
+}
+
+
+define poly_neg(a) {
+	local s,i,y;
+	obj poly y;
+	s = a.p;
+	for (i=0; i< size(s); i++) s[[i]] = -s[[i]];
+	y.p = s;
+	return y;
+}
+
+define poly_conj(a) {
+	local s,i,y;
+	obj poly y;
+	s = a.p;
+	for (i=0; i < size(s); i++) s[[i]] = conj(s[[i]]);
+	y.p = s;
+	return y;
+}
+
+define poly_inv(a) = pol(1)/a;	/* This exists only for a of zero degree */
+
+define poly_add(a,b) {
+	local sa, sb, i, y;
+	obj poly y;
+	sa=findlist(a); sb=findlist(b);
+	if (size(sa) > size(sb)) swap(sa,sb);
+	for (i=0; i< size(sa); i++) sa[[i]] += sb[[i]];
+	while (i < size(sb)) append (sa, sb[[i++]]);
+	while (i > 0 && sa[[--i]]==0) remove (sa);
+	y.p = sa;
+	return y;
+}
+
+define poly_sub(a,b) {
+	 return a + (-b);
+}
+
+define poly_cmp(a,b) {
+	local sa, sb;
+	sa = findlist(a);
+	sb=findlist(b);
+	return  (sa != sb);
+}
+
+define poly_mul(a,b) {
+	local sa,sb,i, j, y;
+	if (ismat(a)) swap(a,b);
+	if (ismat(b)) {
+		y = b;
+		for (i=matmin(b,1); i <= matmax(b,1); i++)
+			for (j = matmin(b,2); j<= matmax(b,2); j++)
+				y[i,j] = a * b[i,j];
+		return y;
+	}
+	obj poly y;
+	sa=findlist(a); sb=findlist(b);
+	y.p = listmul(sa,sb);
+	return y;
+}
+
+define listmul(a,b) {
+	local da,db, s, i, j, u;
+	da=size(a)-1; db=size(b)-1;
+	s=list();
+	if (da >= 0 && db >= 0) {
+		for (i=0; i<= da+db; i++) { u=0;
+			for (j = max(0,i-db); j <= min(i, da); j++)
+			u += a[[j]]*b[[i-j]]; append (s,u);}}
+	return s;
+}
+
+define ev(a,t) {
+	local v, i, j;
+	if (ismat(a)) {
+		v = a;
+		for (i = matmin(a,1); i <= matmax(a,1); i++)
+			for (j = matmin(a,2); j <= matmax(a,2); j++)
+				v[i,j] = ev(a[i,j], t);
+		return v;
+	}
+	if (islist(a)) {
+		v = list();
+		for (i = 0; i < size(a); i++)
+			append(v, ev(a[[i]], t));
+		return v;
+	}
+	if (!ispoly(a)) return a;
+	if (islist(t)) {
+		v = list();
+		for (i = 0; i < size(t); i++)
+			append(v, ev(a, t[[i]]));
+		return v;
+	}
+	if (ismat(t)) return evpm(a.p, t);
+	return evp(a.p, t);
+}
+
+define evp(s,t) {
+	local n,v,i;
+	n = size(s);
+	if (!n) return 0;
+	v = s[[n-1]];
+	for (i = n - 2; i >= 0; i--) v=t * v +s[[i]];
+	return v;
+}
+
+define evpm(s,t) {
+	local m, n, V, i, I;
+	n = size(s);
+	m = matmax(t,1) - matmin(t,1);
+	if (matmax(t,2) - matmin(t,2) != m) quit "Non-square matrix";
+	mat V[m+1, m+1];
+	if (!n) return V;
+	mat I[m+1, m+1];
+	matfill(I, 0, 1);
+	V = s[[n-1]] * I;
+	for (i = n - 2; i >= 0; i--) V = t * V + s[[i]] * I;
+	return V;
+}
+pzero = pol(0);
+x = pol(0,1);
+varname = "x";
+define x(t) = ev(x, t);
+
+define iszero(a) {
+	if (ispoly(a))
+		return !size(a.p);
+	return a == 0;
+}
+
+define isstring(a) = istype(a, " ");
+
+define var(name) {
+	if (!isstring(name)) quit "Argument of var is to be a string";
+	varname = name;
+	return pol(0,1);
+}
+
+define pcoeff(a) {
+		if (isreal(a)) print a:;
+		else print "(":a:")":;
+}
+
+define pterm(a,n) {
+	if (n==0) {
+		pcoeff(a);
+		return;
+	}
+	if (n==1) {
+		if (a!=1) {
+			pcoeff(a);
+			if (ims) print"*":;
+		}
+		print varname:;
+		return;
+	}
+	if (a!=1) {
+		pcoeff(a);
+		if (ims) print"*":;
+	}
+	print varname:"^":n:;
+}
+
+define plist(s) {
+	local i, n;
+	n = size(s);
+	print "( ":;
+	if (order == "up") {
+		for (i=0; i< n-1 ; i++)
+			print s[[i]]:",",:;
+		if (n) print s[[i]],")":;
+		else print "0 )":;
+	}
+	else {
+		if (n) print s[[n-1]]:;
+		for (i = n - 2; i >= 0; i--)
+			print ", ":s[[i]]:;
+		print " )":;
+	}
+}
+
+define deg(a) = size(a.p) - 1;
+
+define polydiv(a,b) {
+	local q, r, d, u, i, m, n, sa, sb, sq;
+	obj poly q, r;
+	sa=findlist(a); sb = findlist(b); sq = list();
+	m=size(sa)-1; n=size(sb)-1;
+	if (n<0) quit "Zero divisor";
+	if (m<n) return list(pzero, a);
+	d = sb[[n]];
+	while ( m >= n) { u = sa[[m]]/d;
+		for (i = 0; i< n; i++) sa[[m-n+i]] -= u*sb[[i]];
+		push(sq,u); remove(sa); m--;
+		while (m>=n && sa[[m]]==0) { m--; remove(sa); push(sq,0)}}
+	while (m>=0 && sa[[m]]==0) { m--; remove(sa);}
+	q.p = sq;  r.p = sa;
+	return list(q, r);}
+
+define poly_mod(a,b)  {
+	local u;
+	u=polydiv(a,b);
+	return u[[1]];
+}
+
+define poly_quo(a,b) {
+	local p;
+	p = polydiv(a,b);
+	return p[[0]];
+}
+
+define ispmult(a,b) = iszero(a % b);
+
+define poly_div(a,b) {
+	if (!ispmult(a,b)) quit "Result not a polynomial";
+	return poly_quo(a,b);
+}
+
+define pgcd(a,b) {
+	local r;
+	if (iszero(a) && iszero(b)) return pzero;
+	while (!iszero(b)) {
+		r = a % b;
+		a = b;
+		b = r;
+	}
+	return monic(a);
+}
+
+define plcm(a,b) = monic( a * b // pgcd(a,b));
+
+define pfgcd(a,b) {
+	local u, v, u1, v1, s, q, r, d, w;
+	u = v1 = pol(1); v = u1 = pol(0);
+	while (size(b.p) > 0) {s = polydiv(a,b);
+		q = s[[0]];
+		a = b; b = s[[1]]; u -= q*u1; v -= -q*v1;
+		swap(u,u1); swap(v,v1);}
+	d=size(a.p)-1; if (d>=0 && (w= 1/a.p[[d]]) !=1)
+		 { a *= w; u *= w; v *= w;}
+	return list(a,u,v);
+}
+
+define monic(a) {
+	local s, c, i, d, y;
+	if (iszero(a)) return pzero;
+	obj poly y;
+	s = findlist(a);
+	d = size(s)-1;
+	for (i=0; i<=d; i++) s[[i]] /= s[[d]];
+	y.p = s;
+	return y;
+}
+
+define coefficient(a,n) = (n < size(a.p)) ? a.p[[n]] : 0;
+
+define D(a, n) {
+	local i,j,v;
+	if (isnull(n)) n = 1;
+	if (!isint(n) || n < 1) quit "Bad order for derivative";
+ 	if (ismat(a)) {
+		v = a;
+		for (i = matmin(a,1); i <= matmax(a,1); i++)
+			for (j = matmin(a,2); j <= matmax(a,2); j++)
+				v[i,j] = D(a[i,j], n);
+		return v;
+	}
+	if (!ispoly(a)) return 0;
+	return Dp(a,n);
+}
+
+define Dp(a,n) {
+	local i, v;
+	if (n > 1) return Dp(Dp(a, n-1), 1);
+ 	obj poly v;
+	v.p=list();
+	for (i=1; i<size(a.p); i++) append (v.p, i*a.p[[i]]);
+	return v;
+}
+
+
+define cgcd(a,b) {
+	if (isreal(a) && isreal(b)) return gcd(a,b);
+	while (a) {
+		b -= bround(b/a) * a;
+		swap(a,b);
+	}
+	if (re(b) < 0) b = -b;
+	if (im(b) > re(b)) b *= -1i;
+	else if (im(b) <= -re(b)) b *= 1i;
+	return b;
+}
+
+define gcdcoeffs(a) {
+	local s,i,g, c;
+	s = a.p;
+	g=0;
+	for (i=0; i < size(s) && g != 1; i++)
+		if (c = s[[i]]) g = cgcd(g, c);
+	return g;
+}
+
+define interp(X, Y, t) = evalfd(makediffs(X,Y), t);
+
+define makediffs(X,Y) {
+	local U, D, d, x, y, i, j, k, m, n, s;
+	U = D = list();
+	n = size(X);
+	if (size(Y) != n) quit"Arguments to be lists of same size";
+	for (i = n-1; i >= 0; i--) {
+		x = X[[i]];
+		y = Y[[i]];
+		m = size(U);
+		if (isnum(y)) {
+			d = y;
+			for (j = 0; j < m; j++) {
+				d = D[[j]] = (D[[j]]-d)/(U[[j]] - x);
+			}
+			push(U, x);
+			push(D, y);
+		}
+		else {
+			s = size(y);
+			for (k = 0; k < s ; k++) {
+				d = y[[k]];
+				for (j = 0; j < m; j++) {
+					d = D[[j]] = (D[[j]] - d)/(U[[j]] - x);
+				}
+			}
+			for (j=s-1; j >=0; j--) {
+				push(U,x);
+				push(D, y[[j]]);
+			}
+		}
+	}
+	return list(U, D);
+}
+
+define evalfd(T, t) {
+	local U, D, n, i, v;
+	if (isnull(t)) t = pol(0,1);
+	U = T[[0]];
+	D = T[[1]];
+	n = size(U);
+	v = D[[n-1]];
+	for (i = n-2; i >= 0; i--)
+		v = v * (t - U[[i]]) + D[[i]];
+	return v;
+}
+
+
+define mdet(A) {
+	local n, i, j, k, I, J;
+	n = matmax(A,1) - (i = matmin(A,1));
+	if (matmax(A,2) - (j = matmin(A,2)) != n)
+		quit "Non-square matrix for mdet";
+	I = J = list();
+	k = n + 1;
+	while (k--) {
+		append(I,i++);
+		append(J,j++);
+	}
+	return M(A, n+1, I, J);
+}
+
+define M(A, n, I, J) {
+	local v, J0, i, j, j1;
+	if (n == 1) return A[ I[[0]], J[[0]] ];
+	v = 0;
+	i = remove(I);
+	for (j = 0; j < n; j++) {
+		J0 = J;
+		j1 = delete(J0, j);
+		v += (-1)^(n-1+j) * A[i, j1] * M(A, n-1, I, J0);
+	}
+	return v;
+}
+
+define mprint(A) {
+	local i,j;
+	if (!ismat(A)) quit "Argument to be a matrix";
+	for (i = matmin(A,1); i <= matmax(A,1); i++) {
+		for (j = matmin(A,2); j <= matmax(A,2); j++)
+			printf("%8.4d ", A[i,j]);
+		printf("\n");
+	}
+}
+
+obj poly a;
+obj poly b;
+obj poly c;
+
+define a(t) = ev(a,t);
+define b(t) = ev(b,t);
+define c(t) = ev(c,t);
+
+a=pol(1,4,4,2,3,1);
+b=pol(5,16,8,1);
+c=pol(1+2i,3+4i,5+6i);
+
+global lib_debug;
+if (lib_debug >= 0) {
+	print "obj poly {p} defined";
+	print "pol() defined";
+	print "poly_print(a) defined";
+	print "poly_add(a, b) defined";
+	print "poly_sub(a, b) defined";
+	print "poly_mul(a, b) defined";
+	print "poly_div(a, b) defined";
+	print "poly_quo(a,b) defined";
+	print "poly_mod(a,b) defined";
+	print "poly_neg(a) defined";
+	print "poly_conj(a) defined";
+	print "poly_cmp(a,b) defined";
+	print "iszero(a) defined";
+	print "plist(a) defined";
+	print "listmul(a,b) defined";
+	print "ev(a,t) defined";
+	print "evp(s,t) defined";
+	print "ispoly(a) defined";
+	print "isstring(a) defined";
+	print "var(name) defined";
+	print "pcoeff(a) defined";
+	print "pterm(a,n) defined";
+	print "deg(a) defined";
+	print "polydiv(a,b) defined";
+	print "D(a,n) defined";
+	print "Dp(a,n) defined";
+	print "pgcd(a,b) defined";
+	print "plcm(a,b) defined";
+	print "monic(a) defined";
+	print "pfgcd(a,b) defined";
+	print "interp(X,Y,x) defined";
+	print "makediffs(X,Y) defined";
+	print "evalfd(T,x) defined";
+	print "mdet(A) defined";
+	print "M(A,n,I,J) defined";
+	print "mprint(A) defined";
+}

lib/prompt.cal

@@ -0,0 +1,102 @@
+/*
+ * Copyright (c) 1995 Ernest Bowen
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * By: Ernest Bowen <ernie@neumann.une.edu.au>
+ */
+/*
+ * Demonstration of some uses of prompt() and eval().
+ *
+ * adder() simulates a simple adding machine: starting with sum = 0,
+ * each number entered in response to the ? prompt is added to sum
+ * and the result displayed.  Operation of adder() is ended by
+ * entering "end", "exit" or "quit"; "end" returns to the level from
+ * which adder() is called, e.g. with:
+ *
+ *		for (;;) adder()
+ *
+ * entering "end" would start a new edition with sum = 0; "quit" and
+ * "exit" return to the top level.
+ *
+ * Each response to ? is read as
+ * a string terminated by a newline; the statements and expressions
+ * in this string are compiled and evaluated as in function evaluation;
+ * thus the string may include variables, assignments, functions, etc.
+ * as in:
+ *
+ *		2 + 3
+ *		x = 2 + 3, x^3
+ *		x^2
+ *		local x = 2; while (x < 100) x *= 2; x % 100
+ *		x
+ *		exp(2, 1e-5)
+ *		sum
+ *		print sum^2;
+ *		3; print sum^2;
+ *
+ * (Here the second line creates x as a global variable; the local
+ * variable x in the fourth line has no effect on the global x.  In
+ * the last three lines, sum is the sum of numbers already entered, so
+ * the third last line doubles the value of sum.  The value returned
+ * by "print sum^2;" is the null value, so the second last line adds
+ * nothing to sum.  The last line returns the value 3, i.e. the last
+ * non-null value found for the expressions separated by semicolons,
+ * so sum will be increased by 3 after the "print sum^2;" command
+ * is executed.  xxx The terminating semicolon is essential in the
+ * last two lines.  A command like eval("print 7;") is acceptable to
+ * calc but eval("print 7") causes an exit from calc. xxx)
+ *
+ * If the value returned is not a number (e.g. the name of a list or matrix,
+ * or if the string has syntax errors as in "2 + ", in which case the
+ * value returned is an error value), the compile error messages and a
+ * request for another number are displayed.
+ *
+ * Calling showvalues(str) assumes str defines a function of x as in:
+ *
+ *	"sin(x)", "x^2 + 3*x", "exp(x, 1e-5)".
+ *
+ * Values of the function so defined are returned for values of x
+ * entered in reponse to the ? prompt.  Operation is terminated by
+ * entering "end", "exit" or "quit".
+ */
+
+define adder() {
+	global sum = 0;
+	local s, t;
+	for (;;) {
+		s = prompt("? ");
+		if (s == "end")
+			break;
+		t = eval(s);
+		if (!isnum(t)) {
+			print "Please enter a number";
+			continue;
+		}
+		sum += t;
+		print "\t":sum;
+	}
+}
+
+global x;
+
+define showvalues(str) {
+	local s;
+	for (;;) {
+		s = prompt("? ");
+		if (s == "end")
+			break;
+		x = eval(s);
+		if (!isnum(x)) {
+			print "Please enter a number";
+			continue;
+		}
+		print "\t":eval(str);
+	}
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+	print "adder() defined";
+	print "showvalues(str) defined";
+}

lib/psqrt.cal

@@ -0,0 +1,56 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Calculate square roots modulo a prime.
+ *
+ * Returns null if number is not prime or if there is no square root.
+ * The smaller square root is always returned.
+ */
+
+define psqrt(u, p)
+{
+	local	p1, q, n, y, r, v, w, t, k;
+
+	p1 = p - 1;
+	r = lowbit(p1);
+	q = p >> r;
+	t = 1 << (r - 1);
+	for (n = 2; ; n++) {
+		if (ptest(n, 1) == 0)
+			continue;
+		y = pmod(n, q, p);
+		k = pmod(y, t, p);
+		if (k == 1)
+			continue;
+		if (k != p1)
+			return;
+		break;
+	}
+	t = pmod(u, (q - 1) / 2, p);
+	v = (t * u) % p;
+	w = (t^2 * u) % p;
+	while (w != 1) {
+		k = 0;
+		t = w;
+		do {
+			k++;
+			t = t^2 % p;
+		} while (t != 1);
+		if (k == r)
+			return;
+		t = pmod(y, 1 << (r - k - 1), p);
+		y = t^2 % p;
+		v = (v * t) % p;
+		w = (w * y) % p;
+		r = k;
+	}
+	return min(v, p - v);
+}
+
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "psqrt(u, p) defined";
+}

lib/quat.cal

@@ -0,0 +1,216 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Routines to handle quaternions of the form:
+ *	a + bi + cj + dk
+ *
+ * Note: In this module, quaternians are manipulated in the form:
+ *	s + v
+ * Where s is a scalar and v is a vector of size 3.
+ */
+
+obj quat {s, v};		/* definition of the quaternion object */
+
+
+define quat(a,b,c,d)
+{
+	local obj quat	x;
+
+	x.s = isnull(a) ? 0 : a;
+	mat x.v[3];
+	x.v[0] = isnull(b) ? 0 : b;
+	x.v[1] = isnull(c) ? 0 : c;
+	x.v[2] = isnull(d) ? 0 : d;
+	return x;
+}
+
+
+define quat_print(a)
+{
+	print "quat(" : a.s : ", " : a.v[0] : ", " : a.v[1] : ", " : a.v[2] : ")" :;
+}
+
+
+define quat_norm(a)
+{
+	return a.s^2 + dp(a.v, a.v);
+}
+
+
+define quat_abs(a, e)
+{
+	return sqrt(a.s^2 + dp(a.v, a.v), e);
+}
+
+
+define quat_conj(a)
+{
+	local obj quat	x;
+
+	x.s = a.s;
+	x.v = -a.v;
+	return x;
+}
+
+
+define quat_add(a, b)
+{
+	local obj quat	x;
+
+	if (!istype(b, x)) {
+		x.s = a.s + b;
+		x.v = a.v;
+		return x;
+	}
+	if (!istype(a, x)) {
+		x.s = a + b.s;
+		x.v = b.v;
+		return x;
+	}
+	x.s = a.s + b.s;
+	x.v = a.v + b.v;
+	if (x.v)
+		return x;
+	return x.s;
+}
+
+
+define quat_sub(a, b)
+{
+	local obj quat	x;
+
+	if (!istype(b, x)) {
+		x.s = a.s - b;
+		x.v = a.v;
+		return x;
+	}
+	if (!istype(a, x)) {
+		x.s = a - b.s;
+		x.v = -b.v;
+		return x;
+	}
+	x.s = a.s - b.s;
+	x.v = a.v - b.v;
+	if (x.v)
+		return x;
+	return x.s;
+}
+
+
+define quat_inc(a)
+{
+	local	x;
+
+	x = a;
+	x.s++;
+	return x;
+}
+
+
+define quat_dec(a)
+{
+	local	x;
+
+	x = a;
+	x.s--;
+	return x;
+}
+
+
+define quat_neg(a)
+{
+	local obj quat	x;
+
+	x.s = -a.s;
+	x.v = -a.v;
+	return x;
+}
+
+
+define quat_mul(a, b)
+{
+	local obj quat	x;
+
+	if (!istype(b, x)) {
+		x.s = a.s * b;
+		x.v = a.v * b;
+	} else if (!istype(a, x)) {
+		x.s = b.s * a;
+		x.v = b.v * a;
+	} else {
+		x.s = a.s * b.s - dp(a.v, b.v);
+		x.v = a.s * b.v + b.s * a.v + cp(a.v, b.v);
+	}
+	if (x.v)
+		return x;
+	return x.s;
+}
+
+
+define quat_div(a, b)
+{
+	local obj quat	x;
+
+	if (!istype(b, x)) {
+		x.s = a.s / b;
+		x.v = a.v / b;
+		return x;
+	}
+	return a * quat_inv(b);
+}
+
+
+define quat_inv(a)
+{
+	local	x, q2;
+
+	obj quat x;
+	q2 = a.s^2 + dp(a.v, a.v);
+	x.s = a.s / q2;
+	x.v = a.v / (-q2);
+	return x;
+}
+
+
+define quat_scale(a, b)
+{
+	local obj quat	x;
+
+	x.s = scale(a.s, b);
+	x.v = scale(a.v, b);
+	return x;
+}
+
+
+define quat_shift(a, b)
+{
+	local obj quat	x;
+
+	x.s = a.s << b;
+	x.v = a.v << b;
+	if (x.v)
+		return x;
+	return x.s;
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "obj quat {s, v} defined";
+    print "quat(a, b, c, d) defined";
+    print "quat_print(a) defined";
+    print "quat_norm(a) defined";
+    print "quat_abs(a, e) defined";
+    print "quat_conj(a) defined";
+    print "quat_add(a, e) defined";
+    print "quat_sub(a, e) defined";
+    print "quat_inc(a) defined";
+    print "quat_dec(a) defined";
+    print "quat_neg(a) defined";
+    print "quat_mul(a, b) defined";
+    print "quat_div(a, b) defined";
+    print "quat_inv(a) defined";
+    print "quat_scale(a, b) defined";
+    print "quat_shift(a, b) defined";
+}

lib/randbitrun.cal

@@ -0,0 +1,119 @@
+/*
+ * randbitrun - check rand bit run lengths
+ *
+ * We will use randbit(1) to generate a stream if single bits.
+ * The odds that we will have n bits the same in a row is 1/2^n.
+ */
+/*
+ * Copyright 1995 by Landon Curt Noll.  All Rights Reserved.
+ *
+ * Permission to use, copy, modify, and distribute this software and
+ * its documentation for any purpose and without fee is hereby granted,
+ * provided that the above copyright, this permission notice, and the
+ * disclaimer below appear in all of the following:
+ *
+ * 	* supporting documentation
+ * 	* source copies
+ * 	* source works derived from this source
+ * 	* binaries derived from this source or from derived source
+ *
+ * LANDON CURT NOLL DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+ * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY 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.
+ */
+
+define randbitrun(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 = 18;		/* 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 = randbit(1);	/* our first number */
+    run = 1;
+
+    /*
+     * compute the run length probabilities
+     *
+     * A bit run length of 'r' occurs with a probability of:
+     *
+     *			1/2^n;
+     */
+    for (i=1; i <= MAX_RUN; ++i) {
+	prob[i] = 1.0/(1<<i);
+    }
+
+    /*
+     * look at a number of random number trials
+     */
+    for (i=0; i < run_cnt; ++i) {
+
+	/* get our current number */
+	last = current;
+	current = randbit(1);
+
+	/* 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 = randbit(1);
+	    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 runbit 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);
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "randbitrun([run_length]) defined";
+}

lib/randmprime.cal

@@ -0,0 +1,137 @@
+/*
+ * randmprime - generate a random prime of the form h*2^n-1
+ *
+ * Copyright (c) 1995 by Landon Curt Noll.  All Rights Reserved.
+ *
+ * Permission to use, copy, modify, and distribute this software and
+ * its documentation for any purpose and without fee is hereby granted,
+ * provided that the above copyright, this permission notice and text
+ * this comment, and the disclaimer below appear in all of the following:
+ *
+ *	supporting documentation
+ *	source copies
+ *	source works derived from this source
+ *	binaries derived from this source or from derived source
+ *
+ * LANDON CURT NOLL DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+ * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY 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.
+ *
+ * chongo was here	/\../\		chongo@toad.com
+ */
+
+/* obtain our required libs */
+read -once "cryrand.cal"
+read -once "lucas.cal"
+
+/*
+ * randmprime - find a random prime of the form h*2^n-1 of a given size
+ *
+ * given:
+ *	bits	minimum bits in prime to return
+ *	seed	random seed for scryrand()
+ *	[dbg]	if given, enable debugging
+ *
+ * returns:
+ *	a prime of the form h*2^n-1
+ */
+define
+randmprime(bits, seed, dbg)
+{
+    local n;		/* n as in h*2^n-1 */
+    local h;		/* h as in h*2^n-1 */
+    local plush;	/* value added to h since the beginning */
+    local init;		/* initial cpu time */
+    local start;	/* cpu time before last test */
+    local stop;		/* cpu time afte last test */
+    local tmp;		/* just a tmp place holder value */
+    local ret;		/* h*2^n-1 that is prime */
+
+    /* firewall */
+    if (param(0) < 2 || param(0) > 3) {
+	quit "bad usage: rndprime(dig, seed [,dbg])";
+    }
+    if (!isint(bits) || bits < 0 || !isint(seed) || seed < 0) {
+	quit "args must be non-negative integers";
+    }
+    if (bits < 1) {
+	bits = 1;
+    }
+    if (param(0) == 2 || dbg < 0) {
+    	dbg = 0;
+    }
+
+    /* seed generator */
+    tmp = scryrand(seed);
+
+    /* determine initial h and n values */
+    n = random(bits>>1, highbit(bits)+bits>>1+1);
+    h = cryrand(n);
+    h += iseven(h);
+    while (highbit(h) >= n) {
+       ++n;
+    }
+    if (dbg >= 1) {
+	print "DEBUG3: initial h =", h;
+	print "DEBUG3: initial n =", n;
+    }
+
+    /*
+     * loop until we find a prime
+     */
+    if (dbg >= 1) {
+	start = runtime();
+	init = runtime();
+	plush = 0;
+	print "DEBUG1: testing (h+" : plush : ")*2^" : n : "-1";
+    }
+    while (lucas(h,n) == 0) {
+
+	/* bump h, and n if needed */
+	if (dbg >= 2) {
+	    stop = runtime();
+	    print "DEBUG2: last test:", stop-start, "   total time:", stop-init;
+	}
+	if (dbg >= 1) {
+	    print "DEBUG1: composite: (h+" : plush : ")*2^" : n : "-1";
+	    plush += 2;
+	}
+	h += 2;
+	while (highbit(h) >= n) {
+	   ++n;
+	}
+	if (dbg >= 1) {
+	    print "DEBUG1: testing (h+" : plush : ")*2^" : n : "-1";
+	    start = stop;
+	}
+    }
+
+    /* found a prime */
+    if (dbg >= 2) {
+	stop = runtime();
+	print "DEBUG2: last test:", stop-start, "   total time:", stop-init;
+	print "DEBUG3: " : h : "*2^" : n : "-1 is prime";
+    }
+    if (dbg >= 1) {
+	print "DEBUG1: prime: (h+" : plush : ")*2^" : n : "-1";
+    }
+    ret = h*2^n-1;
+    if (dbg >= 3) {
+	print "DEBUG3: highbit(h):", highbit(h);
+	print "DEBUG3: digits(h):", digits(h);
+	print "DEBUG3: highbit(n):", highbit(n);
+	print "DEBUG3: digits(2^n):", int(n*ln(10)/ln(2)+1);
+	print "DEBUG3: highbit(h*2^n-1):", highbit(ret);
+	print "DEBUG3: digits(h*2^n)-1:", digits(ret);
+    }
+    return ret;
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "randmprime(bits, seed [,dbg]) defined";
+}

lib/randrun.cal

@@ -0,0 +1,128 @@
+/*
+ * randrun - perform a run test on rand()
+ *
+ * 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.
+ */
+/*
+ * Copyright 1995 by Landon Curt Noll.  All Rights Reserved.
+ *
+ * Permission to use, copy, modify, and distribute this software and
+ * its documentation for any purpose and without fee is hereby granted,
+ * provided that the above copyright, this permission notice, and the
+ * disclaimer below appear in all of the following:
+ *
+ * 	* supporting documentation
+ * 	* source copies
+ * 	* source works derived from this source
+ * 	* binaries derived from this source or from derived source
+ *
+ * LANDON CURT NOLL DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+ * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY 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.
+ */
+
+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);
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "randrun([run_length]) defined";
+}

lib/seedrandom.cal

@@ -0,0 +1,136 @@
+/*
+ * Copyright (c) 1996 Landon Curt Noll
+ *
+ * Permission to use, copy, modify, and distribute this software and
+ * its documentation for any purpose and without fee is hereby granted,
+ * provided that the above copyright, this permission notice and text
+ * this comment, and the disclaimer below appear in all of the following:
+ *
+ *	supporting documentation
+ *	source copies
+ *	source works derived from this source
+ *	binaries derived from this source or from derived source
+ *
+ * LANDON CURT NOLL DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+ * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY 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.
+ *
+ * chongo was here	/\../\		chongo@toad.com
+ */
+
+global lib_debug;	/* 1 => print debug statements */
+
+/*
+ * seedrandom - seed the cryptographically strong Blum generator
+ *
+ * This function will seed the random() generator using a method
+ * similar to method suggested for the paranoid in the zrand.c source
+ * file and random help file.
+ *
+ * given:
+ *	seed1 - a large random value (at least 10^20 and perhaps < 10^93)
+ *	seed2 - a large random value (at least 10^20 and perhaps < 10^93)
+ *	size - min Blum modulus as a power of 2 (at least 100, perhaps > 1024)
+ *	trials - number of ptest() trials (default 25)
+ *
+ * returns:
+ *	the previous random state
+ *
+ * NOTE: The [10^20, 10^93) range comes from [2^64, 2^64*fact(55)) range
+ *	 where seeds are effective for srand().  All we really need to
+ *	 do is to insist that a seed is > 2^64, which the 10^20 limit does.
+ */
+define seedrandom(seed1, seed2, size, trials)
+{
+	local p;		/* first Blum prime */
+	local fp;		/* prime co-factor of p-1 */
+	local sp;	 	/* min bit size of p */
+	local q;		/* second Blum prime */
+	local fq;		/* prime co-factor of q-1 */
+	local sq;	 	/* min bit size of q */
+	local n;		/* Blum modulus */
+	local binsize;		/* smallest power of 2 > n=p*q */
+	local r;		/* initial quadratic residue */
+	local rand_state;	/* the initial rand state */
+	local rand_junk;	/* rand state that is not needed */
+	local old_state;	/* old random state to return */
+	local random_cfg;	/* old srandom configuration value */
+
+	/*
+	 * firewall
+	 */
+	if (!isint(seed1)) {
+		quit "1st arg (seed1) is not an int";
+	}
+	if (!isint(seed2)) {
+		quit "2nd arg (seed2) is not an int";
+	}
+	if (!isint(size)) {
+		quit "3rd arg (size) is not an int";
+	}
+	if (!isint(trials)) {
+		trials = 25;
+	}
+	if (digits(seed1) <= 20) {
+		quit "1st arg (seed1) must be > 10^20 and perhaps < 10^93";
+	}
+	if (digits(seed2) <= 20) {
+		quit "2nd arg (seed2) must be > 10^20 and perhaps < 10^93";
+	}
+	if (size < 100) {
+		/* 3% of 100 is 2.97 < 3 whereas 3% of 100 is 3 */
+		quit "3rd arg (size) needs to be > 66 (perhaps >= 1024)";
+	}
+	if (trials < 1) {
+		quit "4th arg (trials) must be > 0";
+	}
+
+	/*
+	 * determine the search parameters
+	 */
+	++size;		/* convert power of 2 to bit length */
+	sp = int((size/2)-(size*0.03)+1);
+	sq = size - sp;
+
+	/*
+	 * find the first Blum prime
+	 */
+	rand_state = srand(seed1);
+	do {
+		fp = nextcand(2^sp+randbit(sp), trials, 0, 3, 4);
+		p = 2*fp+1;
+	} while (ptest(p,trials) == 0);
+
+	/*
+	 * find the 2nd Blum prime
+	 */
+	rand_junk = srand(seed2);
+	do {
+		fq = nextcand(2^sq+randbit(sq), trials, 0, 3, 4);
+		q = 2*fq+1;
+	} while (ptest(q,trials) == 0);
+
+	/*
+	 * seed the Blum generator
+	 */
+	n = p*q;				/* the Blum modulus */
+ 	binsize = higbbit(n)+1;			/* smallest power of 2 > p*q */
+	r = pmod(rand(1<<ceil(binsize*4/5), 1<<(binsize-2)), 2, n);
+	random_cfg = config("srandom", 0);	/* no checks are needed */
+	old_state = srandom(r, n);
+
+	/*
+	 * restore other states that we altered
+	 */
+	rand_junk = srand(rand_state);
+	rand_junk = config("srandom", random_cfg);
+
+	/*
+	 * return the previous random state
+	 */
+	return old_state;
+}

lib/solve.cal

@@ -0,0 +1,48 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Solve the equation f(x) = 0 to within the desired error value for x.
+ * The function 'f' must be defined outside of this routine, and the low
+ * and high values are guesses which must produce values with opposite signs.
+ */
+
+define solve(low, high, epsilon)
+{
+	local flow, fhigh, fmid, mid, places;
+
+	if (isnull(epsilon))
+		epsilon = epsilon();
+	if (epsilon <= 0)
+		quit "Non-positive epsilon value";
+	places = highbit(1 + int(1/epsilon)) + 1;
+	flow = f(low);
+	if (abs(flow) < epsilon)
+		return low;
+	fhigh = f(high);
+	if (abs(flow) < epsilon)
+		return high;
+	if (sgn(flow) == sgn(fhigh))
+		quit "Non-opposite signs";
+	while (1) {
+		mid = bround(high - fhigh * (high - low) / (fhigh - flow), places);
+		if ((mid == low) || (mid == high))
+			places++;
+		fmid = f(mid);
+		if (abs(fmid) < epsilon)
+			return mid;
+		if (sgn(fmid) == sgn(flow)) {
+			low = mid;
+			flow = fmid;
+		} else {
+			high = mid;
+			fhigh = fmid;
+		}
+	}
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "solve(low, high, epsilon) defined";
+}

lib/sumsq.cal

@@ -0,0 +1,44 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Determine the unique two positive integers whose squares sum to the
+ * specified prime.  This is always possible for all primes of the form
+ * 4N+1, and always impossible for primes of the form 4N-1.
+ */
+
+define ss(p)
+{
+	local a, b, i, p4;
+
+	if (p == 2) {
+		print "1^2 + 1^2 = 2";
+		return;
+	}
+	if ((p % 4) != 1) {
+		print p, "is not of the form 4N+1";
+		return;
+	}
+	if (!ptest(p, min(p-2, 10))) {
+		print p, "is not a prime";
+		return;
+	}
+	p4 = (p - 1) / 4;
+	i = 2;
+	do {
+		a = pmod(i++, p4, p);
+	} while ((a^2 % p) == 1);
+	b = p;
+	while (b^2 > p) {
+		i = b % a;
+		b = a;
+		a = i;
+	}
+	print a : "^2 +" , b : "^2 =" , a^2 + b^2;
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "ss(p) defined";
+}

lib/surd.cal

@@ -0,0 +1,288 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Calculate using quadratic surds of the form: a + b * sqrt(D).
+ */
+
+obj surd {a, b};		/* definition of the surd object */
+
+global surd_type = -1;		/* type of surd (value of D) */
+static obj surd surd__;		/* example surd for testing against */
+
+
+define surd(a,b)
+{
+	local x;
+
+	obj surd x;
+	x.a = a;
+	x.b = b;
+	return x;
+}
+
+
+define surd_print(a)
+{
+	print "surd(" : a.a : ", " : a.b : ")" :;
+}
+
+
+define surd_conj(a)
+{
+	local	x;
+
+	obj surd x;
+	x.a = a.a;
+	x.b = -a.b;
+	return x;
+}
+
+
+define surd_norm(a)
+{
+	return a.a^2 + abs(surd_type) * a.b^2;
+}
+
+
+define surd_value(a, xepsilon)
+{
+	local	epsilon;
+
+	epsilon = xepsilon;
+	if (isnull(epsilon))
+		epsilon = epsilon();
+	return a.a + a.b * sqrt(surd_type, epsilon);
+}
+
+define surd_add(a, b)
+{
+	local obj surd	x;
+
+	if (!istype(b, x)) {
+		x.a = a.a + b;
+		x.b = a.b;
+		return x;
+	}
+	if (!istype(a, x)) {
+		x.a = a + b.a;
+		x.b = b.b;
+		return x;
+	}
+	x.a = a.a + b.a;
+	x.b = a.b + b.b;
+	if (x.b)
+		return x;
+	return x.a;
+}
+
+
+define surd_sub(a, b)
+{
+	local obj surd	x;
+
+	if (!istype(b, x)) {
+		x.a = a.a - b;
+		x.b = a.b;
+		return x;
+	}
+	if (!istype(a, x)) {
+		x.a = a - b.a;
+		x.b = -b.b;
+		return x;
+	}
+	x.a = a.a - b.a;
+	x.b = a.b - b.b;
+	if (x.b)
+		return x;
+	return x.a;
+}
+
+
+define surd_inc(a)
+{
+	local	x;
+
+	x = a;
+	x.a++;
+	return x;
+}
+
+
+define surd_dec(a)
+{
+	local	x;
+
+	x = a;
+	x.a--;
+	return x;
+}
+
+
+define surd_neg(a)
+{
+	local obj surd	x;
+
+	x.a = -a.a;
+	x.b = -a.b;
+	return x;
+}
+
+
+define surd_mul(a, b)
+{
+	local obj surd	x;
+
+	if (!istype(b, x)) {
+		x.a = a.a * b;
+		x.b = a.b * b;
+	} else if (!istype(a, x)) {
+		x.a = b.a * a;
+		x.b = b.b * a;
+	} else {
+		x.a = a.a * b.a + surd_type * a.b * b.b;
+		x.b = a.a * b.b + a.b * b.a;
+	}
+	if (x.b)
+		return x;
+	return x.a;
+}
+
+
+define surd_square(a)
+{
+	local obj surd	x;
+
+	x.a = a.a^2 + a.b^2 * surd_type;
+	x.b = a.a * a.b * 2;
+	if (x.b)
+		return x;
+	return x.a;
+}
+
+
+define surd_scale(a, b)
+{
+	local obj surd	x;
+
+	x.a = scale(a.a, b);
+	x.b = scale(a.b, b);
+	return x;
+}
+
+
+define surd_shift(a, b)
+{
+	local obj surd	x;
+
+	x.a = a.a << b;
+	x.b = a.b << b;
+	if (x.b)
+		return x;
+	return x.a;
+}
+
+
+define surd_div(a, b)
+{
+	local x, y;
+
+	if ((a == 0) && b)
+		return 0;
+	obj surd x;
+	if (!istype(b, x)) {
+		x.a = a.a / b;
+		x.b = a.b / b;
+		return x;
+	}
+	y = b;
+	y.b = -b.b;
+	return (a * y) / (b.a^2 - surd_type * b.b^2);
+}
+
+
+define surd_inv(a)
+{
+	return 1 / a;
+}
+
+
+define surd_sgn(a)
+{
+	if (surd_type < 0)
+		quit "Taking sign of complex surd";
+	if (a.a == 0)
+		return sgn(a.b);
+	if (a.b == 0)
+		return sgn(a.a);
+	if ((a.a > 0) && (a.b > 0))
+		return 1;
+	if ((a.a < 0) && (a.b < 0))
+		return -1;
+	return sgn(a.a^2 - a.b^2 * surd_type) * sgn(a.a);
+}
+
+
+define surd_cmp(a, b)
+{
+	if (!istype(a, surd__))
+		return ((b.b != 0) || (a != b.a));
+	if (!istype(b, surd__))
+		return ((a.b != 0) || (b != a.a));
+	return ((a.a != b.a) || (a.b != b.b));
+}
+
+
+define surd_rel(a, b)
+{
+	local x, y;
+
+	if (surd_type < 0)
+		quit "Relative comparison of complex surds";
+	if (!istype(a, surd__)) {
+		x = a - b.a;
+		y = -b.b;
+	} else if (!istype(b, surd__)) {
+		x = a.a - b;
+		y = a.b;
+	} else {
+		x = a.a - b.a;
+		y = a.b - b.b;
+	}
+	if (y == 0)
+		return sgn(x);
+	if (x == 0)
+		return sgn(y);
+	if ((x < 0) && (y < 0))
+		return -1;
+	if ((x > 0) && (y > 0))
+		return 1;
+	return sgn(x^2 - y^2 * surd_type) * sgn(x);
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "obj surd {a, b} defined";
+    print "surd(a, b) defined";
+    print "surd_print(a) defined";
+    print "surd_conj(a) defined";
+    print "surd_norm(a) defined";
+    print "surd_value(a, xepsilon) defined";
+    print "surd_add(a, b) defined";
+    print "surd_sub(a, b) defined";
+    print "surd_inc(a) defined";
+    print "surd_dec(a) defined";
+    print "surd_neg(a) defined";
+    print "surd_mul(a, b) defined";
+    print "surd_square(a) defined";
+    print "surd_scale(a, b) defined";
+    print "surd_shift(a, b) defined";
+    print "surd_div(a, b) defined";
+    print "surd_inv(a) defined";
+    print "surd_sgn(a) defined";
+    print "surd_cmp(a, b) defined";
+    print "surd_rel(a, b) defined";
+    print "surd_type defined";
+    print "set surd_type as needed";
+}

lib/test1700.cal

@@ -0,0 +1,12 @@
+/*
+ * Copyright (c) 1995 Landon Curt Noll
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * By: Landon Curt Noll
+ *     chongo@toad.com  -or-  ...!{pyramid,sun,uunet}!hoptoad!chongo
+ *
+ * This library is used by the 1700 series of the regress.cal test suite.
+ */
+
+++value;

lib/test2300.cal

@@ -0,0 +1,97 @@
+/*
+ * Copyright (c) 1995 Landon Curt Noll
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * By: Landon Curt Noll
+ *     chongo@toad.com  -or-  ...!{pyramid,sun,uunet}!hoptoad!chongo
+ *
+ * This library is used by the 2300 series of the regress.cal test suite.
+ */
+
+
+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;
+}

lib/test2600.cal

@@ -0,0 +1,516 @@
+/*
+ * Copyright (c) 1996 Ernest Bowen and Landon Curt Noll
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * By: Ernest Bowen and Landon Curt Noll
+ *     ernie@neumann.une.edu.au and chongo@toad.com
+ *
+ * This library is used by the 2600 series of the regress.cal test suite.
+ */
+/*
+ * Stringent tests of some of calc's builtin functions.
+ * Most of the tests are concerned with the accuracy of the value
+ * returned for a function; usually it is expected that
+ * remainder (true value - calculated value) will be less in
+ * absolute value than "epsilon", where this is either a specified
+ * argument eps, or if this is omitted, the current value of epsilon().
+ * In some cases the remainder is to have a particular sign, or to
+ * have absolute value not exceeding eps/2, or in some cases 3 * eps/4.
+ *
+ * Typical of these tests is testpower("power", n, b, eps, verbose).
+ * Here n is the number of numbers a for which power(a, b, eps) is to
+ * be evaluated; the ratio c = (true value - calculated value)/eps
+ * is calculated and if this is not less in absolute value than
+ * 0.75, a "failure" is recorded and the value of a displayed.
+ * On completion of the tests, the minimum and maximum values of
+ * c are displayed.
+ *
+ * The numbers a are usually large "random" integers or sometimes
+ * ratios of such integers.  In some cases the formulae used to
+ * calculate c assume eps is small compared with the value of the
+ * function.  If eps is very small, say 1e-1000, or if the denominator
+ * of b in power(a, b, eps) is large, the computation required for
+ * a test may be very heavy.
+ *
+ * Test funcations are called as:
+ *
+ *	testabc(str, ..., verbose)
+ *
+ * where str is a string that names the test.  This string is printed
+ * without a newline (if verbose > 0), near the beginning of the function.
+ * The verbose parameter controls how verbose the test will be:
+ *
+ *	0 - print nothing
+ *	1 - print str and the error count
+ *	2 - print min and max errors as well
+ *	3 - print everything including individual loop counts
+ *
+ * All functions return the number of errors that they detected.
+ */
+
+global defaultverbose = 1;	/* default verbose value */
+global err;
+
+define testismult(str, n, verbose)
+{
+	local a, b, c, i, m;
+
+	if (isnull(verbose)) verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	for (i = 0; i < n; i++) {
+		if (verbose > 2) print i,:;
+		a = scale(rand(1,1e1000), rand(100));
+		b = scale(rand(1,1e1000), rand(100));
+		c = a * b;
+		if (!ismult(c,a)) {
+			m++;
+			if (verbose > 1) {
+				printf("*** Failure with:\na = %d\nb = %d\n", a,b);
+			}
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+		} else {
+			printf("no errors\n");
+		}
+	}
+	return m;
+}
+
+define testsqrt(str, n, eps, verbose)
+{
+	local a, c, i, x, m, min, max;
+
+	if (isnull(verbose)) verbose = 2;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	min = 1000;
+	max = -1000;
+	if (isnull(eps))
+		eps = epsilon();
+	for (i = 1; i <= n; i++) {
+		if (verbose > 2) print i,:;
+		a = scale(rand(1,1000), rand(100));
+		x = sqrt(a, eps);
+		if (x)
+			c = (a/x - x)/2/eps;
+		else
+			c = a/eps;		/* ??? */
+		if (c < min)
+			min = c;
+		if (c > max)
+			max = c;
+		if (abs(c) > 1) {
+			m++;
+			if (verbose > 1) {
+				printf("*** Failure with:\na = %d\neps = %d\n", a,eps);
+			}
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+			printf("    %s: rem/eps min=%d, max=%d\n",
+			       str, min, max);
+		} else {
+			printf("no errors\n");
+		}
+	}
+	if (verbose > 1) {
+		printf("    %s: rem/eps min=%0.4d, max=%0.4d\n", str, min, max);
+	}
+	return m;
+}
+
+
+define testexp(str, n, eps, verbose)
+{
+	local i, a, c, m, min, max;
+
+	if (isnull(verbose)) verbose = 2;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	if (isnull(eps))
+		eps = epsilon();
+	min = 1000;
+	max = -1000;
+	for (i = 1; i <= n; i++) {
+		if (verbose > 2) print i,:;
+		a = rand(1,1e20)/rand(1,1e20) + rand(50);
+		if (rand(1))
+			a = -a;
+		c = cexp(a, eps);
+		if (c < min)
+			min = c;
+		if (c > max)
+			max = c;
+		if (abs(c) > 0.02) {
+			m++;
+			if (verbose > 1) {
+				printf("*** Failure with:\na = %d\neps = %d\n", a,eps);
+			}
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+			printf("    %s: rem/eps min=%d, max=%d\n",
+			       str, min, max);
+		} else {
+			printf("no errors\n");
+		}
+	}
+	if (verbose > 1) {
+		printf("    %s: rem/eps min=%0.4d, max=%0.4d\n", str, min, max);
+	}
+	return m;
+}
+
+
+define cexp(x,eps)		/* Find relative rem/eps for exp(x, eps) */
+{
+	local eps1, v, v1, c;
+
+	if (isnull(eps))
+		eps = epsilon();
+	eps1 = eps * 1e-6;
+	v = exp(x, eps);
+	v1 = exp(x, eps1);
+	c = round((v1 - v)/v1/eps, 6, 24);
+	return c;
+}
+
+
+define testln(str, n, eps, verbose)
+{
+	local i, a, c, m, min, max;
+
+	if (isnull(verbose)) verbose = 2;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	if (isnull(eps))
+		eps = epsilon();
+	min = 1000;
+	max = -1000;
+	for (i = 1; i <= n; i++) {
+		if (verbose > 2) print i,:;
+		a = rand(1,1e20)/rand(1,1e20) + rand(50);
+		c = cln(a, eps);
+		if (c < min)
+			min = c;
+		if (c > max)
+			max = c;
+		if (abs(c) > 0.5) {
+			m++;
+			if (verbose > 1) {
+				printf("*** Failure with:\na = %d\neps = %d\n", a,eps);
+			}
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+			printf("    %s: rem/eps min=%d, max=%d\n",
+			       str, min, max);
+		} else {
+			printf("no errors\n");
+		}
+	}
+	if (verbose > 1) {
+		printf("    %s: rem/eps min=%0.4d, max=%0.4d\n", str, min, max);
+	}
+	return m;
+}
+
+define cln(a, eps)
+{
+	local eps1, v, v1, c;
+
+	if (isnull(eps))
+		eps = epsilon();
+	eps1 = eps/1e6;
+	v = ln(a, eps);
+	v1 = ln(a, eps1);
+	c = round((v1 - v)/eps, 6, 24);
+	return c;
+}
+
+
+define testpower(str, n, b, eps, verbose)
+{
+	local i, a, c, m, min, max;
+
+	if (isnull(verbose)) verbose = 2;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	if (isnull(eps))
+		eps = epsilon();
+	if (!isnum(b))
+		quit "Second argument (exponent) to be a number";
+	min = 1000;
+	max = -1000;
+	for (i = 1; i <= n; i++) {
+		if (verbose > 2) print i,:;
+		a = rand(1,1e20)/rand(1,1e20);
+		c = cpow(a, b, eps);
+		if (abs(c) > .75) {
+			m++;
+			if (verbose > 1) {
+				printf("*** Failure for a = %d\n", a);
+			}
+		}
+		if (c < min)
+			min = c;
+		if (c > max)
+			max = c;
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+			printf("    %s: rem/eps min=%d, max=%d\n",
+			       str, min, max);
+		} else {
+			printf("no errors\n");
+		}
+	}
+	if (verbose > 1) {
+		printf("    %s: rem/eps min=%0.4d, max=%0.4d\n", str, min, max);
+	}
+	return m;
+}
+
+
+define cpow(a, b, eps)		/* Find rem/eps for power(a,b,eps) */
+{
+	local v, v1, c, n, d, h;
+
+	if (isnull(eps))
+		eps = epsilon();
+	n = num(b);
+	d = den(b);
+
+	v = power(a, b, eps);
+	h = (a^n/v^d - 1) * v/d;
+	c = round(h/eps, 6, 24);
+	return c;
+}
+
+define testgcd(str, n, verbose)
+{
+	local i, a, b, g, m;
+
+	if (isnull(verbose)) verbose = 2;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	for (i = 1; i <= n; i++) {
+		if (verbose > 2) print i,:;
+		a = rand(1,1e1000);
+		b = rand(1,1e1000);
+		g = gcd(a,b);
+		if (!ismult(a,g) || !ismult(b,g) || !g || !isrel(a/g, b/g)) {
+			m++;
+			printf("*** Failure for a = %d, b = %d\n", a, b);
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+		} else {
+			printf("no errors\n");
+		}
+	}
+	return m;
+}
+
+define mkreal() = scale(rand(-1000,1001)/rand(1,1000), rand(-100, 101));
+
+define mkcomplex() = mkreal() + 1i * mkreal();
+
+define mkbigreal()
+{
+	local x;
+
+	x = rand(100, 1000)/rand(1,10);
+	if (rand(2))
+		x = -x;
+	return x;
+}
+
+define mksmallreal() = rand(-10, 11)/rand(100,1000);
+
+define testappr(str, n, verbose)
+{
+	local x, y, z, m, i, p;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	for (i = 1; i <= n; i++) {
+		x = rand(3) ? mkreal(): mkcomplex();
+		y = mkreal();
+		if (verbose > 2)
+			printf("    %d: x = %d, y = %d\n", i, x, y);
+
+		for (z = 0; z < 32; z++) {
+			p = checkappr(x,y,z,verbose);
+			if (p) {
+				printf("*** Failure for x=%d, y=%d, z=%d\n",
+					x, y, z);
+				m++;
+			}
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+		} else {
+			printf("no errors\n");
+		}
+	}
+	return m;
+}
+
+
+define checkappr(x,y,z,verbose)		/* Returns 1 if an error is detected */
+{
+	local a;
+
+	a = appr(x,y,z);
+	if (verbose > 1)
+		printf("\ta = %d\n", a);
+	if (isreal(x))
+		return checkresult(x,y,z,a);
+	if (isnum(x))
+		return checkresult(re(x), y, z, re(a))
+			| checkresult(im(x), y, z, im(a));
+
+	quit "Bad first argument for checkappr()";
+}
+
+define checkresult(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); break;
+	}
+	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);
+	}
+	if (verbose > 1) {
+		if (err) {
+			print "***", err, "error(s) found in test2600";
+		} else {
+			print "no errors in test2600";
+		}
+	}
+	return tnum;
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+	print "global defaultverbose defined";
+	print "global err defined";
+	print "testismult(str,n,verbose) defined";
+	print "testsqrt(str,n,eps,verbose) defined";
+	print "testexp(str,n,eps,verbose) defined";
+	print "testln(str,n,eps,verbose) defined";
+	print "testpower(str,n,b,eps,verbose) defined";
+	print "testgcd(str,n,verbose) defined";
+	print "cpow(x,n,eps) defined";
+	print "cexp(x,eps) defined";
+	print "cln(x,eps) defined";
+	print "mkreal() defined";
+	print "mkcomplex() defined";
+	print "mkbigreal() defined";
+	print "mksmallreal() defined";
+	print "testappr(str,n,verbose) defined";
+	print "checkappr(x,y,z,verbose) defined";
+	print "checkresult(x,y,z,a) defined";
+	print "test2600(verbose,tnum) defined";
+}

lib/test2700.cal

@@ -0,0 +1,331 @@
+/*
+ * Copyright (c) 1996 Ernest Bowen and Landon Curt Noll
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * By: Ernest Bowen and Landon Curt Noll
+ *     ernie@neumann.une.edu.au and chongo@toad.com
+ *
+ * This library is used by the 2700 series of the regress.cal test suite.
+ */
+/*
+ * The following script 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.  The
+ * number generators are essentially as in the script I sent yesterday.
+ * 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.
+ */
+
+global defaultverbose = 1;
+global err;
+
+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;
+
+/*
+ * XXX - Should 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;
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+	print "global defaultverbose defined";
+	print "global err defined";
+	print "mknonnegreal() defined";
+	print "mkposreal() defined";
+	print "mkreal_2700() defined";
+	print "mknonzeroreal() defined";
+	print "mkposfrac() defined";
+	print "mkfrac() defined";
+	print "mksquarereal() defined";
+	print "mknonsquarereal() defined";
+	print "mkcomplex_2700() defined";
+	print "testcsqrt(str,n,verbose) defined";
+	print "checksqrt(x,y,z,v) defined";
+	print "checkavrem(A,B,X,eps) defined";
+	print "checkrounding(s,n,t,u,z) defined";
+	print "iscomsq(x) defined";
+	print "test2700(verbose,tnum) defined";
+}

lib/test3100.cal

@@ -0,0 +1,31 @@
+/*
+ * Copyright (c) 1995 Ernest Bowen and Landon Curt Noll
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * By: Ernest Bowen and Landon Curt Noll
+ *     ernie@neumann.une.edu.au and chongo@toad.com
+ *
+ * This library is used by the 3100 series of the regress.cal test suite.
+ */
+
+obj res {r};
+global md;
+define res_test(a) = !ismult(a.r, md);
+define res_sub(a,b) {local obj res v = {(a.r - b.r) % md}; return v;};
+define res_mul(a,b) {local obj res v = {(a.r * b.r) % md}; return v;};
+define res_neg(a) {local obj res v = {(-a.r) % md}; return v;};
+define res_inv(a) {local obj res v = {minv(a.r, md)}; return v;};
+define res(x) {local obj res v = {x % md}; return v;};
+
+global lib_debug;
+if (lib_debug >= 0) {
+	print "obj res defined";
+	print "global md defined";
+	print "res_test(a) defined";
+	print "res_sub(a, b) defined";
+	print "res_mul(a, b) defined";
+	print "res_neg(a) defined";
+	print "res_inv(a) defined";
+	print "res(x) defined";
+}

lib/test3300.cal

@@ -0,0 +1,134 @@
+/*
+ * Copyright (c) 1995 Ernest Bowen and Landon Curt Noll
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * By: Ernest Bowen and Landon Curt Noll
+ *     ernie@neumann.une.edu.au and chongo@toad.com
+ *
+ * This library is used by the 3300 series of the regress.cal test suite.
+ */
+
+global defaultverbose = 1;	/* default verbose value */
+global err;
+
+define testi(str, n, N, verbose)
+{
+	local A, t, i, j, d1, d2;
+	local m;
+
+	if (isnull(verbose)) verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	if (isnull(N))
+		N = 1e6;
+	mat A[n,n];
+	for (i = 0; i < n; i++)
+		for (j = 0; j < n; j++)
+			A[i,j] = rand(-N, N);
+	t = runtime();
+	d1 = det(A);
+	t = runtime() - t;
+	d2 = det(A^2);
+	if (d2 != d1^2) {
+		if (verbose > 0) {
+			printf("*** Failure for n=%d, N=%d, d1=%d\n", n, N, d1);
+		}
+		return 1;	/* error */
+	} else {
+		if (verbose > 0) {
+			printf("no errors\n");
+		}
+		if (verbose > 1) {
+			printf("ok: n=%d, N=%d, d1=%d, t=%d\n", n, N, d1, t);
+		}
+	}
+	return 0;	/* ok */
+}
+
+define testr(str, n, N, verbose)
+{
+	local A, t, i, j, d1, d2;
+
+	if (isnull(verbose)) verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	if (isnull(N))
+		N = 1e6;
+	mat A[n,n];
+	for (i = 0; i < n; i++)
+		for (j = 0; j < n; j++)
+			A[i,j] = rand(-(N^2), N^2)/rand(1, N);
+	t = runtime();
+	d1 = det(A);
+	t = runtime() - t;
+	d2 = det(A^2);
+	if (d2 != d1^2) {
+		if (verbose > 0) {
+			printf("*** Failure for n=%d, N=%d, d1=%d\n", n, N, d1);
+		}
+		return 1;	/* error */
+	} else {
+		if (verbose > 0) {
+			printf("no errors\n");
+		}
+		if (verbose > 1) {
+			printf("ok: n=%d, N=%d, d1=%d, t=%d\n", n, N, d1, t);
+		}
+	}
+	return 0;	/* ok */
+}
+
+/*
+ * test3300 - perform all of the above tests a bunch of times
+ */
+define test3300(verbose, tnum)
+{
+	local N;	/* test parameter */
+	local i;
+
+	/*
+	 * set test parameters
+	 */
+	if (isnull(verbose)) {
+		verbose = defaultverbose;
+	}
+	N = 1e6;
+	srand(3300e3300);
+
+	/*
+	 * test a lot of stuff
+	 */
+	for (i=0; i < 19; ++i) {
+		err += testi(strcat(str(tnum++), ": testi(", str(i), ")"), \
+			     i, N, verbose);
+	}
+	for (i=0; i < 9; ++i) {
+		err += testr(strcat(str(tnum++), ": testr(", str(i), ")"), \
+			     i, N, verbose);
+	}
+
+	/*
+	 * test results
+	 */
+	if (verbose > 1) {
+		if (err) {
+			print "***", err, "error(s) found in testall";
+		} else {
+			print "no errors in testall";
+		}
+	}
+	return tnum;
+}
+
+global lib_debug;
+
+if (lib_debug >= 0) {
+	print "global defaultverbose defined";
+	print "global err defined";
+	print "testi(str, n, N, verbose) defined";
+	print "testr(str, n, N, verbose) defined";
+	print "test3300(verbose, tnum) defined";
+}

lib/test3400.cal

@@ -0,0 +1,315 @@
+/*
+ * Copyright (c) 1996 Ernest Bowen and Landon Curt Noll
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * By: Ernest Bowen and Landon Curt Noll
+ *     ernie@neumann.une.edu.au and chongo@toad.com
+ *
+ * This library is used by the 3400 series of the regress.cal test suite.
+ */
+/*
+ * 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]]
+ */
+
+global defaultverbose = 1;	/* default verbose value */
+global err;
+
+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;
+}
+
+global lib_debug;
+
+if (lib_debug >= 0) {
+	print "global defaultverbose defined";
+	print "global err defined";
+	print "test3401(str, n, eps, verbose) defined";
+	print "test3402(str, n, eps, verbose) defined";
+	print "test3403(str, n, eps, verbose) defined";
+	print "test3404(str, n, eps, verbose) defined";
+	print "test3405(str, n, eps, verbose) defined";
+	print "test3406(str, n, eps, verbose) defined";
+	print "test3400(verbose, tnum) defined";
+}

lib/test3500.cal

@@ -0,0 +1,286 @@
+/*
+ * Copyright (c) 1996 Ernest Bowen and Landon Curt Noll
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * By: Ernest Bowen and Landon Curt Noll
+ *     ernie@neumann.une.edu.au and chongo@toad.com
+ *
+ * This library is used by the 3500 series of the regress.cal test suite.
+ */
+/*
+ * Stringent tests of the functions frem, fcnt, gcdrem.
+ *
+ * testf(n) gives n tests of frem(x,y) and fcnt(x,y) with randomly
+ * 	integers x and y generated so that x = f * y^k where f, y and
+ *	k are randomly generated.
+ *
+ * testg(n) gives n tests of gcdrem(x,y) with x and y generated as for
+ *	testf(n).
+ *
+ * testh(n,N) gives n tests of g = gcdrem(x,y) where x and y are products of
+ *	powers of small primes some of which are common to both x and y.
+ *	This test uses f = abs(x) and iteratively f = frem(f,p) where
+ *	p varies over the prime divisors of y; the final value for f
+ *	should equal g.  For both x and y the primes are raised to the
+ *	power rand(N); N defaults to 10.
+ *
+ * If verbose is > 1, the numbers x, y and values for some of the
+ *	functions will be displayed.  Numbers used in testf()
+ *	and testg() occasionally have thousands of digits.
+ *
+ */
+
+global defaultverbose = 1;	/* default verbose value */
+global err;
+
+define testfrem(x,y,verbose)
+{
+	local f, n;
+
+	if (isnull(verbose)) verbose = defaultverbose;
+
+	f = frem(x,y);
+	n = fcnt(x,y);
+	if (verbose > 1)
+		printf("frem = %d, fcnt = %d\n\n", f, n);
+	if (abs(x) != f * abs(y)^n)
+		return 1;
+	if (!ismult(x,y) || abs(y) <= 1) {
+		if (f != abs(x))
+			return 2;
+		if (n != 0)
+			return 3;
+		return 0;
+	}
+	if (x == 0) {
+		if (f != 0 || n != 0)
+			return 4;
+		return 0;
+	}
+	if (f < 0 || !isint(f) || n <= 0)
+		return 5;
+	if (ismult(f, y))
+		return 6;
+	if (!ismult(x, y^n))
+		return 7;
+	if (ismult(x, y^(n+1)))
+		return 8;
+	return 0;
+}
+
+define testgcdrem(x,y,verbose)
+{
+	local d, q;
+
+	if (isnull(verbose)) verbose = defaultverbose;
+
+	d = gcdrem(x,y);
+	if (verbose > 1)
+		printf("gcdrem = %d\n\n", d);
+	if (y == 0) {
+		if (d != 1)
+			return 1;
+		return 0;
+	}
+	if (x == 0) {
+		if (d != 0)
+			return 2;
+		return 0;
+	}
+	if (d <= 0)
+		return 3;
+	q = x/d;
+	if (!isint(q))
+		return 4;
+	if (!isrel(d, y))
+		return 5;
+	if (!isrel(d, q))
+		return 6;
+	return 0;
+}
+
+define testf(str,n,verbose)
+{
+	local m, x, y, i, k, y1, f1, f, fail;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	for (i = 0; i < n; i++) {
+		y1 = rand(2^rand(1,6));
+		y = rand(-(2^y1), 1 + 2^y1);
+		f1 = rand(2^rand(1,11));
+		f = rand(-(2^f1), 1+2^f1);
+		k = rand(1,1+2^10);
+		x = f * y^k;
+		if (verbose > 1) {
+			printf("x = %d\n", x);
+			printf("y = %d\n", y);
+		}
+		fail = testfrem(x,y,verbose);
+		if (fail != 0) {
+			printf("*** Failure %d on loop %d\n", fail, i);
+			if (verbose > 1) {
+				printf("    x = %d\n", x);
+				printf("    y = %d\n", y);
+			}
+			m++;
+		}
+	}
+
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+		} else {
+			printf("no errors\n");
+		}
+	}
+	return m;
+}
+
+
+define testg(str,n,verbose)
+{
+	local m, x, y, i, k, y1, f1, f, fail;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	for (i = 0; i < n; i++) {
+		y1 = rand(2^rand(1,6));
+		y = rand(-(2^y1), 1 + 2^y1);
+		f1 = rand(2^rand(1,11));
+		f = rand(-(2^f1), 1+2^f1);
+		k = rand(1,1+2^10);
+		x = f * y^k;
+		if (verbose > 1) {
+			printf("x = %d\n", x);
+			printf("y = %d\n", y);
+		}
+		fail = testgcdrem(x,y,verbose);
+		if (fail != 0) {
+			printf("*** Failure %d on loop %d\n", fail, i);
+			if (verbose > 1) {
+				printf("    x = %d\n", x);
+				printf("    y = %d\n", y);
+			}
+			m++;
+		}
+	}
+
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+		} else {
+			printf("no errors\n");
+		}
+	}
+	return m;
+}
+
+define testh(str,n,N,verbose)
+{
+	local m, i, x, y, f, g;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	if (isnull(N))
+		N = 61;
+	for (i = 0; i < n; i ++) {
+		x = 2^rand(N)*3^rand(N) * 7^rand(N) * 11^rand(N) * 101^rand(N);
+		y = 2^rand(N) * 3^rand(N) * 5^rand(N) * 11^rand(N) * 53^rand(N);
+		if (rand(2)) x = -x;
+		if (rand(2)) y = -y;
+
+		if (verbose > 1) {
+			printf("x = %d\n", x);
+			printf("y = %d\n", y);
+		}
+		f = abs(x);
+		g = gcdrem(x,y);
+		if (ismult(y,2)) f = frem(f,2);
+		if (ismult(y,3)) f = frem(f,3);
+		if (ismult(y,5)) f = frem(f,5);
+		if (ismult(y,11)) f = frem(f,11);
+		if (ismult(y,53)) f = frem(f,53);
+
+		if (f != g) {
+			printf("*** Failure on loop %d\n", i);
+			if (verbose > 1) {
+				printf("    x = %d\n", x);
+				printf("    y = %d\n", y);
+				printf("    f = %d\n", f);
+				printf("    g = %d\n", g);
+			}
+			m++;
+		}
+	}
+
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+		} else {
+			printf("no errors\n");
+		}
+	}
+	return m;
+}
+
+/*
+ * test3500 - perform all of the above tests a bunch of times
+ */
+define test3500(verbose, tnum, n, N)
+{
+	/* set test parameters */
+	if (isnull(verbose)) {
+		verbose = defaultverbose;
+	}
+	if (isnull(tnum)) {
+		tnum = 3501;	/* default test number */
+	}
+	if (isnull(n)) {
+		n = 200;
+	}
+	if (isnull(N)) {
+		N = 61;
+	}
+
+	/*
+	 * test a lot of stuff
+	 */
+	srand(3500e3500);
+	err += testf(strcat(str(tnum++), ": frem/fcnt"), n, verbose);
+	err += testg(strcat(str(tnum++), ": gcdrem"), n, verbose);
+	err += testh(strcat(str(tnum++),": gcdrem #2"), n, N, verbose);
+	if (verbose > 1) {
+		if (err) {
+			print "***", err, "error(s) found in testall";
+		} else {
+			print "no errors in testall";
+		}
+	}
+	return tnum;
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+	print "global defaultverbose defined";
+	print "global err defined";
+	print "testfrem(x, y, verbose) defined";
+	print "testgcdrem(x, y, verbose) defined";
+	print "testf(str, n, verbose) defined";
+	print "testg(str, n, verbose) defined";
+	print "testh(str, n, N, verbose) defined";
+	print "test3500(verbose, n, N) defined";
+}

lib/test4000.cal

@@ -0,0 +1,485 @@
+/*
+ * Copyright (c) 1996 Ernest Bowen and Landon Curt Noll
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * By: Ernest Bowen and Landon Curt Noll
+ *     ernie@neumann.une.edu.au and chongo@toad.com
+ *
+ * This library is used by the 4000 series of the regress.cal test suite.
+ */
+/*
+ * Functions for testing and timing ptest, nextcand, prevcand.
+ *
+ * rlen(N) for N > 0 generates a random N-word positive integer.
+ *
+ * plen(N) for N > 0 generates an almost certainly prime positive
+ *	integer whose word-count is about N.
+ *
+ * clen(N) for N > 0 generates a composite odd N-word integer.
+ *
+ * ptimes(str, N [, n [, count [, skip, [, verbose]]]])
+ *	tests, and finds the runtime, for
+ *	ptest(x, count, skip) for n random almost certainly prime integers x
+ *	with word-count about N; n defaults to ceil(K1/abs(count)/(H1 + N^3)),
+ *	count to COUNT, skip to SKIP.
+ *
+ * ctimes(str, N [, n [, count [, skip, [, verbose]]]])
+ *	tests, and finds the runtime, for
+ *	ptest(x, count, skip) for n random composite integers x with word-count
+ *	about N; n defaults to ceil(K2/(H2 + N^3)), count to COUNT, skip
+ *	to SKIP.
+ *
+ * crtimes(str,a,b,n, [,count [, skip, [, verbose]]])
+ *	tests, and finds the runtime,
+ *	for ptest(x, count, skip) for n random integers x between a and b;
+ *	count defaults to COUNT, skip to SKIP.
+ *
+ * ntimes (str, N [,n, [, count [, skip [, residue [, modulus[,verb]]]]]]) tests
+ *	and finds the runtime for nextcand(...) and prevcand (...) for
+ *	n integers x with word-count about N, etc. n defaults to
+ *	ceil(K3/(H3 + N^3));
+ *
+ * testnextcand(str, N [, n [, count [, skip [, residue [, modulus [, verb]]]]])
+ *	performs tests of nextcand(x, count, skip, residue, modulus)
+ *	for n values of x with word-count N; n defaults to
+ *	ceil(K3/(H3 + N^3)), count to COUNT, skip to SKIP, residue to 0,
+ *	modulus to 1.
+ *
+ * testprevcand(str, N [, n [, count [, skip [, residue [, modulus [, verb]]]]])
+ *	performs tests of prevcand(x, count, skip, residue, modulus)
+ *	for n values of x with word-count N; n defaults to
+ *	ceil(K3/(H3 + N^3)), count to COUNT, skip to SKIP, residue to 0,
+ *	modulus to 1.
+ */
+
+global defaultverbose = 1;      /* default verbose value */
+global err;
+
+/*
+ * test defaults
+ */
+global BASEB = 32;
+global BASE = 2^BASEB;
+global COUNT = 5;
+global SKIP = 0;
+global RESIDUE = 0;
+global MODULUS = 1;
+
+/*
+ * internal test constants
+ */
+global K1 = 2^15;
+global H1 = 40;
+global K2 = 2^17;
+global H2 = 40;
+global K3 = 2^10;
+global H3 = 10;
+
+
+define rlen(N)
+{
+
+	if (!isint(N) || N <= 0)
+		quit "Bad argument for rlen";
+	return rand(BASE^(N-1), BASE^N);
+}
+
+define plen(N) = nextcand(rlen(N), 10, 0);
+
+define clen(N)
+{
+	local n, v;
+
+	do {
+		v = rlen(N);
+		if (iseven(v))
+			v++;
+	}
+	while
+		(ptest(v, 10, 0));
+	return v;
+}
+
+define ptimes(str, N, n, count, skip, verbose)
+{
+	local A, i, t, p, m;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	if (isnull(count))
+		count = COUNT;
+	if (isnull(n)) {
+		n = ceil(K1/abs(count)/(H1 + N^3));
+		if (verbose > 1) {
+			print "n =",n;
+		}
+	}
+	if (isnull(skip))
+		skip = SKIP;
+	mat A[n];
+	for (i = 0; i < n; i++)
+		A[i] = plen(N);
+	t = runtime();
+	for (i = 0; i < n; i++) {
+		p = ptest(A[i], count, skip);
+		if (!p) {
+			if (verbose > 0) {
+				printf("*** Error for x = %d\n", A[i]);
+				m++;
+			}
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+		} else {
+			t = round(runtime() - t, 4);
+			if (verbose > 1) {
+				printf("%d probable primes: time = %d\n", n, t);
+			} else {
+				printf("%d probable primes: passed\n", n);
+			}
+		}
+	}
+	return m;
+}
+
+define ctimes(str, N, n, count, skip, verbose)
+{
+	local A, i, r, t, p, m;
+
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	if (isnull(count))
+		count = COUNT;
+	if (isnull(n)) {
+		n = ceil(K2/(H2 + N^3));
+		if (verbose > 1) {
+			print "n =",n;
+		}
+	}
+	if (isnull(skip))
+		skip = SKIP;
+	mat A[n];
+	for (i = 0; i < n; i++)
+		A[i] = clen(N);
+	t = runtime();
+	for (i = 0; i < n; i++) {
+		p = ptest(A[i], count, skip);
+		if (p) {
+			if (verbose > 0) {
+				printf("*** Error, what should be rare has occurred for x = %d \n", A[i]);
+				m++;
+			}
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+		} else {
+			t = round(runtime() - t, 4);
+			if (verbose > 1) {
+				printf("%d probable primes: time = %d\n", n, t);
+			} else {
+				printf("%d probable primes: passed\n", n);
+			}
+		}
+	}
+	return m;
+}
+
+define crtimes(str, a, b, n, count, skip, verbose)
+{
+	local A, P, i, t, p, m;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	if (b < a)
+		swap(a,b);
+	b++;
+	if (isnull(count))
+		count = COUNT;
+	if (isnull(skip))
+		skip = SKIP;
+	mat A[n];
+	mat P[n];
+	for (i = 0; i < n; i++) {
+		A[i] = rand(a,b);
+		P[i] = ptest(A[i], 20, 0);
+	}
+	t = runtime();
+	for (i = 0; i < n; i++) {
+		p = ptest(A[i], count, skip);
+		if (p != P[i]) {
+			if (verbose > 0) {
+				printf("*** Apparent error for %s x = %d\n",
+					P[i] ? "prime" : "composite", A[i]);
+				++m;
+			}
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)?\n", m);
+		} else {
+			t = round(runtime() - t, 4);
+			if (verbose > 1) {
+				printf("%d probable primes: time = %d\n", n, t);
+			} else {
+				printf("%d probable primes: passed\n", n);
+			}
+		}
+	}
+	return m;
+}
+
+define ntimes(str, N, n, count, skip, residue, modulus, verbose)
+{
+	local A, i, t, p, tnext, tprev;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	if (isnull(count))
+		count = COUNT;
+	if (isnull(n)) {
+		n = ceil(K3/(H3 + N^3));
+		if (verbose > 1) {
+			print "n =",n;
+		}
+	}
+	if (isnull(skip))
+		skip = SKIP;
+	if (isnull(residue))
+		residue = RESIDUE;
+	if (isnull(modulus))
+		modulus = MODULUS;
+
+	mat A[n];
+	for (i = 0; i < n; i++)
+		A[i] = rlen(N);
+	t = runtime();
+	for (i = 0; i < n; i++) {
+		p = nextcand(A[i], count, skip, residue, modulus);
+	}
+	tnext = round(runtime() - t, 4);
+	t = runtime();
+	for (i = 0; i < n; i++) {
+		p = prevcand(A[i], count, skip, residue, modulus);
+	}
+	tprev = round(runtime() - t, 4);
+	if (verbose > 0) {
+		printf("%d evaluations, nextcand: %d, prevcand: %d\n", n, tnext, tprev);
+	}
+}
+
+define testnextcand(str, N, n, count, skip, residue, modulus, verbose)
+{
+	local p, x, y, i, m;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	if (isnull(count))
+		count = COUNT;
+	if (isnull(n)) {
+ 		n = ceil(K3/(H3 + N^3));
+		print "n =",n;
+	}
+	if (isnull(skip))
+		skip = SKIP;
+	if (isnull(residue))
+		residue = RESIDUE;
+	if (isnull(modulus))
+		modulus = MODULUS;
+	for (i = 0; i < n; i++) {
+		x = rlen(N);
+		y = nextcand(x, count, skip, residue, modulus);
+		p = testnext1(x, y, count, skip, residue, modulus);
+		if (p) {
+			m++;
+			if (verbose > 1) {
+				printf("*** Failure %d for x = %d\n", p, x);
+			}
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)?\n", m);
+		} else {
+			printf("%d successful tests\n", n);
+		}
+	}
+	return m;
+}
+
+define testnext1(x, y, count, skip, residue, modulus)
+{
+	if (y <= x)
+		return 1;
+	if (!ptest(y, count, skip))
+		return 2;
+	if (mne(y, residue, modulus))
+		return 3;
+	return 0;
+}
+
+define testprevcand(str, N, n, count, skip, residue, modulus, verbose)
+{
+	local p, x, y, i, m;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	if (isnull(count))
+		count = COUNT;
+	if (isnull(n)) {
+ 		n = ceil(K3/(H3 + N^3));
+		print "n =",n;
+	}
+	if (isnull(skip))
+		skip = SKIP;
+	if (isnull(residue))
+		residue = RESIDUE;
+	if (isnull(modulus))
+		modulus = MODULUS;
+	for (i = 0; i < n; i++) {
+		x = rlen(N);
+		y = prevcand(x, count, skip, residue, modulus);
+		p = testprev1(x, y, count, skip, residue, modulus);
+		if (p) {
+			m++;
+			if (verbose > 1) {
+				printf("*** Failure %d for x = %d\n", p, x);
+			}
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)?\n", m);
+		} else {
+			printf("%d successful tests\n", n);
+		}
+	}
+	return m;
+}
+
+
+define testprev1(x, y, count, skip, residue, modulus)
+{
+	if (y >= x)
+		return 1;
+	if (!ptest(y, count, skip))
+		return 2;
+	if (mne(y, residue, modulus))
+		return 3;
+	return 0;
+}
+
+/*
+ * test4000 - perform all of the above tests a bunch of times
+ */
+define test4000(v, tnum)
+{
+	local n;	/* test parameter */
+
+	/*
+	 * set test parameters
+	 */
+	srand(4000e4000);
+
+	/*
+	 * test a lot of stuff
+	 */
+	err += ptimes(strcat(str(tnum++),": ptimes(1,250)"), 1, 250,,,v);
+	err += ptimes(strcat(str(tnum++),": ptimes(3,50)"), 3, 50,,,v);
+	err += ptimes(strcat(str(tnum++),": ptimes(5,20)"), 5, 20,,,v);
+
+	err += ctimes(strcat(str(tnum++),": ctimes(1,7500)"), 1, 7500,,,v);
+	err += ctimes(strcat(str(tnum++),": ctimes(3,500)"), 3, 500,,,v);
+	err += ctimes(strcat(str(tnum++),": ctimes(5,200)"), 5, 200,,,v);
+
+	err += crtimes(strcat(str(tnum++),": crtimes(2^30,2^31,2500)"),
+		2^30, 2^31, 2500,,,v);
+	err += crtimes(strcat(str(tnum++),": crtimes(2^300,2^301,75)"),
+		2^300, 2^301, 75,,,v);
+
+	err += testprevcand(strcat(str(tnum++),": testprevcand(1,250)"),
+		1, 250, ,,,,v);
+	err += testprevcand(strcat(str(tnum++),": testprevcand(3,25)"),
+		3, 25, ,,,,v);
+	err += testprevcand(strcat(str(tnum++),": testprevcand(5,10)"),
+		5, 10, ,,,,v);
+
+	err += testnextcand(strcat(str(tnum++),": testnextcand(1,250)"),
+		1, 250, ,,,,v);
+	err += testnextcand(strcat(str(tnum++),": testnextcand(3,25)"),
+		3, 25, ,,,,v);
+	err += testnextcand(strcat(str(tnum++),": testnextcand(5,10)"),
+		5, 10, ,,,,v);
+
+	/*
+	 * report results
+	 */
+	if (v > 1) {
+		if (err) {
+			print "***", err, "error(s) found in testall";
+		} else {
+			print "no errors in testall";
+		}
+	}
+	return tnum;
+}
+
+
+global lib_debug;
+
+if (lib_debug >= 0) {
+	print "global defaultverbose";
+	print "global err";
+	print "global BASEB";
+	print "global BASE";
+	print "global COUNT";
+	print "global SKIP";
+	print "global RESIDUE";
+	print "global MODULUS";
+	print "global K1";
+	print "global H1";
+	print "global K2";
+	print "global H2";
+	print "global K3";
+	print "global H3";
+	print "plen(N) defined";
+	print "clen(N) defined";
+	print "ptimes(str, N, n, count, skip, verbose) defined";
+	print "ctimes(str, N, n, count, skip, verbose) defined";
+	print "crtimes(str, a, b, n, count, skip, verbose) defined";
+	print "ntimes(str, N, n, count, skip, residue, mod, verbose) defined";
+	print "testnextcand(str, N, n, cnt, skip, res, mod, verbose) defined";
+	print "testnext1(x, y, count, skip, residue, modulus) defined";;
+	print "testprevcand(str, N, n, cnt, skip, res, mod, verbose) defined";
+	print "testprev1(x, y, count, skip, residue, modulus) defined";
+	print "test4000(verbose, tnum) defined";
+}

lib/test4100.cal

@@ -0,0 +1,493 @@
+/*
+ * Copyright (c) 1996 Ernest Bowen and Landon Curt Noll
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * By: Ernest Bowen and Landon Curt Noll
+ *     ernie@neumann.une.edu.au and chongo@toad.com
+ *
+ * This library is used by the 4100 series of the regress.cal test suite.
+ */
+/*
+ * Some severe tests and timing functions for REDC functions and pmod.
+ *
+ * testall(str,n,N,M,verbose)
+ *	performs n tests using arguments x, y, ...
+ *	randomly selected from [-N, N) or when nonnegative values are
+ *	required, [0, N), and m an odd positive integer in [1,N],
+ *	and where a "small" (say less than 10000) exponent k is to be
+ *	used (when computing x^k % m directly) k is random in [0,M).
+ *	Default values for N and M are 1e20 and 100.
+ *
+ * times(str,N,n,verbose)
+ *	gives times for n evaluations of rcin, etc. with
+ *	N-word arguments; default n is ceil(K1/power(N,1.585).
+ *
+ * powtimes(str, N1,N2,n, verbose)
+ *	gives times for n evaluations of pmod(x,k,m)
+ *	for the three algorithms "small", "normal", "bignum" that
+ *	pmod may use, and equivalent functions rcpow(xr,k,m) for
+ *	"small" or "bignum" cases, where xr = rcin(x,m).  The
+ *	modulus m is a random positive odd N1-word number; x has
+ *	random integer values in [0, m-1]; k has random N2-word values.
+ *	N2 defaults to 1; n defaults to ceil(K2/power(N1,1.585)/N2).
+ *
+ * inittimes(str, N, n, verbose)
+ *	displays the times and tests n evaluations of rcin(x,m)
+ *	and rcout(x,m) where m is a random positive odd N-word integer,
+ *	x is a random integer in [0, m-1]; n defaults to ceil(K1/N^2).
+ *
+ * rlen_4100(N)
+ *	generates a random positive N-word integer.  The global
+ *	BASEB should be set to the word-size for the computer being
+ *	used.  The parameters K1, K2 which control the default n
+ *	should be adjusted to give reasonable runtimes.
+ *
+ * olen(N)
+ *	generates a random odd positive N-word number.
+ *
+ */
+
+global defaultverbose = 1;      /* default verbose value */
+global err;
+
+/*
+ * test defaults
+ */
+global K1 = 2^17;
+global K2 = 2^12;
+global BASEB = 16;
+global BASE = 2^BASEB;
+
+define rlen_4100(N) = rand(BASE^(N-1), BASE^N);
+
+define olen(N)
+{
+	local v;
+
+	v = rlen_4100(N);
+	if (iseven(v))
+		v++;
+	return v;
+}
+
+define test4101(x,y,m,k,z1,z2,verbose)
+{
+	local xr, yr, v, w, oneone;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	xr = rcin(x,m);
+	yr = rcin(y,m);
+	oneone = rcin(rcin(1,m),m);
+
+	if (xr >= m || xr < 0) {
+		if (verbose > 1)
+			printf("Failure 1 for x = %d, m = %d\n", x, m);
+		return 1;
+	}
+	if (rcin(x * y, m) != mod(xr * y, m, 0)) {
+		if (verbose > 1) {
+			printf("Failure 2 for x = %d, y = %d, m = %d\n",
+				x, y, m);
+		}
+		return 2;
+	}
+	if (rcout(xr, m) != x % m) {
+		if (verbose > 1)
+			printf("Failure 3 for x = %d, m = %d\n", x, m);
+		return 3;
+	}
+	if (rcout(rcmul(xr,yr,m),m) != mod(x * y, m, 0)) {
+		if (verbose > 1)
+			printf("Failure 4 for x = %d, y = %d, m = %d\n",
+				 x, y, m);
+		return 4;
+	}
+	if (rcmul(x,yr,m) != mod(x * y, m, 0)) {
+		if (verbose > 1)
+			printf("Failure 5 for x = %d, y = %d, m = %d\n",
+				 x, y, m);
+		return 5;
+	}
+	if (rcin(rcmul(x,y,m),m) != mod(x * y, m, 0)) {
+		if (verbose > 1)
+			printf("Failure 6 for x = %d, y = %d, m = %d\n",
+				x, y, m);
+		return 6;
+	}
+	if (rcout(rcsq(xr,m),m) != mod(x^2, m, 0)) {
+		if (verbose > 1)
+			printf("Failure 7 for x = %d, m = %d\n", x, m);
+		return 7;
+	}
+	if (rcin(rcsq(x,m),m) != mod(x^2, m, 0)) {
+		if (verbose > 1)
+			printf("Failure 8 for x = %d, m = %d\n",
+				x, y, m);
+		return 8;
+	}
+	if (rcout(rcpow(xr,k,m),m) != mod(x^k, m, 0)) {
+		if (verbose > 1)
+			printf("Failure 9 for x = %d, m = %d, k = %d\n",
+				x, m, k);
+		return 9;
+	}
+	if (rcpow(x,k,m) != rcin(rcout(x,m)^k, m)) {
+		if (verbose > 1)
+			printf("Failure 10: x = %d, k = %d, m = %d\n",
+				x, k, m);
+		return 10;
+	}
+	if (rcpow(x, z1 * z2, m) != rcpow(rcpow(x,z1,m), z2, m)) {
+		if (verbose > 1)
+			printf("Failure 11: x = %d, z1 = %d, z2 = %d, m = %d\n",
+			x, z1, z2, m);
+		return 11;
+	}
+	if (xr != rcmul(oneone, x, m)) {
+		if (verbose > 1)
+			printf("Failure 12: x = %d, m = %d\n", x, m);
+		return 12;
+	}
+
+	return 0;
+}
+
+define testall(str,n,N,M,verbose)
+{
+	local i, p, x, y, z1, z2, c, k, m;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	if (isnull(N))
+		N = 1e20;
+	if (isnull(M))
+		M = 100;
+	c = 0;
+	for (i = 0; i < n; i++) {
+		x = rand(-N, N);
+		y = rand(-N, N);
+		z1 = rand(N);
+		z2 = rand(N);
+		c = rand(N);
+		if (iseven(c))
+			c++;
+		k = rand(M);
+		if (verbose > 1)
+			printf("x = %d, y = %d, c = %d, k = %d\n", x, y, c, k);
+		p = test4101(x,y,c,k,z1,z2);
+		if (p) {
+			m++;
+			if (verbose > 0) {
+				printf("*** Failure %d in test %d\n", p, i);
+			}
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+		} else {
+			printf("passed %d tests\n", n);
+		}
+	}
+	return m;
+}
+
+define times(str,N,n,verbose)
+{
+	local m, m2, A, B, C, x, y, t, i, z;
+	local trcin, trcout, trcmul, trcsq;
+	local tmul, tsq, tmod, tquomod;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = olen(N);
+	m2 = m^2;
+	if (isnull(n)) {
+		n = ceil(K1/power(N,1.585));
+		if (verbose > 1)
+			printf("n = %d\n", n);
+	}
+	mat A[n];
+	mat B[n];
+	mat C[n];
+	for (i = 0; i < n; i++) {
+		A[i] = rand(m);
+		B[i] = rand(m);
+		C[i] = rand(m2);
+	}
+	z = rcin(0,m);	/* to initialize redc and maybe lastmod information */
+	t = runtime();
+	for (i = 0; i < n; i++)
+		z = rcin(A[i],m);
+	trcin = round(runtime() - t, 3);
+	t = runtime();
+	for (i = 0; i < n; i++)
+		z = rcout(A[i],m);
+	trcout = round(runtime() - t, 3);
+	t = runtime();
+	for (i = 0; i < n; i++)
+		z = rcmul(A[i],B[i],m);
+	trcmul = round(runtime() - t, 3);
+	t = runtime();
+	for (i = 0; i < n; i++)
+		z = rcsq(A[i],m);
+	trcsq = round(runtime() - t, 3);
+	t = runtime();
+	for (i = 0; i < n; i++)
+		z = A[i] * B[i];
+	tmul = round(runtime() - t, 3);
+	t = runtime();
+	for (i = 0; i < n; i++)
+		z = A[i]^2;
+	tsq = round(runtime() - t, 3);
+	t = runtime();
+	for (i = 0; i < n; i++)
+		z = C[i] % A[i];
+	tmod = round(runtime() - t, 3);
+	t = runtime();
+	for (i = 0; i < n; i++)
+		quomod(C[i], A[i], x, y);
+	tquomod = round(runtime() - t,3);
+
+	if (verbose > 1) {
+		printf("rcin: %d, rcout: %d, rcmul: %d, rcsq: %d\n",
+			trcin, trcout, trcmul, trcsq);
+		printf("%s: mul: %d, sq: %d, mod: %d, quomod: %d\n",
+			str, tmul, tsq, tmod, tquomod);
+	} else if (verbose > 0) {
+		printf("no error(s)\n");
+	}
+	return 0;
+}
+
+
+define powtimes(str, N1, N2, n, verbose)
+{
+	local A, Ar, B, v, i, t, z1, z2, z3, z4, z5, cp, crc;
+	local tsmall, tnormal, tbignum, trcsmall, trcbig, m;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+
+	if (isnull(N2))
+		N2 = 1;
+
+	if (isnull(n)) {
+		n = ceil(K2/power(N1, 1.585)/N2);
+		printf ("n  = %d\n", n);
+	}
+	mat A[n];
+	mat Ar[n];
+	mat B[n];
+	v = olen(N1);
+
+	cp = config("pow2", 1);
+	crc = config("redc2", 1);
+
+	/* initialize redc and lastmod info */
+
+	z1 = z2 = z3 = z4 = z5 = rcin(0,v);
+
+	for (i = 0; i < n; i++) {
+		A[i] = rand(v);
+		Ar[i] = rcin(A[i], v);
+		B[i] = rlen_4100(N2);
+	}
+	t = runtime();
+	for (i = 0; i < n; i++)
+		z1 += pmod(A[i], B[i], v);
+	tbignum = round(runtime() - t, 4);
+	config("pow2", 1e6);
+	t = runtime();
+	for (i = 0; i < n; i++)
+		z2 += pmod(A[i], B[i], v);
+	tnormal = round(runtime() - t, 4);
+	config("redc2",1e6);
+	t = runtime();
+	for (i = 0; i < n; i++)
+		z3 += pmod(A[i], B[i], v);
+	tsmall = round(runtime() - t, 4);
+	t = runtime();
+	for (i = 0; i < n; i++)
+		z4 += rcpow(Ar[i], B[i], v);
+	trcsmall = round(runtime() - t, 4);
+	config("redc2", 1);
+	t = runtime();
+	for (i = 0; i < n; i++)
+		z5 += rcpow(Ar[i], B[i], v);
+	trcbig = round(runtime() - t, 4);
+
+	if (z1 != z2) {
+		++m;
+		if (verbose > 0) {
+			print "*** z1 != z2";
+		}
+	} else if (z1 != z3) {
+		++m;
+		if (verbose > 0) {
+			print "*** z1 != z3";
+		}
+	} else if (z2 != z3) {
+		++m;
+		if (verbose > 0) {
+			print "*** z2 != z3";
+		}
+	} else if (rcout(z4, v) != z1 % v) {
+		++m;
+		if (verbose > 0) {
+			print "*** z4 != z1";
+		}
+	} else if (z4 != z5) {
+		++m;
+		if (verbose > 0) {
+			print "*** z4 != z5";
+		}
+	}
+	config("pow2", cp);
+	config("redc2", crc);
+	if (verbose > 1) {
+	}
+	if (verbose > 1) {
+		printf("Small: %d, normal: %d, bignum: %d\n",
+			tsmall, tnormal, tbignum);
+		printf("%s: rcsmall: %d, rcbig: %d\n",
+			str, trcsmall, trcbig);
+	} else if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)\n", m);
+		} else {
+			printf("passed\n");
+		}
+	}
+	return 0;
+}
+
+define inittimes(str,N,n,verbose)
+{
+	local A, M, B, R, i, t, trcin, trcout, m;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	m = 0;
+	if (isnull(n)) {
+		n = ceil(K1/N^2);
+		if (verbose > 1) {
+			printf ("n = %d\n", n);
+		}
+	}
+	mat A[n];
+	mat M[n];
+	mat B[n];
+	mat R[n];
+	for (i = 0; i < n; i++) {
+		M[i] = olen(N);
+		A[i] = rand(M[i]);
+	}
+	t = runtime();
+	for (i = 0; i < n; i++)
+		R[i] = rcin(A[i], M[i]);
+	trcin = round(runtime() - t, 4);
+	for (i = 0; i < n; i++)
+		B[i] = rcout(R[i], M[i]);
+	trcout = round(runtime() - t, 4);
+	for (i = 0; i < n; i++) {
+		if (B[i] != A[i]) {
+			++m;
+			if (verbose > 0) {
+				print "*** Error!!!";
+			}
+			break;
+		}
+	}
+	if (verbose > 0) {
+		if (m) {
+			printf("*** %d error(s)?\n", m);
+		} else {
+			if (verbose > 1) {
+			    printf("%d successful tests: rcin: %d, rcout: %d\n",
+			      n, trcin, trcout);
+			} else {
+			    printf("%d successful tests: passed\n", n);
+			}
+		}
+	}
+	return m;
+}
+
+/*
+ * test4100 - perform all of the above tests a bunch of times
+ */
+define test4100(v, tnum)
+{
+	local n;	/* test parameter */
+
+	/*
+	 * set test parameters
+	 */
+	srand(4100e4100);
+
+	/*
+	 * test a lot of stuff
+	 */
+	err += testall(strcat(str(tnum++),": testall(10,,500)"), 10,, 500, v);
+	err += testall(strcat(str(tnum++),": testall(200)"), 200,,, v);
+
+	err += times(strcat(str(tnum++),": times(3,3000)"), 3, 3000, v);
+	err += times(strcat(str(tnum++),": times(30,300)"), 30, 300, v);
+	err += times(strcat(str(tnum++),": times(300,30)"), 300, 30, v);
+	err += times(strcat(str(tnum++),": times(1000,3)"), 1000, 3, v);
+
+	err += powtimes(strcat(str(tnum++),": powtimes(100)"),100,,v);
+	err += powtimes(strcat(str(tnum++),": powtimes(250)"),250,,v);
+
+	err += inittimes(strcat(str(tnum++),": inittimes(10)"),10,,v);
+	err += inittimes(strcat(str(tnum++),": inittimes(100,70)"),100,70,v);
+	err += inittimes(strcat(str(tnum++),": inittimes(1000,4)"),1000,4,v);
+
+	/*
+	 * report results
+	 */
+	if (v > 1) {
+		if (err) {
+			print "***", err, "error(s) found in testall";
+		} else {
+			print "no errors in testall";
+		}
+	}
+	return tnum;
+}
+
+global lib_debug;
+
+if (lib_debug >= 0) {
+	print "global defaultverbose";
+	print "global err";
+	print "global K1";
+	print "global K2";
+	print "global BASEB";
+	print "global BASE";
+	print "rlen_4100(N) defined";
+	print "olen(N) defined";
+	print "test4101(x, y, m, k, z1, z2) defined";
+	print "testall(str, n, N, M, verbose) defined";
+	print "times(str, N, n, verbose) defined";
+	print "powtimes(str, N1, N2, n, verbose) defined";
+	print "inittimes(str, N, n, verbose) defined";
+	print "test4100(verbose, tnum) defined";
+}

lib/test4600.cal

@@ -0,0 +1,311 @@
+/*
+ * Copyright (c) 1996 Ernest Bowen and Landon Curt Noll
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * By: Ernest Bowen and Landon Curt Noll
+ *     ernie@neumann.une.edu.au and chongo@toad.com
+ *
+ * This library is used by the 4600 series of the regress.cal test suite.
+ */
+
+
+global defaultverbose = 1		/* default verbose value */
+global err;
+
+/*
+ * test globals
+ */
+global A, f, pos;
+
+define stest(str, verbose)
+{
+	local x;
+
+	/* setup */
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	x = rm("junk4600");
+
+	/*
+	 * do file operations
+	 */
+	f = fopen("junk4600", "w");
+	if (iserror(f)) {
+		print 'failed';
+		print '**** fopen("junk4600", "w") failed';
+		return 1;
+	}
+	if (iserror(fputs(f,
+	 "Fourscore and seven years ago our fathers brought forth\n",
+	 "on this continent a new nation, conceived in liberty and dedicated\n",
+	 "to the proposition that all men are created equal.\n"))) {
+		print 'failed';
+		print '**** fputs(f, "Fourscore ... failed';
+		return 1;
+	}
+	if (iserror(freopen(f, "r"))) {
+		print 'failed';
+		print '**** iserror(freopen(f, "r")) failed';
+		return 1;
+	}
+	if (iserror(rewind(f))) {
+		print 'failed';
+		print '**** iserror(rewind(f)) failed';
+		return 1;
+	}
+	if (search(f, "and") != 10) {
+		print 'failed';
+		print '**** search(f, "and") != 10 failed';
+		return 1;
+	}
+	if (ftell(f) != 13) {
+		print 'failed';
+		print '**** ftell(f) != 13 failed';
+		return 1;
+	}
+	if (search(f, "and") != 109) {
+		print 'failed';
+		print '**** search(f, "and") != 109 failed';
+		return 1;
+	}
+	if (ftell(f) != 112) {
+		print 'failed';
+		print '**** ftell(f) != 112 failed';
+		return 1;
+	}
+	if (!isnull(search(f, "and"))) {
+		print 'failed';
+		print '**** !isnull(search(f, "and")) failed';
+		return 1;
+	}
+	if (ftell(f) != 172) {
+		print 'failed';
+		print '**** ftell(f) != 172 failed';
+		return 1;
+	}
+	if (rsearch(f, "and") != 109) {
+		print 'failed';
+		print '**** rsearch(f, "and") != 109 failed';
+		return 1;
+	}
+	if (ftell(f) != 112) {
+		print 'failed';
+		print '**** ftell(f) != 112 failed';
+		return 1;
+	}
+	if (iserror(fseek(f, -1, 1))) {
+		print 'failed';
+		print '**** iserror(fseek(f, -1, 1)) failed';
+		return 1;
+	}
+	if (rsearch(f, "and") != 10) {
+		print 'failed';
+		print '**** rsearch(f, "and") != 10 failed';
+		return 1;
+	}
+	if (ftell(f) != 13) {
+		print 'failed';
+		print '**** ftell(f) != 13 failed';
+		return 1;
+	}
+	if (iserror(fseek(f, -1, 1))) {
+		print 'failed';
+		print '**** iserror(fseek(f, -1, 1)) failed';
+		return 1;
+	}
+	if (!isnull(rsearch(f, "and"))) {
+		print 'failed';
+		print '**** !isnull(rsearch(f, "and")) failed';
+		return 1;
+	}
+	if (ftell(f) != 0) {
+		print 'failed';
+		print '**** ftell(f) != 0 failed';
+		return 1;
+	}
+	if (iserror(fclose(f))) {
+		print 'failed';
+		print '**** iserror(fclose(f)) failed';
+		return 1;
+	}
+
+	/*
+	 * cleanup
+	 */
+	x = rm("junk4600");
+	if (verbose > 0) {
+		printf("passed\n");
+	}
+	return 0;
+}
+
+define ttest(str, m, n, verbose)
+{
+	local a, s, i, j;
+
+	if (isnull(verbose))
+		verbose = defaultverbose;
+	if (verbose > 0) {
+		print str:":",:;
+	}
+	i = rm("junk4600");
+	f = fopen("junk4600", "w");
+
+	if (isnull(n))
+		n = 4;
+	if (isnull(m))
+		m = 4;
+
+	mat A[m];
+	mat pos[m + 1];
+
+	pos[0] = 0;
+	for (i = 0; i < m; i++) {
+		j = 1 + randbit(n);
+		a = "";
+		while (j-- > 0)
+			a = strcat(a, char(rand(1, 256)));
+		A[i] = a;
+		fputs(f, a);
+		pos[i+1] = ftell(f);
+		if (verbose > 1)
+			printf("A[%d] has length %d\n", i, strlen(a));
+	}
+	if (verbose > 1)
+		printf("File has size %d\n", pos[i]);
+	freopen(f, "r");
+	if (size(f) != pos[i]) {
+		print 'failed';
+		printf("**** Failure 1 for file size\n");
+		return 1;
+	}
+	for (i = 0; i < m; i++) {
+		rewind(f);
+		for (;;) {
+			j = search(f, A[i]);
+			if (isnull(j) || j > pos[i]) {
+				print 'failed';
+				printf("**** Failure 2 for i = %d\n", i);
+				return 1;
+			}
+			if (j == pos[i])
+				break;
+			fseek(f, j + 1, 0);
+
+		}
+		if (ftell(f) != pos[i + 1]) {
+			print 'failed';
+			printf("**** Failure 3 for i = %d\n", i);
+			return 1;
+		}
+	}
+	for (i = m - 1; i >= 0; i--) {
+		fseek(f, 0, 2);
+		for (;;) {
+			j = rsearch(f, A[i]);
+			if (isnull(j) || j < pos[i]) {
+				print 'failed';
+				printf("**** Failure 4 for i = %d\n", i);
+				return 1;
+			}
+			if (j == pos[i])
+				break;
+			fseek(f, -1, 1);
+		}
+		if (ftell(f) != pos[i + 1]) {
+			print 'failed';
+			printf("**** Failure 5 for i = %d\n", i);
+			return 1;
+		}
+	}
+	i = rm("junk4600");
+	if (verbose > 0) {
+		printf("passed\n");
+	}
+	return 0;
+}
+
+define sprint(x)
+{
+	local i, n;
+
+	n = strlen(x);
+	for (i = 1; i <= n; i++) print ord(substr(x, i, 1)),;
+	print;
+}
+
+define findline(f,s)
+{
+
+	if (!isfile(f))
+		quit "First argument to be a file";
+	if (!isstr(s))
+		quit "Second argument to be a string";
+	if (!isnull(search(f,s))) {
+		rsearch(f, "\n");
+		print fgetline(f);
+	}
+}
+
+define findlineold(f,s)
+{
+	local str;
+
+	if (!isfile(f))
+		quit "First argument to be a file";
+	if (!isstr(s))
+		quit "Second argument to be a string";
+
+	while (!isnull(str = fgetline(f)) && strpos(str, s) == 0);
+		print str;
+}
+
+/*
+ * test4600 - perform all of the above tests a bunch of times
+ */
+define test4600(v, tnum)
+{
+	local n;	/* test parameter */
+	local i;
+
+	/*
+	 * set test parameters
+	 */
+	srand(4600e4600);
+
+	/*
+	 * test a lot of stuff
+	 */
+	for (i=0; i < 10; ++i) {
+		err += ttest(strcat(str(tnum++),
+			     ": ttest(",str(i),",",str(i),")"), i, i, v);
+		err += stest(strcat(str(tnum++), ": stest()"), v);
+	}
+
+	/*
+	 * report results
+	 */
+	if (v > 1) {
+		if (err) {
+			print "****", err, "error(s) found in testall";
+		} else {
+			print "no errors in testall";
+		}
+	}
+	return tnum;
+}
+
+global lib_debug;
+
+if (lib_debug >= 0) {
+	print "stest(str [, verbose]) defined";
+	print "ttest([m, [n [,verbose]]]) defined";
+	print "sprint(x) defined";
+	print "findline(f,s) defined";
+	print "findlineold(f,s) defined";
+	print "test4600(verbose, tnum) defined";
+}

lib/unitfrac.cal

@@ -0,0 +1,35 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Represent a fraction as sum of distinct unit fractions.
+ * The output is the unit fractions themselves, and in square brackets,
+ * the number of digits in the numerator and denominator of the value left
+ * to be found.  Numbers larger than 3.5 become very difficult to calculate.
+ */
+
+define unitfrac(x)
+{
+	local	d, di, n;
+
+	if (x <= 0)
+		quit "Non-positive argument";
+	d = 2;
+	do {
+		n = int(1 / x) + 1;
+		if (n > d)
+			d = n;
+		di = 1/d;
+		print '  [': digits(num(x)): '/': digits(den(x)): ']',, di;
+		x -= di;
+		d++;
+	} while ((num(x) > 1) || (x == di) || (x == 1));
+	print '  [1/1]',, x;
+}
+
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "unitfrac(x) defined";
+}

lib/varargs.cal

@@ -0,0 +1,29 @@
+/*
+ * Copyright (c) 1995 David I. Bell
+ * Permission is granted to use, distribute, or modify this source,
+ * provided that this copyright notice remains intact.
+ *
+ * Example program to use 'varargs'.
+ *
+ * Program to sum the cubes of all the specified numbers.
+ */
+
+define sc()
+{
+	local s, i;
+
+	s = 0;
+	for (i = 1; i <= param(0); i++) {
+		if (!isnum(param(i))) {
+			print "parameter",i,"is not a number";
+			continue;
+		}
+		s += param(i)^3;
+	}
+	return s;
+}
+
+global lib_debug;
+if (lib_debug >= 0) {
+    print "sc(a, b, ...) defined";
+}