Big-integer arithmetic in BASIC

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

Motivation

Euclid (flourished ca. 300 B.C.)

The material presented on this page is due to my interest in number theory. It's 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 basic arithmetics in terms of such ‘successive operations’. Instead of adapting one of the many existing BigInt or BigNum libraries for my purpose, it seemed satisfactory to start from scratch and implement the most necessary operations only.

The project grew under its own impetus, eventually including some very useful functions, but retaining the initial, simple lay-out. The library is written in pure QuickBasic (see Download page for FreeBasic, XBasic and Visual Basic versions), supplemented with assembly shift-routines.

Library header file

' *************************************************************************
'Subject: Include file for large-integer arithmetic Basic library.
'Author : Sjoerd.J.Schaper
'Date   : 02-02-2004
'Code   : QuickBasic, PDS, VBdos
'
' ****************************************************************************
'                       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.
'
DEFINT A-B, D-Z: DEFLNG C
DECLARE FUNCTION LargeInit% (k, f AS STRING)
'Initialize the number-arrays, open the primetable and logfile. To be called
'with k = number of largeint's in your program (k<=450 in the $STATIC version)
'and string f = the name of the program. If f = "" no logfile is opened.
'Return value is the maximal largeint length measured in array elements
'(='words'). The decimal number length is about 4.51 * largeint length.
'
DECLARE SUB Letf (p, c)
'Assign signed long integer c to large integer p.
'
DECLARE SUB Readst (p, g AS STRING)
'Assign decimal number string g to largeint p. Except for a prefixed minus
'sign, there's no check on non-digit characters in the string. {t1..t2}
'
DECLARE SUB Rndf (p, k)
'Assign a random positive value of bitlength k to large integer p.
'
DECLARE SUB Term ()
'Close the primetable and logfile.
'
' ****************************************************************************
'                            Conversion functions
' ****************************************************************************
DECLARE FUNCTION Decf% (p)
'Convert largeint p to base MD (see below), return the bytelength of
'decimal(p), call before CnvSt. {t1..t2}
'
DECLARE SUB CnvSt (g AS STRING)
'Convert base MD number to string. First call k = Decf(p),
'create stringbuffer g = STRING$(k, "0") and pass to CnvSt. {t2}
'
DECLARE FUNCTION Logf# (p)
'Return double precision natural logarithm of |largeint p|.
'
DECLARE SUB Ratdec (g AS STRING, p, q)
'Convert the decimal expansion of rational p/q in [0, 1) to string. {t0..t2}
'
' ****************************************************************************
'                              Obtaining output
' ****************************************************************************
DECLARE SUB Printf (f AS STRING, g AS STRING, h AS STRING, sw)
'Print prefix f, number string g and suffix h to screen and the logfile,
'set switch = 1 to attach a CrLf, switch = 2 to also attach length(g).
'
DECLARE SUB Printn (p, f AS STRING, h AS STRING, sw)
'Print largeint p with prefix f and suffix h using the Printf switches.
'{t1..t2}
'
DECLARE SUB Printr (p, q, f AS STRING, k)
'Print the decimal expansion of rational p/q with prefix f and length k.
'{t1..t2}
'
DECLARE SUB Prints (f AS STRING, sw)
'Print string f. Set switch = 1 to attach a CrLf, switch = 2 for two.
'
DECLARE SUB Work (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 (p, q)
'Set large integer p to p - q.
'
DECLARE SUB Add (p, q)
'Set large integer p to p + q.
'
DECLARE SUB Inc (p, a)
'Increment largeint p by signed one-word value a.
'
DECLARE SUB Divd (p, d, q)
'Divide p by d, set p to the remainder and q to the quotient. {t2}
'
DECLARE FUNCTION Divint% (p, a)
'Divide largeint p by |one-word a|,
'return |remainder| and set p to the quotient.
'
DECLARE SUB Ldiv (p, d)
'Set p to quotient p / d. {t1..t2}
'
DECLARE SUB Mult (p, q, r)
'Set large integer r to p * q.
'
DECLARE SUB Lmul (p, q)
'Set p to product p * q. {t2}
'
DECLARE SUB Squ (p, q)
'Set q to the square of p, faster than p * p.
'
DECLARE SUB Isqrt (p, q)
'Set q to the truncated integer part of square root(p). {t0..t2}
'
DECLARE SUB Chs (p)
'Change the sign of largeint p.
'
DECLARE FUNCTION Absf% (p)
'Clear sign bit of largeint p and return the current value.
'
' ****************************************************************************
'                            Copying and comparison
' ****************************************************************************
DECLARE SUB Dup (p, q)
'Copy largeint p to q.
'
DECLARE SUB Swp (p, q)
'Exchange largeint's p and q.
'
DECLARE FUNCTION Cmp% (p, q)
'Compare returns -1 if p < q, 0 if p = q, or 1 if p > q.
'
DECLARE FUNCTION Isf% (p, a)
'Boolean: check if largeint p equals one-word value a.
'
DECLARE FUNCTION Ismin1% (p, m)
'Boolean: check if p = m - 1
'
' ****************************************************************************
'                              Bit manipulation
' ****************************************************************************
DECLARE SUB Boolf (p, q, k)
'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 unaffected.
'
DECLARE SUB Lsft (p, k)
'Left-shift largeint p by k bits; k < 0 shifts |k| full words.
'
DECLARE SUB Rsft (p, k)
'Right-shift largeint p by k bits; k < 0 shifts |k| full words.
'
DECLARE FUNCTION Odd% (p)
'Remove trailing zero bits from largeint p, return (right) shift count.
'
DECLARE FUNCTION Bitl% (p)
'Return the bitlength of largeint p.
'
' ****************************************************************************
'                        Modular arithmetic functions
' ****************************************************************************
DECLARE SUB Moddiv (p, m)
'Compute positive residue p modulo largeint m. {t1..t2}
'
DECLARE SUB Modbal (p, m)
'Balanced residue p modulo m: |p| <= m / 2. {t1..t2}
'
DECLARE FUNCTION Mp2% (p, a)
'Return largeint p modulo a, which must be a one-word power of 2.
'
DECLARE SUB Modmlt (p, q, m)
'Set largeint p to p * q (modulo m). {t1..t2}
'
DECLARE SUB Modsqu (p, m)
'Set largeint p to p ^ 2 (modulo m). {t1..t2}
'
DECLARE SUB Modpwr (p, q, m)
'Set largeint p to base p ^ exponent q (modulo m). {t0..t2, t3 if q < 0}
'
' ****************************************************************************
'                         Number theoretic functions
' ****************************************************************************
DECLARE SUB Powr (p, k)
'Raise largeint p to the power |one-word k|. {t0..t3}
'
DECLARE SUB Fctl (p, a)
'Set largeint p to factorial a(a-1)···2·1. {t1..t2}
'
DECLARE SUB Lcm (p, q, d)
'Set largeint d to the least common multiple (p, q). {t0..t2}
'
DECLARE SUB Gcd (p, q, d)
'Set largeint d to the greatest common divisor (p, q). {t0..t2}
'
DECLARE SUB Euclid (p, q, d)
'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 (p, q, d)
'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 FUNCTION Kronec% (p, q)
'Return Kronecker's quadratic residuosity symbol (p/q) = 0, 1, or -1. {t0..t3}
'
DECLARE FUNCTION IsSqr% (p, q)
'Boolean: return -1 if largeint p is a perfect square, else zero.
'If true, set q to the square root of p. {t0..t2}
'
DECLARE FUNCTION IsPPrm% (p, d, k)
'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& (sw)
'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& ()
'Return the upper limit of primetable PrimFlgs.bin.
'
DECLARE SUB Triald (q, p, c)
'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 (q, p, a)
'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% (p, a, q, k)
'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% (p)
'Return the length of largeint p measured in words.
'
DECLARE FUNCTION Gets% (p)
'Return the sign bit of largeint p: -32768 if p < 0, else zero.
'
DECLARE FUNCTION Getw% (p, k)
'Return word number k of largeint p: a base MB digit (see below).
'
DECLARE SUB Setl (p, a)
'Set the word length of largeint p to a.
'
DECLARE SUB Sets (p, a)
'Set the sign bit of largeint p if a < 0, else clear.
'
DECLARE SUB Setw (p, k, a)
'Set word number k of largeint p to |a| < MB.
'
' ****************************************************************************
'                             Internal functions
' ****************************************************************************
DECLARE SUB Errorh (f AS STRING)
'Print an error message.
'
DECLARE SUB Lftj (p, k)
'Left-justify largeint p, beginning at word number k.
'
DECLARE FUNCTION Hp2% (p)
'Return 2^(highest set bit in the leftmost word of largeint p).
'
' ****************************************************************************
'                 BASIC prototypes for machine language support
' ****************************************************************************
DECLARE SUB SftL (c, BYVAL k)
'Left-shift long integer c by k (modulo 32) bits.
'
DECLARE SUB SftR (c, BYVAL k)
'Right-shift long integer c by k (modulo 32) bits.
'
DECLARE FUNCTION EstQ% (SEG a, c)
'Build a 45-bit dividend prefix, divide by 30-bit prefix c
'and return 15-bit trial quotient.
'
DECLARE SUB Ffix ()
'Fix up the FWAIT bug in the QB FPU emulation library.
'
' ****************************************************************************
'                             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.

Fibonacc.bas

Fibonacci numbers and the golden ratio phi.
A simple recurrent sequence, minimal module sample.

GaussEli.bas

Gauss-Bareiss fraction free elimination.
Bring matrix A into upper triangular form and solve A·X = B

LDioSyst.bas

Solution of an m x n 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.

ExtEucli.bas

Minimal solution of the Diophantine equation ax − by = c, also known as Bézout's identity.

ChineseR.bas

Solution of a system of linear congruences. Apply the Chinese remainder theorem to p simultaneous congruences a_i·x ≡ b_i Mod M_i.

QuadCngr.bas

Roots of a quadratic congruence with prime modulus N, Ressol,
the Shanks-Tonelli algorithm. For r² ≡ -1 mod N, if N ≡ 1 mod 4 also the Diophantine equation x² + y² = N is solved with Cornacchia's algorithm.

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.

FundUnit.bas

Power basis representation of the fundamental unit (p + q √d) / 2 of the real quadratic field K = Z(√d).
Also given are the continued fraction expansion and decimal approximation of the square root.

CnFcRoot.bas

Approximates a real root of integral polynomial f(x)
with Lagrange's 1769 continued fractions method.

RadpAdic.bas

p-Adic roots of an integral polynomial using Hensel's lemma.

PolDisc.bas

Discriminant of an integral poly f via resultant = determinant of Sylvester's matrix for f and f'.

Rational.bas

The continued fraction algorithm.
Simplify ratios and express decimals as common fractions.

ZnStruct.bas

Structure of the multiplicative group of units modulo N:
another Smith normal form application.

PrimRoot.bas

Test if base a belongs to phi(N) modulo N,
input a = 1 to find the least primitive root of N.

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.

PrimFlgs.bas

The Sieve of Erathostenes: a variation of Chartres' Algorithm 311.
Generates an encoded list of prime numbers read by library function Nxtprm().

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 factor 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.

SquFoFac.bas

Shanks's square form factorization (SquFoF) with queue and multipliers:
non-trivial ambiguous quadratic forms of discriminant 4kN yield factors of N.

CFracFac.bas

The 1975 Brillhart-Morrison continued-fraction factorization method (CFrac)
with multipliers, early abort strategy and large prime relations.

DiscrLog.bas

Silver-Pohlig-Hellman method for discrete logarithms, using Pollard rho.

MASH-1.bas

Cryptography: Modular Arithmetic Secure Hash function.

Encoder.bas

Cryptography: exponential encryption with salt (RSA).

Decoder.bas

Decryption using the Chinese remainder theorem.

Genrator.bas

Cryptography: probable primes and crypto key generation.
Prints relevant values to public & private key-files.

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’.

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.

PowrMtrx.bas

Pseudoprime test for N using k-th order recurrence relations,
powering companion matrix M for characteristic polynomial f(x).

LucasFun.bas

Lucas pseudoprime test, quick calculation of Fibonacci giants.

Mersenne.bas

Lucas-Lehmer test for Mersenne numbers 2^ p - 1, with odd prime p.
Toy version for small p.

WilsonTh.bas

Test for primality with Wilson's theorem:
iff p is prime (p - 1)! ≡ -1 (mod p).

Goldbach.bas

The strong Goldbach conjecture:
every even integer greater than 2 is the sum of two primes.

Reciproc.bas

Reciprocal of a small number, print decimal periods of maximal length.

Kratsuba.bas

Karatsuba multiplication, play with a recursive subroutine.

BernPoly.bas

Sums of integer powers with Bernoulli polynomials.

Pi.bas

Calculation of pi with Gauss's equation
π/4 = 12 * arctan(1/18) + 8 * arctan(1/57) - 5 * arctan(1/239).

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.

If you've written some interesting program using the library or have suggestions for improvements, let me know. If you find any bugs, please report to me here.

↑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, großen Ganzzahlen, 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.