home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Power-Programmierung
/
CD1.mdf
/
basic
/
compiler
/
ubasic
/
malm
/
malm.doc
< prev
next >
Wrap
Text File
|
1989-04-25
|
20KB
|
428 lines
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.