Big-integer arithmetic in BASIC

down Motivation

Euclid, litho by Fiorini
Euclid (flourished ca. 300 B.C.)

The material on this page is due to my interest in the theory of numbers, a field whence originate many classic algorithms you can't wait to implement and watch them work their magic. Alas, as soon as things get captivating, you run unto the limits the hardware imposes:

‘[..] there is a bound to the number of symbols which the computer can observe at one moment. If he wishes to observe more, he must use successive operations.’ Alan Turing, On computable numbers [..], 1936.

So I needed a little library to rephrase integer arithmetic in terms of ‘successive operations’. Instead of using an existing BigInt library, I preferred writing my own to obtain a minimum sized-set of concisely implemented operations called by the LargeInt modules listed below.

The project grew under its own impetus, eventually including some rather useful functions. There are versions for FreeBasic, QuickBasic, XBasic and Visual Basic. Please visit the Download page for details.

down Library header file

' *************************************************************************
'Subject: Include file for large-integer arithmetic Basic library.
'Author : Sjoerd.J.Schaper
'Update : 06-06-2009
'Code   : FreeBasic 1.06.0
'
#INCLIB "largeint"
' ****************************************************************************************
'                            Initialization and assignment
' ****************************************************************************************
'The LargeInt class is initialized with a 'LargeInit(size%, name$)' call at
'module level. You have to declare (constant) indices for referencing the
'large integers you're going to use (CONST p = 0, q = 1, ..), so that each
'one is associated with its own distinct pointer. These pointers are passed
'to the procedures. Don't call with same pointers ('Squ p, p'), and beware
'of duplicate references, for they will crash your program.
'
'In the following prototypes,
'letters p, q, r, d, m denote 16-bit (pseudo) pointers to largeint's;
'variables a, k, sw are 16-bit integers, |a| stands for absolute value,
'c is a 32-bit integer; g is a number string, f and h are any string.
'The t#'s within brackets denote the range of scratchpad-registers shared
'by the subroutine, that cannot be passed as arguments.
'
DECLARE FUNCTION LargeInit (ByVal k AS SHORT, ByVal f AS ZSTRING PTR) AS SHORT
'Initialize the number-arrays, open the primetable and logfile.
'To be called with k = number of largeint's in your program and
'string f = the name of the program. If f = "" no logfile is opened.
'Return value is the maximum largeint length measured in array elements ('words'),
'the decimal number length is about 4.51 * largeint length. If overflow occurs,
'then increase number array-size constant Asize in largeint.bas and recompile.
'
DECLARE SUB Letf (ByVal p AS SHORT, ByVal c AS LONG)
'Assign signed long integer c to large integer p.
'
DECLARE SUB Readst (ByVal p AS SHORT, ByVal g AS ZSTRING PTR)
'Assign signed (hexa)decimal number string g to largeint p.
'There's no check on non-digit characters in the string,
'add prefix &H to indicate a hex number string. {t1..t2}
'
DECLARE SUB Rndf (ByVal p AS SHORT, ByVal k AS SHORT)
'Assign a random positive value of bitlength k to largeint p.
'
DECLARE FUNCTION LibErr () AS INTEGER
'Return run-time error code and reset.
'
DECLARE SUB Term ()
'Close the primetable and logfile.
'
' ****************************************************************************************
'                                 Conversion functions
' ****************************************************************************************
DECLARE FUNCTION Logf (ByVal p AS SHORT) AS DOUBLE
'Return double precision natural logarithm of |largeint p|.
'
DECLARE FUNCTION Bufl (ByVal p AS SHORT) AS LONG
'Return the maximum stringlength of decimal p, call before CnvSt.
'
DECLARE SUB CnvSt (ByVal g AS ZSTRING PTR, ByVal p AS SHORT)
'Convert largeint p to a string of decimal digits. First create
'stringbuffer g = STRING(Bufl(p), "0") and pass to CnvSt. {t2}
'
DECLARE SUB RatCnv (ByVal g AS ZSTRING PTR, ByVal p AS SHORT, ByVal q AS SHORT, _
 ByVal a AS SHORT)
'Convert rational p / q to a k-long string of base a-digits, 2 <= a <= 17.
'Pass stringbuffer g = STRING(k, "0"), 9 < k < 32768. {t0..t3}
'
#IFNDEF __VB
' ****************************************************************************************
'                                   Obtaining output
' ****************************************************************************************
DECLARE SUB Printn (ByVal p AS SHORT, ByRef f AS STRING, ByRef h AS STRING, _
 ByVal sw AS SHORT)
'Print largeint p with prefix f and suffix h to screen and the logfile.
'Set switch = 1 to add a CrLf, switch = 2 to also attach length(p). {t2}
'
DECLARE SUB Printr (ByVal p AS SHORT, ByVal q AS SHORT, ByRef f AS STRING, _
 ByVal k AS SHORT, ByVal sw AS SHORT)
'Print the decimal expansion of rational p / q with prefix f and length k.
'Set switch = 1 to add a CrLf. {t0..t3}
'
DECLARE SUB Prints (ByRef f AS STRING, ByVal sw AS SHORT)
'Print string f. Set switch = 1 to add a CrLf, switch = 2 for two.
'
#ENDIF
DECLARE SUB Work (ByRef f AS STRING)
'Set string f to the working directory in ENVIRON("LARGEINT"),
'the stringbuffer must be large enough to hold the path.
'
' ****************************************************************************************
'                              Basic arithmetic functions
' ****************************************************************************************
DECLARE SUB Subt (ByVal p AS SHORT, ByVal q AS SHORT)
'Set large integer p to p - q.
'
DECLARE SUB Add (ByVal p AS SHORT, ByVal q AS SHORT)
'Set large integer p to p + q.
'
DECLARE SUB Inc (ByVal p AS SHORT, ByVal a AS SHORT)
'Increment largeint p by signed one-word value a.
'
DECLARE SUB Divd (ByVal p AS SHORT, ByVal d AS SHORT, ByVal q AS SHORT)
'Divide p by d, set p to the remainder and q to the quotient. {t2}
'
DECLARE FUNCTION Divint (ByVal p AS SHORT, ByVal a AS SHORT) AS SHORT
'Divide largeint p by |one-word a|,
'return |remainder| and set p to the quotient.
'
DECLARE SUB Ldiv (ByVal p AS SHORT, ByVal d AS SHORT)
'Set p to quotient p / d. {t1..t2}
'
DECLARE SUB Mult (ByVal p AS SHORT, ByVal q AS SHORT, ByVal r AS SHORT)
'Set large integer r to p * q.
'
DECLARE SUB Lmul (ByVal p AS SHORT, ByVal q AS SHORT)
'Set p to product p * q. {t2}
'
DECLARE SUB Squ (ByVal p AS SHORT, ByVal q AS SHORT)
'Set q to the square of p, faster than p * p.
'
DECLARE SUB Powr (ByVal p AS SHORT, ByVal c AS LONG)
'Raise largeint p to the power |longword c|. {t0..t2}
'
DECLARE SUB Chs (ByVal p AS SHORT)
'Change the sign of largeint p.
'
DECLARE FUNCTION Absf (ByVal p AS SHORT) AS SHORT
'Clear sign bit of largeint p and return the current value.
'
' ****************************************************************************************
'                                 Copying and comparison
' ****************************************************************************************
DECLARE SUB Dup (ByVal p AS SHORT, ByVal q AS SHORT)
'Copy large integer p to q.
'
DECLARE SUB Swp (ByVal p AS SHORT, ByVal q AS SHORT)
'Exchange largeint's p and q.
'
DECLARE FUNCTION Cmp (ByVal p AS SHORT, ByVal q AS SHORT) AS SHORT
'Compare returns -1 if p < q, 0 if p = q, or 1 if p > q.
'
DECLARE FUNCTION Isf (ByVal p AS SHORT, ByVal a AS SHORT) AS SHORT
'Boolean: check if largeint p equals one-word value a.
'
DECLARE FUNCTION Ismin1 (ByVal p AS SHORT, ByVal m AS SHORT) AS SHORT
'Boolean: check if p = m - 1
'
' ****************************************************************************************
'                                   Bit manipulation
' ****************************************************************************************
DECLARE SUB Boolf (ByVal p AS SHORT, ByVal q AS SHORT, ByVal k AS SHORT)
'Bitwise Boolean functions set largeint p to q Op(k) p, with Op(1) = AND,
'Op(2) = XOR, Op(3) = OR. If k = 0 then p is set to NOT p, and q is ignored.
'In each case, the sign of p is left unchanged.
'
DECLARE SUB Lsft (ByVal p AS SHORT, ByVal k AS SHORT)
'Left-shift largeint p by k bits; k < 0 shifts |k| full words.
'
DECLARE SUB Rsft (ByVal p AS SHORT, ByVal k AS SHORT)
'Right-shift largeint p by k bits; k < 0 shifts |k| full words.
'
DECLARE FUNCTION Odd (ByVal p AS SHORT) AS SHORT
'Remove trailing zero bits from largeint p, return (right) shift count.
'
DECLARE FUNCTION Bitl (ByVal p AS SHORT) AS LONG
'Return the bitlength of largeint p.
'
' ****************************************************************************************
'                             Modular arithmetic functions
' ****************************************************************************************
DECLARE SUB Moddiv (ByVal p AS SHORT, ByVal m AS SHORT)
'Compute positive residue p modulo largeint m. {t1..t2}
'
DECLARE SUB Modbal (ByVal p AS SHORT, ByVal m AS SHORT)
'Balanced residue p modulo m: |p| <= m / 2. {t1..t2}
'
DECLARE FUNCTION Mp2 (ByVal p AS SHORT, ByVal a AS SHORT) AS SHORT
'Return largeint p modulo a, which must be a one-word power of 2.
'
DECLARE SUB Modmlt (ByVal p AS SHORT, ByVal q AS SHORT, ByVal m AS SHORT)
'Set largeint p to p * q (modulo m). {t1..t2}
'
DECLARE SUB Modsqu (ByVal p AS SHORT, ByVal m AS SHORT)
'Set largeint p to p ^ 2 (modulo m). {t1..t2}
'
DECLARE SUB Modpwr (ByVal p AS SHORT, ByVal q AS SHORT, ByVal m AS SHORT)
'Set largeint p to base p ^ exponent q (modulo m). {t0..t2, t3 if q < 0}
'
' ****************************************************************************************
'                             Number theoretic functions
' ****************************************************************************************
DECLARE SUB Fctl (ByVal p AS SHORT, ByVal a AS SHORT)
'Set largeint p to factorial a(a-1)···2·1. {t1..t2}
'
DECLARE SUB Lcm (ByVal p AS SHORT, ByVal q AS SHORT, ByVal d AS SHORT)
'Set largeint d to the least common multiple (p, q). {t0..t2}
'
DECLARE SUB Gcd (ByVal p AS SHORT, ByVal q AS SHORT, ByVal d AS SHORT)
'Set largeint d to the greatest common divisor (p, q). {t0..t2}
'
DECLARE SUB Euclid (ByVal p AS SHORT, ByVal q AS SHORT, ByVal d AS SHORT)
'Apply the extended Euclidean algorithm to largeint's p and q,
'set p to 1/p (modulo q), q to q/gcd(p, q), and d to gcd(p, q). {t0..t2}
'
DECLARE SUB Bezout (ByVal p AS SHORT, ByVal q AS SHORT, ByVal d AS SHORT)
'Solution of the Diophantine equation px + qy = gcd(p, q),
'set largeint p to x, q to y, and d to gcd(p, q). {t0..t3}
'
DECLARE SUB Isqrt (ByVal p AS SHORT, ByVal q AS SHORT, ByVal r AS SHORT)
'Set r to the truncated integer part of square root(p),
'and q to residue p - r ^ 2. {t0..t2}
'
DECLARE FUNCTION IsSqr (ByVal p AS SHORT, ByVal r AS SHORT) AS SHORT
'Boolean: return -1 if largeint p is a perfect square, else zero.
'If true, set r to the square root of p. {t0..t3}
'
DECLARE FUNCTION Kronec (ByVal p AS SHORT, ByVal q AS SHORT) AS SHORT
'Return Kronecker's quadratic residuosity symbol (p/q) = 0, 1, or -1. {t0..t3}
'
DECLARE FUNCTION IsPPrm (ByVal p AS SHORT, ByVal d AS SHORT, ByVal k AS SHORT) AS SHORT
'Check primality of p with the Miller-Rabin test. Return 1 if
'p is definitely prime, -1 if probably prime, 0 if p is composite.
'k is the number of witnesses, k < 0 skips initial trial divisions. {t0..t3}
'
DECLARE FUNCTION Nxtprm (ByRef sw AS SHORT) AS LONG
'Loop through primetable PrimFlgs.bin. Initialize with sw = 0,
'then odd(sw) returns the next prime with each successive call.
'Use sw = 2 to push the current state and reset, sw =-2 to pop.
'
DECLARE FUNCTION PrimCeil () AS LONG
'Return the upper limit of primetable PrimFlgs.bin.
'
DECLARE SUB Triald (ByVal q AS SHORT, ByVal p AS SHORT, ByVal c AS LONG)
'Trial divide p up to limit c. Primedivisors w/multiplicities are
'stored in list {q}, largeint p is set to the remaining cofactor. {t0..t3}
'
' ****************************************************************************************
'                                    Direct access
' ****************************************************************************************
DECLARE SUB EnQ (ByVal q AS SHORT, ByVal p AS SHORT, ByVal a AS SHORT)
'Store largeint p and |one-word a| in sequential list {q}.
'As {q} contains negative record separators, an attempt to
'simply 'Printn q, ...' will crash your program. Instead use:
'
DECLARE FUNCTION ExQ (ByVal p AS SHORT, ByRef a AS SHORT, ByVal q AS SHORT, _
 ByRef k AS SHORT) AS SHORT
'Get successive numbers p and |a| from list {q} at offset k.
'Set k = 0 to initialize reading, returns zero if @end-of-list.
'
DECLARE FUNCTION Getl (ByVal p AS SHORT) AS SHORT
'Return the length of largeint p measured in words.
'
DECLARE FUNCTION Gets (ByVal p AS SHORT) AS SHORT
'Return the sign bit of largeint p: -32768 if p < 0, else zero.
'
DECLARE FUNCTION Getw (ByVal p AS SHORT, ByVal k AS SHORT) AS SHORT
'Return word number k of largeint p: a base MB digit (see below).
'
DECLARE SUB Setl (ByVal p AS SHORT, ByVal a AS SHORT)
'Set the word length of largeint p to a.
'
DECLARE SUB Sets (ByVal p AS SHORT, ByVal a AS SHORT)
'Set the sign bit of largeint p if a < 0, else clear.
'
DECLARE SUB Setw (ByVal p AS SHORT, ByVal k AS SHORT, ByVal a AS SHORT)
'Set word number k of largeint p to |a| < MB.
'
' ****************************************************************************************
'                                  Internal functions
' ****************************************************************************************
DECLARE SUB Lftj (ByVal p AS SHORT, ByVal k AS SHORT)
'Left-adjust largeint p, beginning at word number k.
'
DECLARE FUNCTION Hp2 (ByVal p AS SHORT) AS SHORT
'Return 2^(highest set bit in the leftmost word of largeint p).
'
' ****************************************************************************************
'                                   Global constants
' ****************************************************************************************
CONST LB = 15, MB = &H8000& ' binary storage base
CONST LD = 4, MD = 10000 '    decimal string base < MB
CONST t0 = -1, t1 = -2 '      pointers to swap-registers
CONST t2 = -3, t3 = -4
'
' ****************************************************************************************

down Number theory modules survey

Pierre Fermat, anonymous engraving
Pierre Fermat (1601-1665)

With the library come a few BASIC-modules pertinent to number theory, which may serve as a test suite. Running these modules with the sample input files will demonstrate correctness and performance of the library.

Sequences and series

Fibonacc.bas
Fibonacci numbers and the golden ratio phi.
A simple recurrent sequence, minimal module sample.
Pi.bas
Calculation of pi with Gauss's equation π/4 = 12·arctan(1/18) + 8·arctan(1/57) - 5·arctan(1/239).
Transcnd.bas
Complex transcendental functions using rational Taylor series, includes a basic RPN-parser.

Euclid's algorithm

ExtEucli.bas
Minimal solution of the Diophantine equation ax − by = c,
also known as Bézout's identity.
Rational.bas
The continued fraction algorithm.
Simplify ratios and express decimals as common fractions.
FundUnit.bas
Fundamental unit (p + q·√d)/ 2 of the real quadratic field K = Z(√d).
Also given are the continued fraction and decimal approximation of the square root.

Congruences and residues

QuadCngr.bas
Roots of a quadratic congruence with prime modulus N: the Shanks-Tonelli algorithm Ressol.
For r² ≡−1 mod N the Diophantine equation x² + y² = N is solved with Cornacchia's algorithm.
ChineseR.bas
Solution of a system of linear congruences.
Apply the Chinese remainder theorem to simultaneous congruences ai·x ≡ bi mod Mi.
PrimRoot.bas
Test if base a belongs to phi(N) modulo N,
input a = 1 to find the least primitive root of N.
ZnStruct.bas
Structure of the multiplicative group of units modulo N:
another Smith normal form application.
DiscrLog.bas
Silver-Pohlig-Hellman method for discrete logarithms, using Pollard rho.
Genrator.bas
Cryptography: probable primes and crypto key generation.
Prints relevant values to public and private key-files.
Encoder.bas
Exponential encryption with salt (RSA).
Decoder.bas
Decryption using the Chinese remainder theorem.

Primes and primality

PrimFlgs.bas
The Sieve of Erathostenes: a variation of Chartres' Algorithm 311.
Generates an encoded list of prime numbers read by library function Nxtprm().
Goldbach.bas
The strong Goldbach conjecture: every even integer > 2 is the sum of two primes.
NxtPrime.bas
Find the next probable prime ≥ N with the Miller-Rabin test.
PowrModC.bas
Gaussian integer exponentiation: complex base and modulus, rational exponent.
The Gaussian Miller-Rabin test is performed for exponent ‘Z-1’.
WilsonTh.bas
Test for primality with Wilson's theorem: iff p is prime (p − 1)! ≡−1 mod p.
Mersenne.bas
Lucas-Lehmer test for Mersenne numbers 2^p − 1, with odd prime p.
Toy version for small p.
LucasFun.bas
Lucas pseudoprime test, quick calculation of Fibonacci giants.
PowrMtrx.bas
Pseudoprime test for N using k-th order recurrence relations,
powering companion matrix M for characteristic polynomial f(x).

Integer factorization

TrialDiv.bas
Prime factorization and arithmetic functions.
Calculates Euler totient and Carmichael lambda functions,
also divisor functions μ(a), ω(a), Ω(a), δ(a) and σ(a).
GaussInt.bas
An extension of the previous program, to work with Gaussian integers.
PollaRho.bas
Pollard's Rho Monte Carlo factorization method, Brent's modification.
PollarP1.bas
Pollard's double-step P−1 method finds factor p of N
if p−1 is a product of small primes and a single larger prime.
WilliaP1.bas
Williams's double-step, Lucas series P+1 method finds factor p of N
if all prime factors but one of p±1 are at most B, and a single larger one is at most B1.
EllCrvFr.bas
Lenstra's elliptic curve method (ECM) finds factor p of N
if some random elliptic curve group has a smooth order mod p.
FermatFr.bas
Factorization by the difference of two squares.
Fermat's method is useful only if N has factors close to its square root.
CFracFac.bas
The 1975 Brillhart-Morrison continued-fraction factorization method (CFrac)
with multipliers, early abort strategy and large prime relations.
SquFoFac.bas
Shanks's square form factorization (SquFoF) with queue and multipliers:
non-trivial ambiguous quadratic forms of discriminant 4kN yield factors of N.
Pmpqs.bas
Pomerance's quadratic sieve factoring with multiple polynomials and large prime relations.

Polynomials and roots

CnFcRoot.bas
Approximate real root of an integral poly with Lagrange's 1769 continued fractions method.
PolRoot.bas
The Newton-Raphson method.
Approximates any simple root of an integral poly given the right initial value.
RadpAdic.bas
p-Adic roots of an integral polynomial using Hensel's lemma.
BernPoly.bas
Sums of integer powers with Bernoulli polynomials.
PolDisc.bas
Discriminant of an integral polynomial via resultant
= determinant of Sylvester's matrix for f(x) and its derivative f'(x).
QuadForm.bas
Integral binary quadratic form f =(a,±b,c) of discriminant D.
The reduced forms and two representations a = fred(x,y) are also given.
Reprsent.bas
Solutions of the Diophantine equation ax² + bxy + cy² = N.
PowrQfb.bas
Powering binary quadratic forms: Gaussian composition with Shanks's Nucomp algorithm.
Input exponent k = 0 to compute class numbers with Pollard-Teske discrete log methods.

Linear systems and lattices

GaussEli.bas
Gauss-Bareiss fraction free elimination.
Bring matrix A into upper triangular form and solve A·X = B
LDioSyst.bas
Solution of a linear Diophantine system A·x = b using Hermite and Smith normal forms.
LLL_Hnf.bas
Alternate version of the previous program,
using Hnf via integral LLL-lattice basis reduction for shorter solutions.
LLL_int.bas
Find the minimal polynomial of a complex algebraic number,
or LLL-reduce the Gram matrix of a set of linearly independent basis vectors.
C.F. Gauss, litho by S. Bendixen (1828)
Carl Friedrich Gauss (1777-1855)

Many of these functions are now included in my VB large integer/Bigint RPN calculator.
 
Download

View the ReadMe or Changelog.

Mirror site (auf deutsch)

Computing Ackermann's function in Basic.

This site has an annex: Finite field arithmetic in BASIC.

↑Top   Node
 

References:

Disclaimer:

This BASIC code for large integers (big numbers, grote gehele getallen, grands nombres entiers, grandi numeri interi, números enteros grandes, inteiros grandes, μεγάλοι ακέραιοι αριθμοί, больших целых чисел, الأعداد الصحيحة الكبيرة &c.) is distributed free of charge in the hope that it will be useful, but without warranty of any kind, even the implied warranty of merchantability or fitness for a particular purpose. Permission is granted for modifying your copy of the code, or forming a non-commerical application based on it, and copy and distribute such modifications, provided that the modified content carries a reference to the originator of the source code and notices stating that it has been changed and by whom.

Valid xhtml 1.0! css valide! SourceForge.net logo Copyright © december 2003 by Sjoerd J. Schaper