14 March 1991 Procedures ( File : ProcedureList ) All procedures/functions written by Donald E. G. Malm and copyright by him. For each procedure/function we normally give its Pascal heading, together with other pertinent information. The UBASIC versions are like the Pascal versions, using the same variable names. The procedures are listed alphabetically, with the separation of the Pascal routines into the three files BODY 1,2, and 3 indicated. (The UBASIC subroutines are in individual files.) There are minor differences between the UBASIC and Pascal versions, chiefly because UBASIC has no Boolean or character variables. ************************** BODY1 ************************************* PROCEDURE BWppt1(n : number; VAR f : integer); { Baillie-Wagstaff Lucas pseudoprime test. See their paper. We assume n is odd and > 11. D is chosen by method "A". f is 1 to indicate that n is a probable prime, 0 for composite, and -1 if n is a square. Requires unit MultiP. Calls Add, Sub, Mul, Dvd, Dvdby2, Equal, Isqrt, Jacobi, and Modd. } { 6 June 1990 } The article is "LUCAS PSEUDOPRIMES", Robert Baillie and Samuel S. Wagstaff, Jr., Math. Comp., v. 35 (1980) pp. 1391-1417. PROCEDURE BWppt2(n : number; VAR f : integer); { Baillie-Wagstaff strong Lucas pseudoprime test. See their paper. We assume n is odd and > 11. D is chosen by method "A*". f is 1 to indicate that n is a probable prime, 0 for composite and -1 if n is a square. Requires unit MultiP. Calls Add, Sub, Mul, Dvd, Dvdby2, Isqrt, Equal, Modd, Jacobi, and Spwr. } { 6 June 1990 } UBASIC: There is a subroutine Cgcd(A,B,&C,&Ca,&Cb). This is part of the unit MultiP in the Pascal versions. PROCEDURE Chinese(n:longint; VAR A,B,M:Chn; VAR soln, modulus : number); { Algorithm for Chinese remaindering of the system A[i] X = B[i] (mod M[i]), 1 = 1,..,n. It is NOT assumed that the moduli M[i] are pairwise relatively prime. Soln holds the solution. If there is no solution, modulus = 0, otherwise modulus is the modulus of uniqueness of the solution. Variables A, B, and M are input variables; they are VAR only for efficiency. The type Chn is imported from the unit MultiP. Requires MultiP. Calls Fcgcd, Dvd, Mul, Sub, Add, and Mod. Uses the type Chn. } { 4 August. 1990 } PROCEDURE Confrat(VAR C : tenr; leng : longint; VAR a,b : number); { From the finite continued fraction C this procedure calculates the rational a/b that corresponds. Reverses the procedure Ratconf. Parameter C is VAR only for efficiency; it is an input parameter. leng is the last index used in C. Note that it is assumed that leng >= 0. Requires unit MultiP, including type tenr. Calls Mul and Add. } { 7 Jan. 1990 } PROCEDURE Elcfctr(n,c1,c2,a : number; VAR f : number); { Elliptic curve method to factorize n, returning a factor in f. The parameters c1,c2, and a are chosen randomly - more than one choice may be necessary. Requires unit MultiP. Calls Add, Sub, Mul, Modd, and Gcd. } { 14 March 1991 } PROCEDURE Fermat(a : number; ubnd : longint; VAR f,g : number); { Fermat's method of factoring a. ubnd limits the number of iterations. The factors (or -1 if ubnd is reached) are placed in f and g. a must be positive and odd. Requires Unit MultiP. Calls Isqrt, Mul, Add, Sub, Equal, and Dvdby2. } { 30 Dec. 1989 } UBASIC: There is a subroutine Fcgcd(A,B,&C,&Ca). This is part of the unit MultiP in the Pascal versions. PROCEDURE Gpcftoq(VAR CF : tenr; leng1,leng2 : longint; VAR A,B,C : number); { General periodic continued fraction to quadratic routine. CF contains the continued fraction - indexes 0 through leng1 contain the first part before the period while indexes leng1+1 through leng2 contain the complete period. It is assumed that the continued fraction represents a positive number. CF is an input variable, made VAR for efficiency. A, B, and C return the coefficients of the quadratic satisfied by the continued fraction z : A z^2 + Bz + C = 0. There is no error checking to detect improper input. Requires unit MultiP, including type tenr. Calls Mul, Add, Sub, Gcd, and Dvd. } { 8 Jan. 1990 } PROCEDURE Lambda(VAR P, M : triala; k : longint; VAR Lam : number); { Returns in Lam the value of Carmichael's function - the minimum universal exponent. P, M, and k have the same meaning as in Sigma, and P and M are input varaibles - VAR only for efficiency. Defines by convention Lambda(1) = 0. Requires unit MultiP, including type triala. Calls Mul, Sub, Dvd, and Gcd. } { 2 May 1990 } PROCEDURE LDiosys(r,c : integer; VAR A,V : NumMat; VAR soln : Boolean; VAR rank : integer); { Solves the linear Diophantine system represented by the (augmented) matrix A. Matrix V contains a record of the manipulations. A is r by c+1, while V is c by c. Soln is TRUE if there is a solution, otherwise FALSE. Rank is the rank of the system. Contains the function Pivotind. Requires unit MultiP. Calls Sub, Mul, Dvd, Less, and Modd. Type NumMat is supplied by MultiP. } { 16 June 1990 } This is the algorithm described in "A COMPUTER LABORATORY MANUAL FOR NUMBER THEORY, by Donald E. G. Malm, COMPRess, 1980. PROCEDURE Lehmer(p,q : number; VAR r : number); { D. H. Lehmer's method of solving x^2 = q (mod p), assuming that p is an odd prime and q is a quadratic residue modulo p. The solution is returned in r. Requires unit MultiP. Calls Add, Sub, Mul, Dvd, Dvdby2, Spwr, Equal, Less, Jacobi, and Modd. } { 7 June 1990 } The algorithm is described in Lehmer's article "COMPUTER TECHNOLOGY APPLIED TO THE THEORY OF NUMBERS", in "STUDIES IN NUMBER THEORY", MAA 1969, pp. 117 - 151. FUNCTION Mobius(VAR M : triala; k : longint) : longint; { Returns the Mobius function. M and k have the same meaning as in Tau. Requires unit MultiP for the type triala. Calls no procedures nor functions.} { 22 April 1990 } UBASIC supplies this as a function in the language. ******************************** BODY2 ************************************ PROCEDURE Mppt(n,Q : number; VAR f : Boolean); { Implements Malm's probabilistic test for primality (or compositeness). Assumes n is odd, not a square, and that (Q|n)=1. Returns f = TRUE if n is probably prime, and f = FALSE if n is composite. Requires unit MultiP. Calls Add,Sub, Mul, Dvdby2, Modd, Gcd, Less, Jacobi, and Equal. } { 7 June 1990 } This procedure Mppt and the following one Mst are very effective unpublished probable prime tests devised by the author. PROCEDURE Mst(n : number; VAR f : Boolean); { Test for probable primality or compositeness devised by Malm. Returns f = TRUE if n is a probable prime, and f = FALSE if n is composite. Requires unit MultiP. Calls Add, Sub, Mul, Dvdby2, Modd, Isqrt, Gcd, Equal, Jacobi,and Less. Contains a subprocedure Lehmerc, which is the procedure Lehmer adapted to input that is not necessarily prime. } { 21 August 1990 } FUNCTION NxtPrm(a: integer): longint; { Returns the smallest prime greater than a, or else 0 if a < 0 or a >= 10,000. Returns an integer, not a number. Requires unit MultiP for the array prm[]. Calls no procedures nor functions. } { 20 April 1990 } UBASIC supplies this as a part of the language. PROCEDURE Pell(d : number; VAR x,y : number; VAR f : longint); { Uses continued fractions to compute the smallest nontrivial positive solution to x^2 - d y^2 = f, where f is 1 or -1. Note that f is computed by the program; it is not your choice. f = 0 indicates improper input. Requires unit MultiP. Calls Isqrt, Mul, Sub, Add, Dvd, and Equal. } { 10 January 1990 } Procedure Peralta(p,x : number; VAR s : number); { Implements Peralta's (second) probabilistic algorithm to compute s for which s^2 = x (mod p), assuming that p is an odd prime and that x is a quadratic residue modulo p. Requires unit MultiP. Calls Add, Sub, Mul, Dvd, Dvdby2, Modd, Lcs, and Spwr. } { 16 May 1990. } This algorithm is given in "A SIMPLE AND FAST PROBABILISTIC ALGORITHM FOR COMPUTING SQUARE ROOTS MODULO A PRIME NUMBER", by Rene C. Peralta, IEEE Trans. Info. Thry.,v. IT-32, pp. 846-848. PROCEDURE Perrin(n : number; VAR mc,mb,ma,a,b,c : number); { Given n, returns the signature of n (mod n) in the six output parameters. Perrin-Shanks-Adams sequence. } { 18 May 1990 } This algorithm is described in "STRONG PRIMALITY TESTS THAT ARE NOT SUFFICIENT", by William Adams and Daniel Shanks, Math. Comp. v. 39, 1982, pp. 255-300. PROCEDURE Polpm1(n,b : number; ubnd : longint; VAR f,g : number); { Pollard p-1 method of factoring n using base b. ubnd is the maximum number of iterations. f and g hold the factors (or -1 if the method fails or 0 if n < 2.). n and b should be relatively prime. In practice, b is small and n has been determined to have no small factors. Needs Unit MultiP. Calls Less, Spwr, Sub, Gcd, and Dvd. ubnd must be less than Base. } { 30 Dec. 1989 } The algorithm is described in "PRIME NUMBERS AND COMPUTER METHODS FOR FACTORIZATION", by Hans Riesel, Birkhauser 1985. It is also in the book by D. Bressoud listed above. PROCEDURE Primen1(n:number; VAR P:triala; leng:longint; VAR f:longint); { Attempts to prove n (assumed odd) is prime, given a partial prime factorization of n-1. The array P is an input parameter made VAR for efficiency, which holds the prime factors of n-1. Leng is the last index used in P[]. If n is prime then f= 1 is returned,while f= 0 is returned for no information and f= -1 for n composite. Requires unit MultiP. Calls Sub, Dvd, Spwr, Gcd, and Equal. Uses the array prm[]. } { 24 May 1990 } See Riesel or Bressoud. Function Prmdiv(a : number) : longint; { Returns the smallest prime factor of a, or else 0. Assumes that Base is 10^8. Finds the smallest prime factor of all numbers less than Base, and of those larger numbers that have a prime factor less than 10,000. Note that Prmdiv is always less than Base (which is 10^8). Requires unit MultiP, including array prm[]. Calls Less and Modd.} { 20 April 1990 } UBASIC supplies Prmdiv in the language. FUNCTION Prroot(n : number ):longint; { Returns the smallest positive primitive root of n. n > 0 is assumed to be prime and n-1 must be factorable by the function Prmdiv. Returns 0 for failure or error. Requires unit MultiP. Calls Prmdiv, Less, Equal, Sub, Dvd, Mul, Add, and Spwr. Note that Prmdiv is not in the unit MultiP. } { 2 May 1990 } PROCEDURE Quadtocf(m,d,q : number; VAR CF : tenr; VAR lengs,lengf : longint; VAR s : Boolean); { Calculates the continued fraction for the positive quadratic irrational (m + Ăd)/q and places it in the array CF. s = TRUE indicates the array was long enough to hold the continued fraction. s = FALSE indicates either that error or that d is a perfect square. lengs and lengf are indexes of the start and end respectively of the periodic part. Type tenr is imported from the unit. Calls Isqrt, Equal, Mul, Sub, Dvd, and Add. } { 6 May 1990 } PROCEDURE Ratconf(a,b : number; VAR C : tenr; VAR leng : longint); { Calculates the continued fraction of the rational a/b, terms are in C and the last index of C that is used is in leng. (leng = -1 indicates improper input.) Note that indexing starts with 0. It is assumed that b > 0. No error-checking is done to assure that the array isn't overrun, so range-checking should be on. The type tenr is imported from MultiP. Requires unit MultiP, including the type tenr. Calls Dvd. } { 7 Jan. 1990 } PROCEDURE Rennet(VAR CF : tenr; leng : longint; VAR A,B,C : number); { Given a purely periodic continued fraction in CF, with leng as the last index used, produces the coefficients A, B, and C of a quadratic that the continued fraction satisfies: Az^2 + Bz + ‚ =0. It is assumed that indices 0 through leng comprise the integer part plus the complete period of the continued fraction, and that the cont. fraction represents a positive number. There is no error checking for proper input. The parameter CF is VAR only for efficiency - it is an input parameter.If leng < 0 then the output is A=B=C=0. Requires unit MultiP, including type tenr. Calls Mul, Add, Gcd, Dvd, and Sub. } { 7 Jan. 1990 } ********************************** BODY3 ************************************ PROCEDURE Ressol(p,a : number; VAR r : number); { Shank's algorithm to solve x^2 = a (mod p). The solution is returned in r. It is assumed that p is an odd prime, and that a is a quadratic residue modulo p. Requires unit MUltiP. Calls Equal, Add, Sub, Mul, Dvdby2, Modd, Jacobi, and Spwr. } { 15 May 1990 } This algorithm is described in "FIVE NUMBER-THEORETIC ALGORITHMS", by Daniel Shanks, in PROC. SECOND MANITOBA CONFERENCE ON NUMERICAL MATH., 1972, pp. 51-70. PROCEDURE Rho(n,c : number; ubnd : longint; VAR f,g : number); { Brent's modification of Pollard's rho method of factoring n. c is the constant in the iteration formula x --> x*x + c. ubnd is the maximum number of iterations, while f and g hold the factors ( or -1 if factorization fails, or 0 if n <= 0). n must be positive. The user may want to change the constant cc which governs how often the gcd is calculated. Needs Unit MultiP. Calls Add, Mul, Modd, Sub, Gcd, and Dvd. } { 30 Dec. 1989 } See Riesel or Bressoud. PROCEDURE Sformf(n : number; VAR f : number); { The square form factoring method of Shanks. This implementation uses two continued fractions and rejects small denominators. The constant zzz is the upper limit for the number of iterations. n must be positive. Requires the unit MultiP. Calls Isqrt, Mul, Sub, Add, Dvd, Less, Equal, and Dvdby2. } { 6 May 1990 } This factoring method is described in Riesel. PROCEDURE Sigma(VAR P,M : triala; k : longint; VAR sig : number); { For the number n (which is not a parameter) its prime factor and multiplicity arrays and the last index used in them, namely k, are input parameters. The value of Sigma (the sum of the divisors of n) is computed using the standard formula, and returned in sig. P and M are input variables which are VAR only for efficiency. Requires unit MultiP. Calls Mul and Add.} { 2 May 1990 } FUNCTION Spspt(n,b : number) : Boolean; { Strong pseudoprime test of n using witness b. n >= 2 is required. Needs Unit MultiP. Calls Dvdby2, Spwr, Lmod, Sub, Less, Equal, and Mul. } { 20 May 1990 } This is described in the article "THE PSEUDOPRIMES TO 25*10^9" by Carl Pomerance, J. L. Selfridge, and Samuel S. Wagstaff, Jr., Math. of Comp. v.35 (1980) pp1003-1026. PROCEDURE Squf(n : number; VAR f : number); { Implements the square forms factorization method of Shanks, using the numerators of the convergents instead of a second continued fraction. These numerators are subject to overflow. n must be positive. Requires the unit MultiP. Calls Isqrt, Mul, Sub, Add, Dvd, Equal, and Gcd. } { 6 May 1990 } This is probably not as fast as the other implementation (Sformf). PROCEDURE Tau(VAR M : triala; k : longint; VAR ta : number); { For the number n (which is not a parameter) its prime factor multiplicity array M and the last index used in M, namely k, are input parameters. The value of Tau (the number of divisors of n) is computed using the standard formula. M is VAR only for efficiency. Requires unit MultiP. Calls Mul. } { 2 May 1990 } PROCEDURE Tdvd(a: number; VAR P,M: triala; VAR k: longint; VAR b: number); { Trial-divides a, using the primes in the prime list. The prime factors and their multiplicities are stored in P[] and M[] respectively. The variable b holds the remaining unfactored part (or 0 if a = 0). k is the last index used in the arrays P and M. Note that a is allowed to be negative. Requires unit MultiP, including the array prm and type triala. Calls Dvdby2 and Dvd. Note that a prime factor found by Tdvd must be less than Base (10^8) (Pascal versions). } { 20 April 1990 } Pascal version: If a number has all prime divisors less than Base, and if all but one of them are less than 10^5, then Tdvd will completely factor it. UBASIC version: If a number has all prime divisors less than 2^34, and if all but one are less than 2^17, then Tdvd will completely factor it. PROCEDURE Tenner(n : number; VAR A : tenr; VAR leng:longint; VAR f:Boolean); { Tenner's algorithm for the continued fraction of Ăn. A contains the integer part (in A[0]) and a full period. leng is the last index used. Improper input is indicated by setting leng = -1. If the period is too long to fit into A, then f = FALSE, else f = TRUE. The type tenr is imported from the unit MultiP. The constant ub = 100, making it easy to change the text of this procedure to accomodate a larger array. Requires MultiP, including the type tenr. Calls Isqrt, Mul, Sub, Add, Dvd, and Equal. } { 3 Jan. 1990 } PROCEDURE Tenners(n:number; VAR A:tenr; VAR leng:longint; VAR f,odd:Boolean); { Tenner's algorithm for the continued fraction of Ăn. A contains the integer part (in A[0]) and a half period. leng is the last index used. Improper input is indicated by setting leng = -1. If the period is too long to fit into A, then f = FALSE, else f = TRUE. odd = TRUE if the symmetric part of the period is of odd length, else odd = FALSE and the symmetric part is of even length. The type tenr is imported from the unit MultiP. The constant ub = 100, making it easy to change the text of this procedure to accomodate a larger array. Requires MultiP, including the type tenr. Calls Isqrt, Mul, Sub, Add, Dvd, and Equal. } { 4 Jan. 1990 } Since Tenners stores only half as much as Tenner, it doesn't need as long an array. PROCEDURE TestPrime(n : number; VAR ans : longint); { This primality test is valid for n less than 50 (U.S.) billion. TestPrime first uses a base-2 strong pseudoprime test, then the Shanks-Adams-Perrin sequence test, and then checks the six composites less than 5*10^10 which pass both tests. Small numbers are tested separately by looking in prm[]. Returns 1 for prime, 0 for composite (or 1), and -1 for error: too big or 0. Requires unit MultiP. Calls Sub, Less, Dvdby2, Spwr, Equal, Mul, Lmod, Modd, and Add. } { 20 May 1990 } See "FAST PRIMALITY TESTS FOR NUMBERS LESS THAN 50*10^9", by G. C. Kurtz, Daniel Shanks, and H. C. Williams, Math. of Comp. v.46 (1986) pp.691-701. PROCEDURE Totient(VAR P,M : triala; k : longint; VAR tot : number); { Returns the Euler totient function in the variable tot. P, M, and k have the same meaning as in Sigma. Requires unit MultiP, including type triala. Calls Mul and Sub. } { 2 May 1990 } This is supplied in the language UBASIC as a function EUL(N). PROCEDURE TwoSum(p : number; VAR a, b : number); { Expresses p = a^2 + b^2 for prime p = 1 (mod 4). Note that primality is not checked. Requires the unit MultiP. Calls Modd, Equal, Add, Sub, Spwr, Dvd, Dvdby2, Less, Jacobi, and Isqrt. } { 2 May 1990 } See "NOTE ON REPRESENTING A PRIME AS A SUM OF TWO SQUARES", by John Brillhart, Math. of Comp. v.26 (1972) pp. 1011-1013. PROCEDURE Wilpp1(n,P : number; ub : longint; VAR f : number); { William's p+1 method for factoring n, returning the result in f. If the method fails, f = 1 (or n). ub is the maximum number of cycles. The constant cc determines how often the gcd is calculated. Requires unit MultiP. Calls Sub, Gcd, Mul, Modd, and Equal. } { 7 June 1990 } This algorithm is described in Bressoud.