home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Usenet 1994 October
/
usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso
/
unix
/
volume27
/
calc-2.9.0
/
part18
< prev
next >
Wrap
Text File
|
1993-12-07
|
61KB
|
2,060 lines
Newsgroups: comp.sources.unix
From: dbell@canb.auug.org.au (David I. Bell)
Subject: v27i145: calc-2.9.0 - arbitrary precision C-like programmable calculator, Part18/19
References: <1.755316719.21314@gw.home.vix.com>
Sender: unix-sources-moderator@gw.home.vix.com
Approved: vixie@gw.home.vix.com
Submitted-By: dbell@canb.auug.org.au (David I. Bell)
Posting-Number: Volume 27, Issue 145
Archive-Name: calc-2.9.0/part18
#!/bin/sh
# this is part 18 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc2.9.0/lib/lucas.cal continued
#
CurArch=18
if test ! -r s2_seq_.tmp
then echo "Please unpack part 1 first!"
exit 1; fi
( read Scheck
if test "$Scheck" != $CurArch
then echo "Please unpack part $Scheck next!"
exit 1;
else exit 0; fi
) < s2_seq_.tmp || exit 1
echo "x - Continuing file calc2.9.0/lib/lucas.cal"
sed 's/^X//' << 'SHAR_EOF' >> calc2.9.0/lib/lucas.cal
X quit "bad args: testval must be an integer";
X }
X if (!isint(v1)) {
X quit "bad args: v1 must be an integer";
X }
X if (testval <= 0) {
X quit "bogus arg: testval is <= 0";
X }
X if (v1 <= 0) {
X quit "bogus arg: v1 is <= 0";
X }
X
X /*
X * enforce the h mod rules
X */
X if (h%2 == 0) {
X quit "h must not be even";
X }
X
X /*
X * enforce the h > 0 and n >= 2 rules
X */
X if (h <= 0 || n < 1) {
X quit "reduced args violate the rule: 0 < h < 2^n";
X }
X hbits = highbit(h);
X if (hbits >= n) {
X quit "reduced args violate the rule: 0 < h < 2^n";
X }
X
X /*
X * build up u2 based on the reversed bits of h
X */
X /* setup for bit loop */
X r = v1;
X s = (r^2 - 2);
X
X /*
X * deal with small h as a special case
X *
X * The h value is odd > 0, and it needs to be
X * at least 2 bits long for the loop below to work.
X */
X if (h == 1) {
X ldebug("gen_u0", "quick h == 1 case");
X return r%testval;
X }
X
X /* cycle from second highest bit to second lowest bit of h */
X for (i=hbits-1; i > 0; --i) {
X
X /* bit(i) is 1 */
X if (isset(h,i)) {
X
X /* compute v(2n+1) = v(r+1)*v(r)-v1 */
X r = (r*s - v1) % testval;
X
X /* compute v(2n+2) = v(r+1)^2-2 */
X s = (s^2 - 2) % testval;
X
X /* bit(i) is 0 */
X } else {
X
X /* compute v(2n+1) = v(r+1)*v(r)-v1 */
X s = (r*s - v1) % testval;
X
X /* compute v(2n) = v(r)^-2 */
X r = (r^2 - 2) % testval;
X }
X }
X
X /* we know that h is odd, so the final bit(0) is 1 */
X r = (r*s - v1) % testval;
X
X /* compute the final u2 return value */
X return r;
X}
X
X/*
X * Trial tables used by gen_v1()
X *
X * When h mod 3 == 0, one needs particular values of D, a and b (see gen_v1
X * documentation) in order to find a value of v(1).
X *
X * This table defines 'quickmax' possible tests to be taken in ascending
X * order. The v1_qval[x] refers to a v(1) value from Ref1, Table 1. A
X * related D value is found in d_qval[x]. All D values expect d_qval[1]
X * are also taken from Ref1, Table 1. The case of D == 21 as listed in
X * Ref1, Table 1 can be changed to D == 7 for the sake of the test because
X * of {note 6}.
X *
X * It should be noted that the D values all satisfy the selection values
X * as outlined in the gen_v1() function comments. That is:
X *
X * D == P*(2^f)*(3^g)
X *
X * where f == 0 and g == 0, P == D. So we simply need to check that
X * one of the following two cases are true:
X *
X * P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1
X * P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1
X *
X * In all cases, the value of r is:
X *
X * r == Q*(2^j)*(3^k)*(z^2)
X *
X * where Q == 1. No further processing is needed to compute v(1) when r
X * is of this form.
X */
Xquickmax = 8;
Xmat d_qval[quickmax];
Xmat v1_qval[quickmax];
Xd_qval[0] = 5; v1_qval[0] = 3; /* a=1 b=1 r=4 */
Xd_qval[1] = 7; v1_qval[1] = 5; /* a=3 b=1 r=12 D=21 */
Xd_qval[2] = 13; v1_qval[2] = 11; /* a=3 b=1 r=4 */
Xd_qval[3] = 11; v1_qval[3] = 20; /* a=3 b=1 r=2 */
Xd_qval[4] = 29; v1_qval[4] = 27; /* a=5 b=1 r=4 */
Xd_qval[5] = 53; v1_qval[5] = 51; /* a=53 b=1 r=4 */
Xd_qval[6] = 17; v1_qval[6] = 66; /* a=17 b=1 r=1 */
Xd_qval[7] = 19; v1_qval[7] = 74; /* a=38 b=1 r=2 */
X
X/*
X * gen_v1 - compute the v(1) for a given h*2^n-1 if we can
X *
X * This function assumes:
X *
X * n > 2 (n==2 has already been eliminated)
X * h mod 2 == 1
X * h < 2^n
X * h*2^n-1 mod 3 != 0 (h*2^n-1 has no small factors, such as 3)
X *
X * The generation of v(1) depends on the value of h. There are two cases
X * to consider, h mod 3 != 0, and h mod 3 == 0.
X *
X ***
X *
X * Case 1: (h mod 3 != 0)
X *
X * This case is easy and always finds v(1).
X *
X * In Ref1, page 869, one finds that if: (or see Ref2, page 131-132)
X *
X * h mod 6 == +/-1
X * h*2^n-1 mod 3 != 0
X *
X * which translates, gives the functions assumptions, into the condition:
X *
X * h mod 3 != 0
X *
X * If this case condition is true, then:
X *
X * u(0) = (2+sqrt(3))^h + (2-sqrt(3))^h (see Ref1, page 869)
X * = (2+sqrt(3))^h + (2+sqrt(3))^(-h)
X *
X * and since Ref1, Theorem 5 states:
X *
X * u(0) = alpha^h + alpha^(-h)
X * r = abs(2^2 - 1^2*3) = 1
X *
X * and the bottom of Ref1, page 872 states:
X *
X * v(x) = alpha^x + alpha^(-x)
X *
X * If we let:
X *
X * alpha = (2+sqrt(3))
X *
X * then
X *
X * u(0) = v(h)
X *
X * so we simply return
X *
X * v(1) = alpha^1 + alpha^(-1)
X * = (2+sqrt(3)) + (2-sqrt(3))
X * = 4
X *
X ***
X *
X * Case 2: (h mod 3 == 0)
X *
X * This case is not so easy and finds v(1) in most all cases. In this
X * version of this program, we will simply return -1 (failure) if we
X * hit one of the cases that fall thru the cracks. This does not happen
X * often, so this is not too bad.
X *
X * Ref1, Theorem 5 contains the following definitions:
X *
X * r = abs(a^2 - b^2*D)
X * alpha = (a + b*sqrt(D))^2/r
X *
X * where D is 'square free', and 'alpha = epsilon^s' (for some s>0) are units
X * in the quadratic field K(sqrt(D)).
X *
X * One can find possible values for a, b and D in Ref1, Table 1 (page 872).
X * (see the file lucas_tbl.cal)
X *
X * Now Ref1, Theorem 5 states that if:
X *
X * L(D, h*2^n-1) = -1 [condition 1]
X * L(r, h*2^n-1) * (a^2 - b^2*D)/r = -1 [condition 2]
X *
X * where L(x,y) is the Legendre symbol (see below), then:
X *
X * u(0) = alpha^h + alpha^(-h)
X *
X * The bottom of Ref1, page 872 states:
X *
X * v(x) = alpha^x + alpha^(-x)
X *
X * thus since:
X *
X * u(0) = v(h)
X *
X * so we want to return:
X *
X * v(1) = alpha^1 + alpha^(-1)
X *
X * Therefore we need to take a given (D,a,b), determine if the two conditions
X * are true, and return the related v(1).
X *
X * Before we address the two conditions, we need some background information
X * on two symbols, Legendre and Jacobi. In Ref 2, pp 278, 284-285, we find
X * the following definitions of J(a,p) and L(a,n):
X *
X * The Legendre symbol L(a,p) takes the value:
X *
X * L(a,p) == 1 => a is a quadratic residue of p
X * L(a,p) == -1 => a is NOT a quadratic residue of p
X *
X * when
X *
X * p is prime
X * p mod 2 == 1
X * gcd(a,p) == 1
X *
X * The value x is a quadratic residue of y if there exists some integer z
X * such that:
X *
X * z^2 mod y == x
X *
X * The Jacobi symbol J(x,y) takes the value:
X *
X * J(x,y) == 1 => y is not prime, or x is a quadratic residue of y
X * J(x,y) == -1 => x is NOT a quadratic residue of y
X *
X * when
X *
X * y mod 2 == 1
X * gcd(x,y) == 1
X *
X * In the following comments on Legendre and Jacobi identities, we shall
X * assume that the arguments to the symbolic are valid over the symbol
X * definitions as stated above.
X *
X * In Ref2, pp 280-284, we find that:
X *
X * L(a,p)*L(b,p) == L(a*b,p) {A3.5}
X * J(x,y)*J(z,y) == J(x*z,y) {A3.14}
X * L(a,p) == L(p,a) * (-1)^((a-1)*(p-1)/4) {A3.8}
X * J(x,y) == J(y,x) * (-1)^((x-1)*(y-1)/4) {A3.17}
X *
X * The equality L(a,p) == J(a,p) when: {note 0}
X *
X * p is prime
X * p mod 2 == 1
X * gcd(a,p) == 1
X *
X * It can be shown that (see Ref3):
X *
X * L(a,p) == L(a mod p, p) {note 1}
X * L(z^2, p) == 1 {note 2}
X *
X * From Ref2, table 32:
X *
X * p mod 8 == +/-1 implies L(2,p) == 1 {note 3}
X * p mod 12 == +/-1 implies L(3,p) == 1 {note 4}
X *
X * Since h*2^n-1 mod 8 == -1, for n>2, note 3 implies:
X *
X * L(2, h*2^n-1) == 1 (n>2) {note 5}
X *
X * Since h=3*A, h*2^n-1 mod 12 == -1, for A>0, note 4 implies:
X *
X * L(3, h*2^n-1) == 1 {note 6}
X *
X * By use of {A3.5}, {note 2}, {note 5} and {note 6}, one can show:
X *
X * L((2^g)*(3^l)*(z^2), h*2^n-1) == 1 (g>=0,l>=0,z>0,n>2) {note 7}
X *
X * Returning to the testing of conditions, take condition 1:
X *
X * L(D, h*2^n-1) == -1 [condition 1]
X *
X * In order for J(D, h*2^n-1) to be defined, we must ensure that D
X * is not a factor of h*2^n-1. This is done by pre-screening h*2^n-1 to
X * not have small factors and selecting D less than that factor check limit.
X *
X * By use of {note 7}, we can show that when we choose D to be:
X *
X * D is square free
X * D = P*(2^f)*(3^g) (P is prime>2)
X *
X * The square free condition implies f = 0 or 1, g = 0 or 1. If f and g
X * are both 1, P must be a prime > 3.
X *
X * So given such a D value:
X *
X * L(D, h*2^n-1) == L(P*(2^g)*(3^l), h*2^n-1)
X * == L(P, h*2^n-1) * L((2^g)*(3^l), h*2^n-1) {A3.5}
X * == L(P, h*2^n-1) * 1 {note 7}
X * == L(h*2^n-1, P)*(-1)^((h*2^n-2)*(P-1)/4) {A3.8}
X * == L(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) {note 1}
X * == J(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) {note 0}
X *
X * When does J(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) take the value of -1,
X * thus satisfy [condition 1]? The answer depends on P. Now P is a prime>2,
X * thus P mod 4 == 1 or -1.
X *
X * Take P mod 4 == 1:
X *
X * P mod 4 == 1 implies (-1)^((h*2^n-2)*(P-1)/4) == 1
X *
X * Thus:
X *
X * L(D, h*2^n-1) == L(h*2^n-1 mod P, P) * (-1)^((h*2^n-2)*(P-1)/4)
X * == L(h*2^n-1 mod P, P)
X * == J(h*2^n-1 mod P, P)
X *
X * Take P mod 4 == -1:
X *
X * P mod 4 == -1 implies (-1)^((h*2^n-2)*(P-1)/4) == -1
X *
X * Thus:
X *
X * L(D, h*2^n-1) == L(h*2^n-1 mod P, P) * (-1)^((h*2^n-2)*(P-1)/4)
X * == L(h*2^n-1 mod P, P) * -1
X * == -J(h*2^n-1 mod P, P)
X *
X * Therefore [condition 1] is met if, and only if, one of the following
X * to cases are true:
X *
X * P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1
X * P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1
X *
X * Now consider [condition 2]:
X *
X * L(r, h*2^n-1) * (a^2 - b^2*D)/r == -1 [condition 2]
X *
X * We select only a, b, r and D values where:
X *
X * (a^2 - b^2*D)/r == -1
X *
X * Therefore in order for [condition 2] to be met, we must show that:
X *
X * L(r, h*2^n-1) == 1
X *
X * If we select r to be of the form:
X *
X * r == Q*(2^j)*(3^k)*(z^2) (Q == 1, j>=0, k>=0, z>0)
X *
X * then by use of {note 7}:
X *
X * L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
X * == L((2^j)*(3^k)*(z^2), h*2^n-1)
X * == 1 {note 2}
X *
X * and thus, [condition 2] is met.
X *
X * If we select r to be of the form:
X *
X * r == Q*(2^j)*(3^k)*(z^2) (Q is prime>2, j>=0, k>=0, z>0)
X *
X * then by use of {note 7}:
X *
X * L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
X * == L(Q, h*2^n-1) * L((2^j)*(3^k)*(z^2), h*2^n-1) {A3.5}
X * == L(Q, h*2^n-1) * 1 {note 2}
X * == L(h*2^n-1, Q) * (-1)^((h*2^n-2)*(Q-1)/4) {A3.8}
X * == L(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) {note 1}
X * == J(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) {note 0}
X *
X * When does J(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) take the value of 1,
X * thus satisfy [condition 2]? The answer depends on Q. Now Q is a prime>2,
X * thus Q mod 4 == 1 or -1.
X *
X * Take Q mod 4 == 1:
X *
X * Q mod 4 == 1 implies (-1)^((h*2^n-2)*(Q-1)/4) == 1
X *
X * Thus:
X *
X * L(D, h*2^n-1) == L(h*2^n-1 mod Q, Q) * (-1)^((h*2^n-2)*(Q-1)/4)
X * == L(h*2^n-1 mod Q, Q)
X * == J(h*2^n-1 mod Q, Q)
X *
X * Take Q mod 4 == -1:
X *
X * Q mod 4 == -1 implies (-1)^((h*2^n-2)*(Q-1)/4) == -1
X *
X * Thus:
X *
X * L(D, h*2^n-1) == L(h*2^n-1 mod Q, Q) * (-1)^((h*2^n-2)*(Q-1)/4)
X * == L(h*2^n-1 mod Q, Q) * -1
X * == -J(h*2^n-1 mod Q, Q)
X *
X * Therefore [condition 2] is met by selecting D = Q*(2^j)*(3^k)*(z^2),
X * where Q is prime>2, j>=0, k>=0, z>0; if and only if one of the following
X * to cases are true:
X *
X * Q mod 4 == 1 and J(h*2^n-1 mod Q, Q) == 1
X * Q mod 4 == -1 and J(h*2^n-1 mod Q, Q) == -1
X *
X ***
X *
X * In conclusion, we can compute v(1) by attempting to do the following:
X *
X * h mod 3 != 0
X *
X * we return:
X *
X * v(1) == 4
X *
X * h mod 3 == 0
X *
X * define:
X *
X * r == abs(a^2 - b^2*D)
X * alpha == (a + b*sqrt(D))^2/r
X *
X * we return:
X *
X * v(1) = alpha^1 + alpha^(-1)
X *
X * if and only if we can find a given a, b, D that obey all the
X * following selection rules:
X *
X * D is square free
X *
X * D == P*(2^f)*(3^g) (P is prime>2, f,g == 0 or 1)
X *
X * (a^2 - b^2*D)/r == -1
X *
X * r == Q*(2^j)*(3^k)*(z^2) (Q==1 or Q is prime>2, j>=0, k>=0, z>0)
X *
X * one of the following is true:
X * P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1
X * P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1
X *
X * if Q is prime, then one of the following is true:
X * Q mod 4 == 1 and J(h*2^n-1 mod Q, Q) == 1
X * Q mod 4 == -1 and J(h*2^n-1 mod Q, Q) == -1
X *
X * If we cannot find a v(1) quickly enough, then we will give up
X * testing h*2^n-1. This does not happen too often, so this hack
X * is not too bad.
X *
X ***
X *
X * input:
X * h h as in h*2^n-1
X * n n as in h*2^n-1
X *
X * output:
X * returns v(1), or -1 is there is no quick way
X */
Xdefine
Xgen_v1(h, n)
X{
X local d; /* the 'D' value to try */
X local val_mod; /* h*2^n-1 mod 'D' */
X local i;
X
X /*
X * check for case 1
X */
X if (h % 3 != 0) {
X /* v(1) is easy to compute */
X return 4;
X }
X
X /*
X * We will try all 'D' values until we find a proper v(1)
X * or run out of 'D' values.
X */
X for (i=0; i < quickmax; ++i) {
X
X /* grab our 'D' value */
X d = d_qval[i];
X
X /* compute h*2^n-1 mod 'D' quickly */
X val_mod = (h*pmod(2,n%(d-1),d)-1) % d;
X
X /*
X * if 'D' mod 4 == 1, then
X * (h*2^n-1) mod 'D' can not be a quadratic residue of 'D'
X * else
X * (h*2^n-1) mod 'D' must be a quadratic residue of 'D'
X */
X if (d%4 == 1) {
X /* D mod 4 == 1, so check for J(D, h*2^n-1) == -1 */
X if (jacobi(val_mod, d) == -1) {
X /* it worked, return the related v(1) value */
X return v1_qval[i];
X }
X } else {
X /* D mod 4 == -1, so check for J(D, h*2^n-1) == 1 */
X if (jacobi(val_mod, d) == 1) {
X /* it worked, return the related v(1) value */
X return v1_qval[i];
X }
X }
X }
X
X /*
X * This is an example of a more complex proof construction.
X * The code above will not be able to find the v(1) for:
X *
X * 81*2^81-1
X *
X * We will check with:
X *
X * v(1)=81 D=6557 a=79 b=1 r=316
X *
X * Now, D==79*83 and r=79*2^2. If we show that:
X *
X * J(h*2^n-1 mod 79, 79) == -1
X * J(h*2^n-1 mod 83, 83) == 1
X *
X * then we will satisfy [condition 1]. Observe:
X *
X * 79 mod 4 == -1 implies (-1)^((h*2^n-2)*(79-1)/4) == -1
X * 83 mod 4 == -1 implies (-1)^((h*2^n-2)*(83-1)/4) == -1
X *
X * J(D, h*2^n-1) == J(83, h*2^n-1) * J(79, h*2^n-1)
X * == J(h*2^n-1, 83) * (-1)^((h*2^n-2)*(83-1)/4) *
X * J(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4)
X * == J(h*2^n-1 mod 83, 83) * -1 *
X * J(h*2^n-1 mod 79, 79) * -1
X * == 1 * -1 *
X * -1 * -1
X * == -1
X *
X * We will also satisfy [condition 2]. Observe:
X *
X * (a^2 - b^2*D)/r == (79^2 - 1^1*6557)/316
X * == -1
X *
X * L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
X * == L(79, h*2^n-1) * L(2^2, h*2^n-1)
X * == L(79, h*2^n-1) * 1
X * == L(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4)
X * == L(h*2^n-1, 79) * -1
X * == L(h*2^n-1 mod 79, 79) * -1
X * == J(h*2^n-1 mod 79, 79) * -1
X * == -1 * -1
X * == 1
X */
X if (jacobi( ((h*pmod(2,n%(79-1),79)-1)%79), 79 ) == -1 &&
X jacobi( ((h*pmod(2,n%(83-1),83)-1)%83), 83 ) == 1) {
X /* return the associated v(1)=81 */
X return 81;
X }
X
X /* no quick and dirty v(1), so return -1 */
X return -1;
X}
X
X/*
X * ldebug - print a debug statement
X *
X * input:
X * funct name of calling function
X * str string to print
X */
Xdefine
Xldebug(funct, str)
X{
X if (lib_debug > 0) {
X print "DEBUG:", funct:":", str;
X }
X return;
X}
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X print "lucas(h, n) defined";
X}
SHAR_EOF
echo "File calc2.9.0/lib/lucas.cal is complete"
chmod 0644 calc2.9.0/lib/lucas.cal || echo "restore of calc2.9.0/lib/lucas.cal fails"
set `wc -c calc2.9.0/lib/lucas.cal`;Sum=$1
if test "$Sum" != "27896"
then echo original size 27896, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/lucas_chk.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/lucas_chk.cal &&
X/*
X * Copyright (c) 1993 Landon Curt Noll
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * By: Landon Curt Noll
X * chongo@toad.com -or- ...!{pyramid,sun,uunet}!hoptoad!chongo
X *
X *
X * primes of the form h*2^n-1 for 1<=h<200 and 1<=n<1000
X *
X * For all 0 <= i < prime_cnt, h_p[i]*2^n_p[i]-1 is prime.
X *
X * These values were taken from:
X *
X * "Prime numbers and Computer Methods for Factorization", by Hans Riesel,
X * Birkhauser, 1985, pp 384-387.
X *
X * This routine assumes that the file "lucas.cal" has been loaded.
X *
X * NOTE: There are several errors in Riesel's table that have been corrected
X * in this file:
X *
X * 193*2^87-1 is prime
X * 193*2^97-1 is NOT prime
X * 199*2^211-1 is prime
X * 199*2^221-1 is NOT prime
X */
X
Xstatic prime_cnt = 1145; /* number of primes in the list */
X
X/* h = prime parameters */
Xstatic mat h_p[prime_cnt] = {
X 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, /* element 0 */
X 1, 1, 1, 1, 3, 3, 3, 3, 3, 3,
X 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
X 3, 3, 3, 3, 3, 3, 3, 3, 3, 5,
X 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
X 5, 5, 5, 5, 5, 5, 7, 7, 7, 7,
X 7, 7, 7, 7, 9, 9, 9, 9, 9, 9,
X 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
X 9, 9, 9, 11, 11, 11, 11, 11, 11, 11,
X 11, 11, 11, 13, 13, 13, 13, 13, 13, 15,
X 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, /* 100 */
X 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
X 15, 15, 17, 17, 17, 17, 17, 17, 17, 17,
X 17, 17, 17, 17, 17, 17, 17, 17, 17, 17,
X 17, 17, 19, 19, 19, 19, 19, 19, 19, 19,
X 19, 19, 19, 19, 19, 19, 19, 19, 19, 19,
X 19, 19, 21, 21, 21, 21, 21, 21, 21, 21,
X 21, 21, 21, 21, 21, 21, 21, 21, 23, 23,
X 23, 23, 23, 23, 23, 23, 23, 25, 25, 25,
X 25, 25, 25, 25, 25, 25, 25, 25, 25, 25,
X 25, 25, 25, 27, 27, 27, 27, 27, 27, 27, /* 200 */
X 27, 27, 27, 27, 27, 27, 27, 27, 27, 27,
X 27, 27, 27, 27, 27, 27, 27, 29, 29, 29,
X 29, 29, 31, 31, 31, 31, 31, 31, 31, 31,
X 31, 31, 31, 31, 31, 31, 31, 31, 31, 31,
X 33, 33, 33, 33, 33, 33, 33, 33, 33, 33,
X 33, 33, 33, 33, 33, 33, 33, 33, 33, 33,
X 33, 33, 33, 33, 35, 35, 35, 35, 35, 35,
X 35, 35, 35, 35, 35, 35, 35, 35, 35, 35,
X 35, 37, 39, 39, 39, 39, 39, 39, 39, 39,
X 39, 41, 41, 41, 41, 41, 41, 41, 41, 41, /* 300 */
X 41, 41, 41, 41, 43, 43, 43, 43, 43, 45,
X 45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
X 45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
X 45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
X 45, 45, 45, 45, 45, 47, 47, 47, 47, 49,
X 49, 49, 49, 49, 49, 49, 49, 49, 49, 49,
X 49, 49, 49, 49, 49, 49, 51, 51, 51, 51,
X 51, 51, 51, 51, 51, 51, 51, 51, 51, 51,
X 51, 53, 53, 53, 53, 53, 53, 53, 53, 53,
X 53, 55, 55, 55, 55, 55, 55, 55, 55, 55, /* 400 */
X 55, 55, 55, 55, 55, 55, 55, 55, 55, 55,
X 57, 57, 57, 57, 57, 57, 57, 57, 57, 57,
X 57, 57, 57, 57, 57, 57, 57, 57, 59, 59,
X 59, 59, 59, 59, 61, 61, 61, 61, 61, 61,
X 61, 61, 61, 61, 61, 61, 61, 61, 61, 61,
X 61, 63, 63, 63, 63, 63, 63, 63, 63, 63,
X 63, 63, 63, 63, 63, 63, 63, 63, 63, 63,
X 63, 63, 63, 63, 65, 65, 65, 65, 65, 65,
X 65, 65, 65, 65, 65, 65, 65, 65, 65, 65,
X 65, 65, 67, 67, 67, 67, 67, 67, 67, 67, /* 500 */
X 69, 69, 69, 69, 69, 69, 69, 69, 69, 69,
X 69, 69, 69, 69, 69, 69, 69, 69, 69, 69,
X 69, 69, 69, 69, 69, 69, 69, 69, 69, 69,
X 69, 69, 71, 71, 71, 73, 73, 73, 73, 73,
X 73, 75, 75, 75, 75, 75, 75, 75, 75, 75,
X 75, 75, 75, 75, 75, 75, 75, 75, 75, 75,
X 75, 75, 75, 75, 75, 75, 75, 77, 77, 77,
X 77, 77, 77, 77, 77, 77, 77, 77, 77, 79,
X 79, 79, 79, 79, 79, 79, 79, 79, 79, 79,
X 81, 81, 81, 81, 81, 81, 81, 81, 81, 81, /* 600 */
X 81, 81, 81, 83, 83, 83, 83, 83, 83, 83,
X 83, 83, 83, 83, 83, 83, 83, 83, 83, 83,
X 83, 83, 83, 83, 83, 85, 85, 85, 85, 85,
X 85, 85, 85, 85, 87, 87, 87, 87, 87, 87,
X 87, 87, 87, 87, 87, 87, 87, 87, 87, 87,
X 87, 87, 87, 87, 87, 87, 89, 89, 89, 89,
X 89, 89, 89, 89, 89, 91, 91, 91, 91, 91,
X 91, 91, 91, 91, 91, 91, 91, 91, 91, 91,
X 91, 91, 91, 91, 91, 91, 91, 93, 93, 93,
X 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, /* 700 */
X 93, 93, 93, 93, 93, 95, 95, 95, 95, 95,
X 95, 95, 95, 95, 95, 97, 97, 97, 97, 97,
X 99, 99, 99, 99, 99, 99, 99, 99, 99, 99,
X 99, 99, 99, 99, 99, 99, 101, 101, 101, 101,
X 103, 103, 103, 103, 103, 103, 103, 103, 103, 103,
X 103, 103, 103, 105, 105, 105, 105, 105, 105, 105,
X 105, 105, 105, 105, 105, 105, 105, 105, 105, 105,
X 105, 105, 107, 107, 107, 107, 107, 107, 107, 107,
X 107, 107, 107, 107, 107, 107, 109, 109, 109, 109,
X 109, 113, 113, 113, 113, 113, 113, 113, 113, 113, /* 800 */
X 113, 115, 115, 115, 115, 115, 115, 115, 115, 115,
X 115, 115, 115, 115, 115, 115, 115, 119, 119, 119,
X 119, 119, 119, 119, 119, 121, 121, 121, 121, 121,
X 121, 121, 121, 121, 121, 121, 121, 125, 125, 125,
X 125, 125, 125, 127, 127, 131, 131, 131, 131, 131,
X 131, 131, 131, 131, 131, 133, 133, 133, 133, 133,
X 133, 133, 133, 133, 133, 133, 133, 133, 137, 137,
X 137, 137, 139, 139, 139, 139, 139, 139, 139, 139,
X 139, 139, 139, 139, 139, 139, 139, 139, 139, 139,
X 139, 139, 139, 139, 139, 139, 139, 139, 139, 143, /* 900 */
X 143, 143, 143, 143, 143, 143, 143, 143, 143, 143,
X 143, 143, 143, 143, 143, 143, 143, 143, 143, 143,
X 143, 143, 143, 145, 145, 145, 145, 145, 145, 145,
X 145, 145, 145, 145, 149, 149, 149, 149, 149, 149,
X 149, 151, 151, 151, 155, 155, 155, 155, 155, 155,
X 155, 155, 155, 155, 155, 155, 157, 157, 157, 157,
X 157, 157, 157, 157, 157, 161, 161, 161, 161, 161,
X 161, 161, 161, 161, 161, 161, 161, 161, 161, 161,
X 163, 163, 163, 163, 167, 167, 167, 167, 167, 167,
X 167, 167, 167, 167, 167, 167, 169, 169, 169, 169, /* 1000 */
X 169, 169, 169, 169, 169, 169, 169, 169, 173, 173,
X 173, 173, 173, 173, 173, 173, 173, 173, 173, 173,
X 173, 173, 173, 173, 175, 175, 175, 175, 175, 175,
X 175, 175, 175, 175, 175, 175, 175, 175, 175, 175,
X 179, 179, 179, 181, 181, 181, 181, 181, 181, 181,
X 181, 181, 181, 181, 181, 181, 181, 181, 181, 181,
X 181, 181, 181, 181, 181, 181, 181, 181, 185, 185,
X 185, 185, 185, 185, 185, 185, 185, 185, 187, 187,
X 187, 187, 187, 191, 193, 193, 193, 193, 193, 193,
X 193, 197, 197, 197, 197, 197, 197, 197, 197, 197, /* 1100 */
X 197, 197, 197, 197, 197, 197, 197, 197, 197, 199,
X 199, 199, 199, 199, 199, 199, 199, 199, 199, 199,
X 199, 199, 199, 199, 199, 199, 199, 199, 199, 199,
X 199, 199, 199, 199, 199
X};
X
X
X/* n (exponent) prime parameters */
Xstatic mat n_p[prime_cnt] = {
X 2, 3, 5, 7, 13, 17, 19, 31, 61, 89, /* element 0 */
X 107, 127, 521, 607, 1, 2, 3, 4, 6, 7,
X 11, 18, 34, 38, 43, 55, 64, 76, 94, 103,
X 143, 206, 216, 306, 324, 391, 458, 470, 827, 2,
X 4, 8, 10, 12, 14, 18, 32, 48, 54, 72,
X 148, 184, 248, 270, 274, 420, 1, 5, 9, 17,
X 21, 29, 45, 177, 1, 3, 7, 13, 15, 21,
X 43, 63, 99, 109, 159, 211, 309, 343, 415, 469,
X 781, 871, 939, 2, 26, 50, 54, 126, 134, 246,
X 354, 362, 950, 3, 7, 23, 287, 291, 795, 1,
X 2, 4, 5, 10, 14, 17, 31, 41, 73, 80, /* 100 */
X 82, 116, 125, 145, 157, 172, 202, 224, 266, 289,
X 293, 463, 2, 4, 6, 16, 20, 36, 54, 60,
X 96, 124, 150, 252, 356, 460, 612, 654, 664, 698,
X 702, 972, 1, 3, 5, 21, 41, 49, 89, 133,
X 141, 165, 189, 293, 305, 395, 651, 665, 771, 801,
X 923, 953, 1, 2, 3, 7, 10, 13, 18, 27,
X 37, 51, 74, 157, 271, 458, 530, 891, 4, 6,
X 12, 46, 72, 244, 264, 544, 888, 3, 9, 11,
X 17, 23, 35, 39, 75, 105, 107, 155, 215, 335,
X 635, 651, 687, 1, 2, 4, 5, 8, 10, 14, /* 200 */
X 28, 37, 38, 70, 121, 122, 160, 170, 253, 329,
X 362, 454, 485, 500, 574, 892, 962, 4, 16, 76,
X 148, 184, 1, 5, 7, 11, 13, 23, 33, 35,
X 37, 47, 115, 205, 235, 271, 409, 739, 837, 887,
X 2, 3, 6, 8, 10, 22, 35, 42, 43, 46,
X 56, 91, 102, 106, 142, 190, 208, 266, 330, 360,
X 382, 462, 503, 815, 2, 6, 10, 20, 44, 114,
X 146, 156, 174, 260, 306, 380, 654, 686, 702, 814,
X 906, 1, 3, 24, 105, 153, 188, 605, 795, 813,
X 839, 2, 10, 14, 18, 50, 114, 122, 294, 362, /* 300 */
X 554, 582, 638, 758, 7, 31, 67, 251, 767, 1,
X 2, 3, 4, 5, 6, 8, 9, 14, 15, 16,
X 22, 28, 29, 36, 37, 54, 59, 85, 93, 117,
X 119, 161, 189, 193, 256, 308, 322, 327, 411, 466,
X 577, 591, 902, 928, 946, 4, 14, 70, 78, 1,
X 5, 7, 9, 13, 15, 29, 33, 39, 55, 81,
X 95, 205, 279, 581, 807, 813, 1, 9, 10, 19,
X 22, 57, 69, 97, 141, 169, 171, 195, 238, 735,
X 885, 2, 6, 8, 42, 50, 62, 362, 488, 642,
X 846, 1, 3, 5, 7, 15, 33, 41, 57, 69, /* 400 */
X 75, 77, 131, 133, 153, 247, 305, 351, 409, 471,
X 1, 2, 4, 5, 8, 10, 20, 22, 25, 26,
X 32, 44, 62, 77, 158, 317, 500, 713, 12, 16,
X 72, 160, 256, 916, 3, 5, 9, 13, 17, 19,
X 25, 39, 63, 67, 75, 119, 147, 225, 419, 715,
X 895, 2, 3, 8, 11, 14, 16, 28, 32, 39,
X 66, 68, 91, 98, 116, 126, 164, 191, 298, 323,
X 443, 714, 758, 759, 4, 6, 12, 22, 28, 52,
X 78, 94, 124, 162, 174, 192, 204, 304, 376, 808,
X 930, 972, 5, 9, 21, 45, 65, 77, 273, 677, /* 500 */
X 1, 4, 5, 7, 9, 11, 13, 17, 19, 23,
X 29, 37, 49, 61, 79, 99, 121, 133, 141, 164,
X 173, 181, 185, 193, 233, 299, 313, 351, 377, 540,
X 569, 909, 2, 14, 410, 7, 11, 19, 71, 79,
X 131, 1, 3, 5, 6, 18, 19, 20, 22, 28,
X 29, 39, 43, 49, 75, 85, 92, 111, 126, 136,
X 159, 162, 237, 349, 381, 767, 969, 2, 4, 14,
X 26, 58, 60, 64, 100, 122, 212, 566, 638, 1,
X 3, 7, 15, 43, 57, 61, 75, 145, 217, 247,
X 3, 5, 11, 17, 21, 27, 81, 101, 107, 327, /* 600 */
X 383, 387, 941, 2, 4, 8, 10, 14, 18, 22,
X 24, 26, 28, 36, 42, 58, 64, 78, 158, 198,
X 206, 424, 550, 676, 904, 5, 11, 71, 113, 115,
X 355, 473, 563, 883, 1, 2, 8, 9, 10, 12,
X 22, 29, 32, 50, 57, 69, 81, 122, 138, 200,
X 296, 514, 656, 682, 778, 881, 4, 8, 12, 24,
X 48, 52, 64, 84, 96, 1, 3, 9, 13, 15,
X 17, 19, 23, 47, 57, 67, 73, 77, 81, 83,
X 191, 301, 321, 435, 867, 869, 917, 3, 4, 7,
X 10, 15, 18, 19, 24, 27, 39, 60, 84, 111, /* 700 */
X 171, 192, 222, 639, 954, 2, 6, 26, 32, 66,
X 128, 170, 288, 320, 470, 1, 9, 45, 177, 585,
X 1, 4, 5, 7, 8, 11, 19, 25, 28, 35,
X 65, 79, 212, 271, 361, 461, 10, 18, 54, 70,
X 3, 7, 11, 19, 63, 75, 95, 127, 155, 163,
X 171, 283, 563, 2, 3, 5, 6, 8, 9, 25,
X 32, 65, 113, 119, 155, 177, 299, 335, 426, 462,
X 617, 896, 10, 12, 18, 24, 28, 40, 90, 132,
X 214, 238, 322, 532, 858, 940, 9, 149, 177, 419,
X 617, 8, 14, 74, 80, 274, 334, 590, 608, 614, /* 800 */
X 650, 1, 3, 11, 13, 19, 21, 31, 49, 59,
X 69, 73, 115, 129, 397, 623, 769, 12, 16, 52,
X 160, 192, 216, 376, 436, 1, 3, 21, 27, 37,
X 43, 91, 117, 141, 163, 373, 421, 2, 4, 44,
X 182, 496, 904, 25, 113, 2, 14, 34, 38, 42,
X 78, 90, 178, 778, 974, 3, 11, 15, 19, 31,
X 59, 75, 103, 163, 235, 375, 615, 767, 2, 18,
X 38, 62, 1, 5, 7, 9, 15, 19, 21, 35,
X 37, 39, 41, 49, 69, 111, 115, 141, 159, 181,
X 201, 217, 487, 567, 677, 765, 811, 841, 917, 2, /* 900 */
X 4, 6, 8, 12, 18, 26, 32, 34, 36, 42,
X 60, 78, 82, 84, 88, 154, 174, 208, 256, 366,
X 448, 478, 746, 5, 13, 15, 31, 77, 151, 181,
X 245, 445, 447, 883, 4, 16, 48, 60, 240, 256,
X 304, 5, 221, 641, 2, 8, 14, 16, 44, 46,
X 82, 172, 196, 254, 556, 806, 1, 5, 33, 121,
X 125, 305, 445, 473, 513, 2, 6, 18, 22, 34,
X 54, 98, 122, 146, 222, 306, 422, 654, 682, 862,
X 3, 31, 63, 303, 4, 6, 8, 10, 16, 32,
X 38, 42, 52, 456, 576, 668, 1, 5, 11, 17, /* 1000 */
X 67, 137, 157, 203, 209, 227, 263, 917, 2, 4,
X 6, 16, 32, 50, 76, 80, 96, 104, 162, 212,
X 230, 260, 480, 612, 1, 3, 9, 21, 23, 41,
X 47, 57, 69, 83, 193, 249, 291, 421, 433, 997,
X 8, 68, 108, 3, 5, 7, 9, 11, 17, 23,
X 31, 35, 43, 47, 83, 85, 99, 101, 195, 267,
X 281, 363, 391, 519, 623, 653, 673, 701, 2, 6,
X 10, 18, 26, 40, 46, 78, 230, 542, 1, 17,
X 21, 53, 253, 226, 3, 15, 27, 63, 87, 135,
X 543, 2, 16, 20, 22, 40, 82, 112, 178, 230, /* 1100 */
X 302, 304, 328, 374, 442, 472, 500, 580, 694, 1,
X 5, 7, 15, 19, 23, 25, 27, 43, 65, 99,
X 125, 141, 165, 201, 211, 331, 369, 389, 445, 461,
X 463, 467, 513, 583, 835
X};
X
X
Xglobal lib_debug; /* from lucas.cal */
X
X
X/*
X * lucas_chk - check the lucas function on known primes
X *
X * This function tests entries in the above h_p, n_p table
X * when n_p is below a given limit.
X *
X * input:
X * high_n skip tests on n_p[i] > high_n
X *
X * returns:
X * 1 all is ok
X * 0 something went wrong
X */
Xdefine
Xlucas_chk(high_n)
X{
X local i; /* index */
X local result; /* 0 => non-prime, 1 => prime, -1 => bad test */
X local error; /* number of errors and bad tests found */
X
X /*
X * firewall
X */
X if (!isint(high_n)) {
X ldebug("test_lucas", "high_n is non-int");
X quit "FATAL: bad args: high_n must be an integer";
X }
X
X /*
X * scan thru the above prime table
X */
X error = 0;
X for (i=0; i < prime_cnt; ++i) {
X
X /* skip primes where h>=2^n */
X if (highbit(h_p[i]) >= n_p[i]) {
X if (lib_debug > 0) {
X print "h>=2^n skip:", h_p[i]:"*2^":n_p[i]:"-1";
X }
X continue;
X }
X
X /* test the prime if it is small enough */
X if (n_p[i] <= high_n) {
X
X /* test the table value */
X result = lucas(h_p[i], n_p[i]);
X
X /* report the test */
X if (result == 0) {
X print "ERROR, bad primality test of",\
X h_p[i]:"*2^":n_p[i]:"-1";
X ++error;
X } else if (result == 1) {
X print h_p[i]:"*2^":n_p[i]:"-1 is prime";
X } else if (result == -1) {
X print "ERROR, failed to compute v(1) for",\
X h_p[i]:"*2^":n_p[i]:"-1";
X ++error;
X } else {
X print "ERROR, bogus return value:", result;
X ++error;
X }
X }
X }
X
X /* return the full status */
X if (error == 0) {
X print "lucas_chk(":high_n:") passed";
X return 1;
X } else if (error == 1) {
X print "lucas_chk(":high_n:") failed", error, "test";
X return 0;
X } else {
X print "lucas_chk(":high_n:") failed", error, "tests";
X return 0;
X }
X}
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X print "lucas_chk(high_n) defined";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/lucas_chk.cal || echo "restore of calc2.9.0/lib/lucas_chk.cal fails"
set `wc -c calc2.9.0/lib/lucas_chk.cal`;Sum=$1
if test "$Sum" != "13089"
then echo original size 13089, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/lucas_tbl.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/lucas_tbl.cal &&
X/*
X * Copyright (c) 1993 Landon Curt Noll
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * By: Landon Curt Noll
X * chongo@toad.com -or- ...!{pyramid,sun,uunet}!hoptoad!chongo
X *
X *
X * Lucasian criteria for primality
X *
X * The following table is taken from:
X *
X * "Lucasian Criteria for the Primality of N=h*2^n-1", by Hans Riesel,
X * Mathematics of Computation, Vol 23 #108, p 872.
X *
X * The index of the *_val[] arrays correspond to the v(1) values found
X * in the table. That is, for v(1) == x:
X *
X * D == d_val[x]
X * a == a_val[x]
X * b == b_val[x]
X * r == r_val[x] (r == abs(a^2 - b^2*D))
X *
X *
X * Note that when *_val[i] is not a number, the related v(1) value
X * is not found in Table 1.
X */
X
Xtrymax = 100;
Xmat d_val[trymax+1];
Xmat a_val[trymax+1];
Xmat b_val[trymax+1];
Xmat r_val[trymax+1];
X/* v1= 0 INVALID */
X/* v1= 1 INVALID */
X/* v1= 2 INVALID */
Xd_val[ 3]= 5; a_val[ 3]= 1; b_val[ 3]=1; r_val[ 3]=4;
Xd_val[ 4]= 3; a_val[ 4]= 1; b_val[ 4]=1; r_val[ 4]=2;
Xd_val[ 5]= 21; a_val[ 5]= 3; b_val[ 5]=1; r_val[ 5]=12;
Xd_val[ 6]= 2; a_val[ 6]= 1; b_val[ 6]=1; r_val[ 6]=1;
X/* v1= 7 INVALID */
Xd_val[ 8]= 15; a_val[ 8]= 3; b_val[ 8]=1; r_val[ 8]=6;
Xd_val[ 9]= 77; a_val[ 9]= 7; b_val[ 9]=1; r_val[ 9]=28;
Xd_val[10]= 6; a_val[10]= 2; b_val[10]=1; r_val[10]=2;
Xd_val[11]= 13; a_val[11]= 3; b_val[11]=1; r_val[11]=4;
Xd_val[12]= 35; a_val[12]= 5; b_val[12]=1; r_val[12]=10;
Xd_val[13]= 165; a_val[13]=11; b_val[13]=1; r_val[13]=44;
X/* v1=14 INVALID */
Xd_val[15]= 221; a_val[15]=13; b_val[15]=1; r_val[15]=52;
Xd_val[16]= 7; a_val[16]= 3; b_val[16]=1; r_val[16]=2;
Xd_val[17]= 285; a_val[17]=15; b_val[17]=1; r_val[17]=60;
X/* v1=18 INVALID */
Xd_val[19]= 357; a_val[19]=17; b_val[19]=1; r_val[19]=68;
Xd_val[20]= 11; a_val[20]= 3; b_val[20]=1; r_val[20]=2;
Xd_val[21]= 437; a_val[21]=19; b_val[21]=1; r_val[21]=76;
Xd_val[22]= 30; a_val[22]= 5; b_val[22]=1; r_val[22]=5;
X/* v1=23 INVALID */
Xd_val[24]= 143; a_val[24]=11; b_val[24]=1; r_val[24]=22;
Xd_val[25]= 69; a_val[25]= 9; b_val[25]=1; r_val[25]=12;
Xd_val[26]= 42; a_val[26]= 6; b_val[26]=1; r_val[26]=6;
Xd_val[27]= 29; a_val[27]= 5; b_val[27]=1; r_val[27]=4;
Xd_val[28]= 195; a_val[28]=13; b_val[28]=1; r_val[28]=26;
Xd_val[29]= 93; a_val[29]= 9; b_val[29]=1; r_val[29]=12;
Xd_val[30]= 14; a_val[30]= 4; b_val[30]=1; r_val[30]=2;
Xd_val[31]= 957; a_val[31]=29; b_val[31]=1; r_val[31]=116;
Xd_val[32]= 255; a_val[32]=15; b_val[32]=1; r_val[32]=30;
Xd_val[33]=1085; a_val[33]=31; b_val[33]=1; r_val[33]=124;
X/* v1=34 INVALID */
Xd_val[35]=1221; a_val[35]=33; b_val[35]=1; r_val[35]=132;
Xd_val[36]= 323; a_val[36]=17; b_val[36]=1; r_val[36]=34;
Xd_val[37]=1365; a_val[37]=35; b_val[37]=1; r_val[37]=140;
Xd_val[38]= 10; a_val[38]= 3; b_val[38]=1; r_val[38]=1;
Xd_val[39]=1517; a_val[39]=37; b_val[39]=1; r_val[39]=148;
Xd_val[40]= 399; a_val[40]=19; b_val[40]=1; r_val[40]=38;
Xd_val[41]=1677; a_val[41]=39; b_val[41]=1; r_val[41]=156;
Xd_val[42]= 110; a_val[42]=10; b_val[42]=1; r_val[42]=10;
Xd_val[43]= 205; a_val[43]=15; b_val[43]=1; r_val[43]=20;
Xd_val[44]= 483; a_val[44]=21; b_val[44]=1; r_val[44]=42;
Xd_val[45]=2021; a_val[45]=43; b_val[45]=1; r_val[45]=172;
Xd_val[46]= 33; a_val[46]= 6; b_val[46]=1; r_val[46]=3;
X/* v1=47 INVALID */
Xd_val[48]= 23; a_val[48]= 5; b_val[48]=1; r_val[48]=2;
Xd_val[49]=2397; a_val[49]=47; b_val[49]=1; r_val[49]=188;
Xd_val[50]= 39; a_val[50]= 6; b_val[50]=1; r_val[50]=3;
Xd_val[51]= 53; a_val[51]= 7; b_val[51]=1; r_val[51]=4;
X/* v1=52 INVALID */
Xd_val[53]=2805; a_val[53]=51; b_val[53]=1; r_val[53]=204;
Xd_val[54]= 182; a_val[54]=13; b_val[54]=1; r_val[54]=13;
Xd_val[55]=3021; a_val[55]=53; b_val[55]=1; r_val[55]=212;
Xd_val[56]= 87; a_val[56]= 9; b_val[56]=1; r_val[56]=6;
Xd_val[57]=3245; a_val[57]=55; b_val[57]=1; r_val[57]=220;
Xd_val[58]= 210; a_val[58]=14; b_val[58]=1; r_val[58]=14;
Xd_val[59]=3477; a_val[59]=57; b_val[59]=1; r_val[59]=228;
Xd_val[60]= 899; a_val[60]=29; b_val[60]=1; r_val[60]=58;
Xd_val[61]= 413; a_val[61]=21; b_val[61]=1; r_val[61]=28;
X/* v1=62 INVALID */
Xd_val[63]=3965; a_val[63]=61; b_val[63]=1; r_val[63]=244;
Xd_val[64]=1023; a_val[64]=31; b_val[64]=1; r_val[64]=62;
Xd_val[65]= 469; a_val[65]=21; b_val[65]=1; r_val[65]=28;
Xd_val[66]= 17; a_val[66]= 4; b_val[66]=1; r_val[66]=1;
Xd_val[67]=4485; a_val[67]=65; b_val[67]=1; r_val[67]=260;
Xd_val[68]=1155; a_val[68]=33; b_val[68]=1; r_val[68]=66;
Xd_val[69]=4757; a_val[69]=67; b_val[69]=1; r_val[69]=268;
Xd_val[70]= 34; a_val[70]= 6; b_val[70]=1; r_val[70]=2;
Xd_val[71]=5037; a_val[71]=69; b_val[71]=1; r_val[71]=276;
Xd_val[72]=1295; a_val[72]=35; b_val[72]=1; r_val[72]=70;
Xd_val[73]= 213; a_val[73]=15; b_val[73]=1; r_val[73]=12;
Xd_val[74]= 38; a_val[74]= 6; b_val[74]=1; r_val[74]=2;
Xd_val[75]=5621; a_val[75]=73; b_val[75]=1; r_val[75]=292;
Xd_val[76]=1443; a_val[76]=37; b_val[76]=1; r_val[76]=74;
Xd_val[77]= 237; a_val[77]=15; b_val[77]=1; r_val[77]=12;
Xd_val[78]= 95; a_val[78]=10; b_val[78]=1; r_val[78]=5;
X/* v1=79 INVALID */
Xd_val[80]=1599; a_val[80]=39; b_val[80]=1; r_val[80]=78;
Xd_val[81]=6557; a_val[81]=79; b_val[81]=1; r_val[81]=316;
Xd_val[82]= 105; a_val[82]=10; b_val[82]=1; r_val[82]=5;
Xd_val[83]= 85; a_val[83]= 9; b_val[83]=1; r_val[83]=4;
Xd_val[84]=1763; a_val[84]=41; b_val[84]=1; r_val[84]=82;
Xd_val[85]=7221; a_val[85]=83; b_val[85]=1; r_val[85]=332;
Xd_val[86]= 462; a_val[86]=21; b_val[86]=1; r_val[86]=21;
Xd_val[87]=7565; a_val[87]=85; b_val[87]=1; r_val[87]=340;
Xd_val[88]= 215; a_val[88]=15; b_val[88]=1; r_val[88]=10;
Xd_val[89]=7917; a_val[89]=87; b_val[89]=1; r_val[89]=348;
Xd_val[90]= 506; a_val[90]=22; b_val[90]=1; r_val[90]=22;
Xd_val[91]=8277; a_val[91]=89; b_val[91]=1; r_val[91]=356;
Xd_val[92]= 235; a_val[92]=15; b_val[92]=1; r_val[92]=10;
Xd_val[93]=8645; a_val[93]=91; b_val[93]=1; r_val[93]=364;
Xd_val[94]= 138; a_val[94]=12; b_val[94]=1; r_val[94]=6;
Xd_val[95]=9021; a_val[95]=93; b_val[95]=1; r_val[95]=372;
Xd_val[96]= 47; a_val[96]= 7; b_val[96]=1; r_val[96]=2;
Xd_val[97]=1045; a_val[97]=33; b_val[97]=1; r_val[97]=44;
X/* v1=98 INVALID */
Xd_val[99]=9797; a_val[99]=97; b_val[99]=1; r_val[99]=388;
Xd_val[100]= 51; a_val[100]= 7; b_val[100]=1; r_val[100]=2;
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X print "d_val[100] defined";
X print "a_val[100] defined";
X print "b_val[100] defined";
X print "r_val[100] defined";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/lucas_tbl.cal || echo "restore of calc2.9.0/lib/lucas_tbl.cal fails"
set `wc -c calc2.9.0/lib/lucas_tbl.cal`;Sum=$1
if test "$Sum" != "6646"
then echo original size 6646, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/mersenne.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/mersenne.cal &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Perform a primality test of 2^p-1, for prime p>1.
X */
X
Xdefine mersenne(p)
X{
X local u, i, p_mask;
X
X /* firewall */
X if (! isint(p))
X quit "p is not an integer";
X
X /* two is a special case */
X if (p == 2)
X return 1;
X
X /* if p is not prime, then 2^p-1 is not prime */
X if (! ptest(p,10))
X return 0;
X
X /* calculate 2^p-1 for later mods */
X p_mask = 2^p - 1;
X
X /* lltest: u(i+1) = u(i)^2 - 2 mod 2^p-1 */
X u = 4;
X for (i = 2; i < p; ++i) {
X u = u^2 - 2;
X u = u&p_mask + u>>p;
X if (u > p_mask)
X u = u&p_mask + 1;
X }
X
X /* 2^p-1 is prime iff u(p) = 0 mod 2^p-1 */
X return (u == p_mask);
X}
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X print "mersenne(p) defined";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/mersenne.cal || echo "restore of calc2.9.0/lib/mersenne.cal fails"
set `wc -c calc2.9.0/lib/mersenne.cal`;Sum=$1
if test "$Sum" != "833"
then echo original size 833, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/mod.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/mod.cal &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Routines to handle numbers modulo a specified number.
X * a (mod N)
X */
X
Xobj mod {a}; /* definition of the object */
X
Xglobal mod_value = 100; /* modulus value (value of N) */
X
X
Xdefine mod(a)
X{
X local obj mod x;
X
X if (!isreal(a) || !isint(a))
X quit "Bad argument for mod function";
X x.a = a % mod_value;
X return x;
X}
X
X
Xdefine mod_print(a)
X{
X if (digits(mod_value) <= 20)
X print a.a, "(mod", mod_value : ")" :;
X else
X print a.a, "(mod N)" :;
X}
X
X
Xdefine mod_one()
X{
X return mod(1);
X}
X
X
Xdefine mod_cmp(a, b)
X{
X if (isnum(a))
X return (a % mod_value) != b.a;
X if (isnum(b))
X return (b % mod_value) != a.a;
X return a.a != b.a;
X}
X
X
Xdefine mod_rel(a, b)
X{
X if (isnum(a))
X a = mod(a);
X if (isnum(b))
X b = mod(b);
X if (a.a < b.a)
X return -1;
X return a.a != b.a;
X}
X
X
Xdefine mod_add(a, b)
X{
X local obj mod x;
X
X if (isnum(b)) {
X if (!isint(b))
X quit "Adding non-integer";
X x.a = (a.a + b) % mod_value;
X return x;
X }
X if (isnum(a)) {
X if (!isint(a))
X quit "Adding non-integer";
X x.a = (a + b.a) % mod_value;
X return x;
X }
X x.a = (a.a + b.a) % mod_value;
X return x;
X}
X
X
Xdefine mod_sub(a, b)
X{
X return a + (-b);
X}
X
X
Xdefine mod_neg(a)
X{
X local obj mod x;
X
X x.a = mod_value - a.a;
X return x;
X}
X
X
Xdefine mod_mul(a, b)
X{
X local obj mod x;
X
X if (isnum(b)) {
X if (!isint(b))
X quit "Multiplying by non-integer";
X x.a = (a.a * b) % mod_value;
X return x;
X }
X if (isnum(a)) {
X if (!isint(a))
X quit "Multiplying by non-integer";
X x.a = (a * b.a) % mod_value;
X return x;
X }
X x.a = (a.a * b.a) % mod_value;
X return x;
X}
X
X
Xdefine mod_square(a)
X{
X local obj mod x;
X
X x.a = a.a^2 % mod_value;
X return x;
X}
X
X
Xdefine mod_inc(a)
X{
X local x;
X
X x = a;
X if (++x.a == mod_value)
X x.a = 0;
X return x;
X}
X
X
Xdefine mod_dec(a)
X{
X local x;
X
X x = a;
X if (--x.a < 0)
X x.a = mod_value - 1;
X return x;
X}
X
X
Xdefine mod_inv(a)
X{
X local obj mod x;
X
X x.a = minv(a.a, mod_value);
X return x;
X}
X
X
Xdefine mod_div(a, b)
X{
X local c, x, y;
X
X obj mod x, y;
X if (isnum(a))
X a = mod(a);
X if (isnum(b))
X b = mod(b);
X c = gcd(a.a, b.a);
X x.a = a.a / c;
X y.a = b.a / c;
X return x * inverse(y);
X}
X
X
Xdefine mod_pow(a, b)
X{
X local x, y, z;
X
X obj mod x;
X y = a;
X z = b;
X if (b < 0) {
X y = inverse(a);
X z = -b;
X }
X x.a = pmod(y.a, z, mod_value);
X return x;
X}
X
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X print "obj mod {a} defined";
X print "mod(a) defined";
X print "mod_print(a) defined";
X print "mod_one(a) defined";
X print "mod_cmp(a, b) defined";
X print "mod_rel(a, b) defined";
X print "mod_add(a, b) defined";
X print "mod_sub(a, b) defined";
X print "mod_mod(a, b) defined";
X print "mod_square(a) defined";
X print "mod_inc(a) defined";
X print "mod_dec(a) defined";
X print "mod_inv(a) defined";
X print "mod_div(a, b) defined";
X print "mod_pow(a, b) defined";
X print "mod_value defined";
X print "set mod_value as needed";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/mod.cal || echo "restore of calc2.9.0/lib/mod.cal fails"
set `wc -c calc2.9.0/lib/mod.cal`;Sum=$1
if test "$Sum" != "3005"
then echo original size 3005, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/nextprim.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/nextprim.cal &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Function to find the next prime (probably).
X */
X
Xdefine nextprime(n, tries)
X{
X if (isnull(tries))
X tries = 20;
X if (iseven(n))
X n++;
X while (ptest(n, tries) == 0)
X n += 2;
X return n;
X}
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X print "nextprime(n, tries) defined";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/nextprim.cal || echo "restore of calc2.9.0/lib/nextprim.cal fails"
set `wc -c calc2.9.0/lib/nextprim.cal`;Sum=$1
if test "$Sum" != "440"
then echo original size 440, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/pell.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/pell.cal &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Solve Pell's equation; Returns the solution X to: X^2 - D * Y^2 = 1.
X * Type the solution to pells equation for a particular D.
X */
X
Xdefine pell(D)
X{
X local X, Y;
X
X X = pellx(D);
X if (isnull(X)) {
X print "D=":D:" is square";
X return;
X }
X Y = isqrt((X^2 - 1) / D);
X print X : "^2 - " : D : "*" : Y : "^2 = " : X^2 - D*Y^2;
X}
X
X
X/*
X * Function to solve Pell's equation
X * Returns the solution X to:
X * X^2 - D * Y^2 = 1
X */
Xdefine pellx(D)
X{
X local R, Rp, U, Up, V, Vp, A, T, Q1, Q2, n;
X local mat ans[2,2];
X local mat tmp[2,2];
X
X R = isqrt(D);
X Vp = D - R^2;
X if (Vp == 0)
X return;
X Rp = R + R;
X U = Rp;
X Up = U;
X V = 1;
X A = 0;
X n = 0;
X ans[0,0] = 1;
X ans[1,1] = 1;
X tmp[0,1] = 1;
X tmp[1,0] = 1;
X do {
X T = V;
X V = A * (Up - U) + Vp;
X Vp = T;
X A = U // V;
X Up = U;
X U = Rp - U % V;
X tmp[0,0] = A;
X ans *= tmp;
X n++;
X } while (A != Rp);
X Q2 = ans[[1]];
X Q1 = isqrt(Q2^2 * D + 1);
X if (isodd(n)) {
X T = Q1^2 + D * Q2^2;
X Q2 = Q1 * Q2 * 2;
X Q1 = T;
X }
X return Q1;
X}
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X print "pell(D) defined";
X print "pellx(D) defined";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/pell.cal || echo "restore of calc2.9.0/lib/pell.cal fails"
set `wc -c calc2.9.0/lib/pell.cal`;Sum=$1
if test "$Sum" != "1247"
then echo original size 1247, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/pi.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/pi.cal &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Calculate pi within the specified epsilon using the quartic convergence
X * iteration.
X */
X
Xdefine qpi(epsilon)
X{
X local niter, yn, ym, tm, an, am, t, tn, sqrt2, epsilon2, count, digits;
X local bits, bits2;
X
X if (isnull(epsilon))
X epsilon = epsilon();
X digits = digits(1/epsilon);
X if (digits <= 8) { niter = 1; epsilon = 1e-8; }
X else if (digits <= 40) { niter = 2; epsilon = 1e-40; }
X else if (digits <= 170) { niter = 3; epsilon = 1e-170; }
X else if (digits <= 693) { niter = 4; epsilon = 1e-693; }
X else {
X niter = 4;
X t = 693;
X while (t < digits) {
X ++niter;
X t *= 4;
X }
X }
X epsilon2 = epsilon/(digits/10 + 1);
X digits = digits(1/epsilon2);
X sqrt2 = sqrt(2, epsilon2);
X bits = abs(ilog2(epsilon)) + 1;
X bits2 = abs(ilog2(epsilon2)) + 1;
X yn = sqrt2 - 1;
X an = 6 - 4 * sqrt2;
X tn = 2;
X for (count = 0; count < niter; count++) {
X ym = yn;
X am = an;
X tn *= 4;
X t = sqrt(sqrt(1-ym^4, epsilon2), epsilon2);
X yn = (1-t)/(1+t);
X an = (1+yn)^4*am-tn*yn*(1+yn+yn^2);
X yn = bround(yn, bits2);
X an = bround(an, bits2);
X }
X return (bround(1/an, bits));
X}
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X print "qpi(epsilon) defined";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/pi.cal || echo "restore of calc2.9.0/lib/pi.cal fails"
set `wc -c calc2.9.0/lib/pi.cal`;Sum=$1
if test "$Sum" != "1311"
then echo original size 1311, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/pollard.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/pollard.cal &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Factor using Pollard's p-1 method.
X */
X
Xdefine factor(N, B, ai, af)
X{
X local a, k, i, d;
X
X if (isnull(B))
X B = 1000;
X if (isnull(ai))
X ai = 2;
X if (isnull(af))
X af = ai + 20;
X k = lcmfact(B);
X d = lfactor(N, B);
X if (d > 1)
X return d;
X for (a = ai; a <= af; a++) {
X i = pmod(a, k, N);
X d = gcd(i - 1, N);
X if ((d > 1) && (d != N))
X return d;
X }
X return 1;
X}
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X print "factor(N, B, ai, af) defined";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/pollard.cal || echo "restore of calc2.9.0/lib/pollard.cal fails"
set `wc -c calc2.9.0/lib/pollard.cal`;Sum=$1
if test "$Sum" != "620"
then echo original size 620, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/poly.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/poly.cal &&
X/*
X * A collection of functions designed for calculations involving
X * polynomials in one variable (by Ernest W. Bowen).
X *
X * On starting the program the independent variable has identifier x
X * and name "x", i.e. the user can refer to it as x, the
X * computer displays it as "x". The name of the independent
X * variable is stored as varname, so, for example, varname = "alpha"
X * will change its name to "alpha". At any time, the independent
X * variable has only one name. For some purposes, a name like
X * "sin(t)" or "(a + b)" or "\lambda" might be useful;
X * names like "*" or "-27" are legal but might give expressions
X * that are difficult to intepret.
X *
X * Polynomial expressions may be constructed from numbers and the
X * independent variable and other polynomials by the algebraic
X * operations +, -, *, ^, and if the result is a polynomial /.
X * The operations // and % are defined to have the quotient and
X * remainder meanings as usually defined for polynomials.
X *
X * When polynomials are assigned to idenfifiers, it is convenient to
X * think of the polynomials as values. For example, p = (x - 1)^2
X * assigns to p a polynomial value in the same way as q = (7 - 1)^2
X * would assign to q a number value. As with number expressions
X * involving operations, the expression used to define the
X * polynomial is usually lost; in the above example, the normal
X * computer display for p will be x^2 - 2x + 1. Different
X * identifiers may of course have the same polynomial value.
X *
X * The polynomial we think of as a_0 + a_1 * x + ... + a_n * x^n,
X * for number coefficients a_0, a_1, ... a_n may also be
X * constructed as pol(a_0, a_1, ..., a_n). Note that here the
X * coefficients are to be in ascending power order. The independent
X * variable is pol(0,1), so to use t, say, as an identifier for
X * this, one may assign t = pol(0,1). To simultaneously specify
X * an identifier and a name for the independent variable, there is
X * the instruction var, used as in identifier = var(name). For
X * example, to use "t" in the way "x" is initially, one may give
X * the instruction t = var("t").
X *
X * There are four parameters pmode, order, iod and ims for controlling
X * the format in which polynomials are displayed.
X * The parameter pmode may have values "alg" or "list": the
X * former gives a display as an algebraic formula, while
X * the latter only lists the coefficients. Whether the terms or
X * coefficients are in ascending or descending power order is
X * controlled by order being "up" or "down". If the
X * parameter iod (for integer-only display), the polynomial
X * is expressed in terms of a polynomial whose coefficients are
X * integers with gcd = 1, the leading coefficient having positive
X * real part, with where necessary a leading multiplying integer,
X * a Gaussian integer multiplier if the coefficients are complex
X * with a common complex factor, and a trailing divisor integer.
X * If a non-zero value is assigned to the parameter ims,
X * multiplication signs will be inserted where appropriate;
X * this may be useful if the expression is to be copied to a
X * program or a string to be used with eval.
X *
X * For evaluation of polynomials the standard function is ev(p, t).
X * If p is a polynomial and t anything for which the relevant
X * operations can be performed, this returns the value of p
X * at t. The function ev(p, t) also accepts lists or matrices
X * as possible values for p; each element of p is then evaluated
X * at t. For other p, t is ignored and the value of p is returned.
X * If an identifier, a, say, is used for the polynomial, list or
X * matrix p, the definition
X * define a(t) = ev(a, t);
X * permits a(t) to be used for the value of a at t as if the
X * polynomial, list or matrix were a function. For example,
X * if a = 1 + x^2, a(2) will return the value 5, just as if
X * define a(t) = 1 + t^2;
X * had been used. However, when the polynomial definition is
X * used, changing the polynomial a will change a(t) to the value
X * of the new polynomial at t. For example,
X * after
X * L = list(x, x^2, x^3, x^4);
X define a(t) = ev(a, t);
X * the loop
X * for (i = 0; i < 4; i++)
X * print ev(L[[i]], 5);
X * may be replaced by
X * for (i = 0; i < 4; i++) {
X * a = L[[i]];
X * print a(5);
X * }
X *
X * Matrices with polynomial elements may be added, subtracted and
X * multiplied as long as the usual rules for compatibility are
X * observed. Also, matrices may be multiplied by polynomials,
X * i.e. if p is a polynomial and A a matrix whose elements
X * may be numbers or polynomials, p * A returns the matrix of
X * the same shape as A with each element multiplied by p.
X * Square matrices may also be 'substituted for the variable' in
X * polynomials, e.g. if A is an m x m matrix, and
X * p = x^2 + 3 * x + 2, ev(p, A) returns the same as
X * A^2 + 3 * A + 2 * I, where I is the unit m x m matrix.
X *
X * On starting this program, three demonstration polynomials a, b, c
X * have been defined. The functions a(t), b(t), c(t) corresponding
X * to a, b, c, and x(t) corresponding to x, have also been
X * defined, so the usual function notation can be used for
X * evaluations of a, b, c and x. For x, as long as x identifies
X * the independent variable, x(t) should return the value of t,
X * i.e. it acts as an identity function.
X *
X * Functions defined include:
X *
X * monic(a) returns the monic multiple of a, i.e., if a != 0,
X * the multiple of a with leading coefficient 1
X * conj(a) returns the complex conjugate of a
X * ispmult(a,b) returns 1 or 0 according as a is or is not
X * a polynomial multiple of b
X * pgcd(a,b) returns the monic gcd of a and b
X * pfgcd(a,b) returns a list of three polynomials (g, u, v)
X * where g = pgcd(a,b) and g = u * a + v * b.
X * plcm(a,b) returns the monic lcm of a and b
X *
X * interp(X,Y,t) returns the value at t of the polynomial given
X * by Newtonian divided difference interpolation, where
X * X is a list of x-values, Y a list of corresponding
X * y-values. If t is omitted, the interpolating
X * polynomial is returned. A y-value may be replaced by
X * list (y, y_1, y_2, ...), where y_1, y_2, ... are
X * the reduced derivatives at the corresponding x;
X * i.e. y_r is the r-th derivative divided by fact(r).
X * mdet(A) returns the determinant of the square matrix A,
X * computed by an algorithm that does not require
X * inverses; the built-in det function usually fails
X * for matrices with polynomial elements.
X * D(a,n) returns the n-th derivative of a; if n is omitted,
X * the first derivative is returned.
X *
X * A first-time user can see what the initially defined polynomials
X * a, b and c are, and experiment with the algebraic operations
X * and other functions that have been defined by giving
X * instructions like:
X * a
X * b
X * c
X * (x^2 + 1) * a
X * a^27
X * a * b
X * a % b
X * a // b
X * a(1 + x)
X * a(b)
X * conj(c)
X * g = pgcd(a, b)
X * g
X * a / g
X * D(a)
X * mat A[2,2] = {1 + x, x^2, 3, 4*x}
X * mdet(A)
X * D(A)
X * A^2
X * define A(t) = ev(A, t)
X * A(2)
X * A(1 + x)
X * define L(t) = ev(L, t)
X * L = list(x, x^2, x^3, x^4)
X * L(5)
X * a(L)
X * interp(list(0,1,2,3), list(2,3,5,7))
X * interp(list(0,1,2), list(0,list(1,0),2))
X *
X * One check on some of the functions is provided by the Cayley-Hamilton
X * theorem: if A is any m x m matrix and I the m x m unit matrix,
X * and x is pol(0,1),
X * ev(mdet(x * I - A), A)
X * should return the zero m x m matrix.
X */
X
Xobj poly {p};
X
Xdefine pol() {
X local u,i,s;
X obj poly u;
X s = list();
X for (i=1; i<= param(0); i++) append (s,param(i));
X i=size(s) -1;
X while (i>=0 && s[[i]]==0) {i--; remove(s)}
X u.p = s;
X return u;
X}
X
Xdefine ispoly(a) {
X local y;
X obj poly y;
X return istype(a,y);
X}
X
Xdefine findlist(a) {
X if (ispoly(a)) return a.p;
X if (a) return list(a);
X return list();
X}
X
Xpmode = "alg"; /* The other acceptable pmode is "list" */
Xims = 0; /* To be non-zero if multiplication signs to be inserted */
Xiod = 0; /* To be non-zero for integer-only display */
Xorder = "down" /* Determines order in which coefficients displayed */
X
Xdefine poly_print(a) {
X local f, g, t;
X if (size(a.p) == 0) {
X print 0:;
X return;
X }
X if (iod) {
X g = gcdcoeffs(a);
X t = a.p[[size(a.p) - 1]] / g;
X if (re(t) < 0) { t = -t; g = -g;}
X if (g != 1) {
X if (!isreal(t)) {
X if (im(t) > re(t)) g *= 1i;
X else if (im(t) <= -re(t)) g *= -1i;
X }
X if (isreal(g)) f = g;
X else f = gcd(re(g), im(g));
X if (num(f) != 1) {
X print num(f):;
X if (ims) print"*":;
X }
X if (!isreal(g)) {
X printf("(%d)", g/f);
X if (ims) print"*":;
X }
X if (pmode == "alg") print"(":;
X polyprint(1/g * a);
X if (pmode == "alg") print")":;
X if (den(f) > 1) print "/":den(f):;
X return;
X }
X }
X polyprint(a);
X}
X
Xdefine polyprint(a) {
X local s,n,i,c;
X s = a.p;
X n=size(s) - 1;
X if (pmode=="alg") {
X if (order == "up") {
X i = 0;
X while (!s[[i]]) i++;
X pterm (s[[i]], i);
X for (i++ ; i <= n; i++) {
X c = s[[i]];
X if (c) {
X if (isreal(c)) {
X if (c > 0) print" + ":;
X else {
X print" - ":;
X c = -c;
X }
X }
X else print " + ":;
X pterm(c,i);
X }
X }
X return;
X }
X if (order == "down") {
X pterm(s[[n]],n);
X for (i=n-1; i>=0; i--) {
X c = s[[i]];
X if (c) {
X if (isreal(c)) {
X if (c > 0) print" + ":;
X else {
X print" - ":;
X c = -c;
X }
X }
X else print " + ":;
X pterm(c,i);
X }
X }
X return;
X }
X quit "order to be up or down";
X }
X if (pmode=="list") {
X plist(s);
X return;
X }
X print pmode,:"is unknown mode";
X}
X
X
Xdefine poly_neg(a) {
X local s,i,y;
X obj poly y;
X s = a.p;
X for (i=0; i< size(s); i++) s[[i]] = -s[[i]];
X y.p = s;
X return y;
X}
X
Xdefine poly_conj(a) {
X local s,i,y;
X obj poly y;
X s = a.p;
X for (i=0; i < size(s); i++) s[[i]] = conj(s[[i]]);
X y.p = s;
X return y;
X}
X
Xdefine poly_inv(a) = pol(1)/a; /* This exists only for a of zero degree */
X
Xdefine poly_add(a,b) {
X local sa, sb, i, y;
X obj poly y;
X sa=findlist(a); sb=findlist(b);
X if (size(sa) > size(sb)) swap(sa,sb);
X for (i=0; i< size(sa); i++) sa[[i]] += sb[[i]];
X while (i < size(sb)) append (sa, sb[[i++]]);
X while (i > 0 && sa[[--i]]==0) remove (sa);
X y.p = sa;
X return y;
X}
X
Xdefine poly_sub(a,b) {
X return a + (-b);
X}
X
Xdefine poly_cmp(a,b) {
X local sa, sb;
X sa = findlist(a);
X sb=findlist(b);
X return (sa != sb);
X}
X
Xdefine poly_mul(a,b) {
X local sa,sb,i, j, y;
X if (ismat(a)) swap(a,b);
X if (ismat(b)) {
X y = b;
X for (i=matmin(b,1); i <= matmax(b,1); i++)
X for (j = matmin(b,2); j<= matmax(b,2); j++)
X y[i,j] = a * b[i,j];
X return y;
X }
X obj poly y;
X sa=findlist(a); sb=findlist(b);
X y.p = listmul(sa,sb);
X return y;
X}
X
Xdefine listmul(a,b) {
X local da,db, s, i, j, u;
X da=size(a)-1; db=size(b)-1;
X s=list();
X if (da >= 0 && db >= 0) {
X for (i=0; i<= da+db; i++) { u=0;
X for (j = max(0,i-db); j <= min(i, da); j++)
X u += a[[j]]*b[[i-j]]; append (s,u);}}
X return s;
X}
X
Xdefine ev(a,t) {
X local v, i, j;
X if (ismat(a)) {
X v = a;
X for (i = matmin(a,1); i <= matmax(a,1); i++)
X for (j = matmin(a,2); j <= matmax(a,2); j++)
X v[i,j] = ev(a[i,j], t);
X return v;
X }
X if (islist(a)) {
X v = list();
X for (i = 0; i < size(a); i++)
X append(v, ev(a[[i]], t));
X return v;
X }
X if (!ispoly(a)) return a;
X if (islist(t)) {
X v = list();
X for (i = 0; i < size(t); i++)
X append(v, ev(a, t[[i]]));
X return v;
X }
X if (ismat(t)) return evpm(a.p, t);
X return evp(a.p, t);
X}
X
Xdefine evp(s,t) {
X local n,v,i;
X n = size(s);
X if (!n) return 0;
X v = s[[n-1]];
X for (i = n - 2; i >= 0; i--) v=t * v +s[[i]];
X return v;
X}
X
Xdefine evpm(s,t) {
X local m, n, V, i, I;
X n = size(s);
X m = matmax(t,1) - matmin(t,1);
X if (matmax(t,2) - matmin(t,2) != m) quit "Non-square matrix";
X mat V[m+1, m+1];
X if (!n) return V;
X mat I[m+1, m+1];
X matfill(I, 0, 1);
X V = s[[n-1]] * I;
X for (i = n - 2; i >= 0; i--) V = t * V + s[[i]] * I;
X return V;
X}
Xpzero = pol(0);
Xx = pol(0,1);
Xvarname = "x";
Xdefine x(t) = ev(x, t);
X
Xdefine iszero(a) {
X if (ispoly(a))
X return !size(a.p);
X return a == 0;
X}
X
Xdefine isstring(a) = istype(a, " ");
X
Xdefine var(name) {
X if (!isstring(name)) quit "Argument of var is to be a string";
X varname = name;
X return pol(0,1);
X}
X
Xdefine pcoeff(a) {
X if (isreal(a)) print a:;
X else print "(":a:")":;
X}
X
Xdefine pterm(a,n) {
SHAR_EOF
echo "End of part 18"
echo "File calc2.9.0/lib/poly.cal is continued in part 19"
echo "19" > s2_seq_.tmp
exit 0
