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