Power-Programmierung/CD1.mdf/basic/compiler/ubasic/malm/malm.doc

Text File | 1989-04-25 | 19.7 KB | 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.