Big-integer arithmetic in BASIC

→Library header file
→Number theory modules survey
→Download & links

Motivation

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.

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
'
' ****************************************************************************************

Number theory modules survey

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 a_i·x ≡ b_i mod M_i.

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 = f_red(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.

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:

G.H. Hardy and E.M. Wright, An Introduction to the Theory of Numbers, Oxford, 1978.
H. Cohen, A Course in Computational Algebraic Number Theory, Berlin, 1993.
A. Menezes, P. van Oorschot and S. Vanstone, Handbook of Applied Cryptography, Boca Raton, 1996.
A History of Algorithms - from the Pebble to the Microchip, J.-L. Chabert, ed., Berlin, 1998.
F. Beukers, Getaltheorie voor beginners, Utrecht, 1999.
Plus many useful papers and abstracts gathered from all around the World Wide Web.

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.