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 only implement the most necessary operations.

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

(See Download page for Visual Basic/win, FreeBasic and XBasic versions.)

' --------------------------------------------------------------------------
'       Number array n(uj, ui) upper bounds ( $STATIC version only )
' --------------------------------------------------------------------------
'CONST ui = 5, uj = 1418 '              long lib
CONST ui = 35, uj = 402 '              medium lib
'CONST ui = 172, uj = 91 '              wide lib
' --------------------------------------------------------------------------
'                             Global constants
' --------------------------------------------------------------------------
CONST LB = 15, MB = 32768, MH = 16384 ' binary storage base
CONST M1 = 32767, M2 = 1073709056 '     masks MB - 1, (MB - 1)* MB
CONST LD = 8, MD = 100000000 '          decimal in/out base, MY < MB
CONST MY = 10000&, Fx = 2147483629 '    fixed limit
CONST t0 = -1, t1 = -2 '                pointers to swap-columns
CONST t2 = -3, t3 = -4
'
' **************************************************************************
'                      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 planning to use ('CONST p = 0, q = 1, ..'), so that
'each number 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,
'the letters p, q, r, d, m  denote 16-bit pointers to large integers;
'variables a, k, sw  are 16-bit integers; c is a 32-bit integer;
'g is a number string; f and h are any string.
'
DECLARE FUNCTION LargeInit% (k AS INTEGER, 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 <= ui in the $STATIC version) and string f = the name of the program.
'If f = "" no logfile is opened. Return value is the largeint word length.
'The maximal decimal number length is approximately 4.51 * word length.
'
DECLARE SUB Letf (p AS INTEGER, c AS LONG)
'Assign long integer c to large integer p.
'
DECLARE SUB Rndf (p AS INTEGER, k AS INTEGER)
'Assign to p a random positive number with bitlength Abs(k).
'
DECLARE SUB Term ()
'Close the primetable and logfile.
'
' **************************************************************************
'                           Conversion functions
' **************************************************************************
DECLARE SUB Readst (p AS INTEGER, g AS STRING)
'Convert decimal number string g to large integer p.
'Except for a prefixed minus sign, there's no check
'on non-digit characters in the string.
'
DECLARE FUNCTION Decf% (p AS INTEGER)
'Convert large integer p to base MD in array ln(),
'the bytelength of p is returned. Call before CnvSt.
'
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.
'
DECLARE FUNCTION Logf# (p AS INTEGER)
'Double precision natural logarithm of large integer Abs(p).
'
DECLARE SUB Ratdec (g AS STRING, p AS INTEGER, q AS INTEGER)
'Return the decimal expansion of rational p / q in [0, 1) in string g.
'
' **************************************************************************
'                             Obtaining output
' **************************************************************************
DECLARE SUB Printf (f AS STRING, g AS STRING, h AS STRING, sw AS INTEGER)
'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 AS INTEGER, f AS STRING, h AS STRING, sw AS INTEGER)
'Print large integer p with prefix f and suffix h using Printf switches.
'
DECLARE SUB Printr (p AS INTEGER, q AS INTEGER, r AS INTEGER, f AS STRING, _
 h AS STRING)
'Print rational p / q. Pass decimal length in number r with 'Setl r, len'.
'
DECLARE SUB Prints (f AS STRING, sw AS INTEGER)
'Print string f. Set switch = 1 to attach a CrLf, switch = 2 for two.
'
' **************************************************************************
'                         Basic arithmetic functions
' **************************************************************************
DECLARE FUNCTION Absf% (p AS INTEGER)
'Set the sign of number p to 1, returns the current value.
'
DECLARE SUB Add (p AS INTEGER, q AS INTEGER)
'Add p and q, the result is in p.
'
DECLARE SUB Chs (p AS INTEGER)
'Change the sign of number p.
'
DECLARE SUB Dcr (p AS INTEGER, a AS INTEGER)
'Decrease large integer p with small positive a.
'
DECLARE SUB Divd (p AS INTEGER, d AS INTEGER, q AS INTEGER)
'Divide p by d. The quotient is in q, the remainder in p.
'
DECLARE SUB Inc (p AS INTEGER, a AS INTEGER)
'Increase large integer p with small positive a.
'
DECLARE SUB Mult (p AS INTEGER, q AS INTEGER, r AS INTEGER)
'Multiply p and q. The result is in p by default,
' 'Swp p, r' puts it in r.
'
DECLARE SUB Powr (p AS INTEGER, k AS INTEGER)
'Raise p to the power k, the result is in p.
'
DECLARE SUB Squ (p AS INTEGER, q AS INTEGER)
'Return the square of p in q.
'
DECLARE SUB Subt (p AS INTEGER, q AS INTEGER)
'Subtract q from p, the result is in p.
'To retain p, use 'Subt q, p: Chs q'.
'
' **************************************************************************
'                          Comparison and copying
' **************************************************************************
DECLARE FUNCTION Cmp% (p AS INTEGER, q AS INTEGER)
'Compare returns -1 if p < q, 0 if p = q, or 1 if p > q.
'
DECLARE FUNCTION Isf% (p AS INTEGER, a AS INTEGER)
'Boolean: check if large integer p equals one-word value a.
'
DECLARE FUNCTION Ismin1% (p AS INTEGER, m AS INTEGER)
'Boolean: check if p = m - 1
'
DECLARE SUB Copyf (p AS INTEGER, q AS INTEGER)
'Copy number p to q.
'
DECLARE SUB Swp (p AS INTEGER, q AS INTEGER)
'Exchange numbers p and q.
'
' **************************************************************************
'                             Bit manipulation
' **************************************************************************
DECLARE FUNCTION Bitl% (p AS INTEGER)
'Return the bitlength of large integer p.
'
DECLARE SUB Boolf (p AS INTEGER, q AS INTEGER, k AS INTEGER)
'Bitwise Boolean functions return p:= q Op(k) p,
'with Op(1) = AND, Op(2) = XOR, Op(3) = OR.
'If k = 0 then NOT p is returned and q ignored.
'In each case, the sign of p is unaffected.
'
DECLARE SUB Lsft (p AS INTEGER, k AS INTEGER)
'Large integer p is shifted left by k bits.
'
DECLARE SUB Rsft (p AS INTEGER, k AS INTEGER)
'Large integer p is shifted right by k bits.
'
' **************************************************************************
'                        Modular arithmetic functions
' **************************************************************************
DECLARE SUB Modbal (p AS INTEGER, m AS INTEGER)
'Balanced residue p modulo m: Abs(p) <= m \ 2.
'
DECLARE SUB Moddiv (p AS INTEGER, m AS INTEGER)
'Compute positive residue p modulo m.
'
DECLARE SUB Modmlt (p AS INTEGER, q AS INTEGER, m AS INTEGER)
'Multiply p and q, then reduce modulo m. The result is in p.
'
DECLARE SUB Modpwr (p AS INTEGER, q AS INTEGER, m AS INTEGER)
'Modular exponentiation: base p, exponent q, modulus m.
'The result is in p.
'
DECLARE SUB Modsqu (p AS INTEGER, m AS INTEGER)
'Square p, then reduce modulo m.
'
DECLARE FUNCTION Mp2% (p AS INTEGER, a AS INTEGER)
'Return large integer p modulo one-word power-of-2.
'
' **************************************************************************
'                         Number theoretic functions
' **************************************************************************
DECLARE SUB Bezout (p AS INTEGER, q AS INTEGER, d AS INTEGER)
'Solution of the Diophantine equation px + qy = gcd(p, q).
'Returns positive x in p, y in q and gcd(p, q) in d.
'
DECLARE SUB Euclid (p AS INTEGER, q AS INTEGER, d AS INTEGER)
'The extended Euclidean algorithm computes
'p^ -1 Mod q in p, q/ gcd in q and gcd(p, q) in d.
'
DECLARE SUB Fctl (p AS INTEGER, a AS INTEGER)
'Computes factorial a(a-1)···2·1 in large integer p.
'
DECLARE SUB Gcd (p AS INTEGER, q AS INTEGER, d AS INTEGER)
'Returns greatest common divisor(p, q) in d.
'
DECLARE FUNCTION IsPPrm% (p AS INTEGER, d AS INTEGER, k AS INTEGER)
'Boolean: check if p is probably prime with the Miller-Rabin test.
'k is the number of witnesses, d returns a small divisor, if any.
'
DECLARE FUNCTION IsSqr% (p AS INTEGER, q AS INTEGER)
'Boolean: check if p is square. If true, the root is in q.
'
DECLARE SUB Isqrt (p AS INTEGER, q AS INTEGER)
'Integer square root of p with Heron's algorithm. Result is in q.
'
DECLARE FUNCTION Kronec% (p AS INTEGER, q AS INTEGER, d AS INTEGER)
'Kronecker's quadratic residuosity symbol (p/q) = 0, 1, or -1.
'Returns an odd gcd(p, q) in d, and 2 if it's even.
'
DECLARE SUB Lcm (p AS INTEGER, q AS INTEGER, d AS INTEGER)
'Returns least common multiple(p, q) in d.
'
DECLARE FUNCTION Nxtprm& (sw AS INTEGER)
'Loop through primetable PrimFlgs.bin. Initialize with sw = 0
'and get the next prime with sw <> 0 for each successive call.
'
DECLARE SUB Triald (q AS INTEGER, p AS INTEGER, c AS LONG)
'Trial divide p up to limit c. Primedivisors w/multiplicities
'are stored in list {q}, the remaining cofactor is returned in p.
'
' **************************************************************************
'                              Direct access
' **************************************************************************
DECLARE SUB EnQ (q AS INTEGER, p AS INTEGER, a AS INTEGER)
'Store small p and one-word Abs(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 AS INTEGER, a AS INTEGER, q AS INTEGER, _
 k AS INTEGER)
'Get successive numbers p and Abs(a) from list {q} at offset k.
'Set k = 0 to initialize reading, returns zero if end-of-list.
'
DECLARE FUNCTION Gete% (p AS INTEGER, k AS INTEGER)
'Return array element k of number p. This is a base MB digit.
'
DECLARE FUNCTION Getl% (p AS INTEGER)
'Return the length of number p measured in array elements.
'
DECLARE FUNCTION Gets% (p AS INTEGER)
'Return the sign of number p, Gets = -1 if p < 0, 1 otherwise.
'
DECLARE SUB Sete (p AS INTEGER, k AS INTEGER, a AS INTEGER)
'Set binary array element k of number p to Abs(a).
'
DECLARE SUB Setl (p AS INTEGER, a AS INTEGER)
'Set the length of number array p to Abs(a).
'
DECLARE SUB Sets (p AS INTEGER, a AS INTEGER)
'Set the sign of number p to Sgn(a + 0.5).
'
' **************************************************************************
'                            Internal functions
' **************************************************************************
DECLARE SUB Errorh (f AS STRING)
'Print an error message.
'
DECLARE FUNCTION Hp2% (p AS INTEGER)
'Return 2^ (highest set bit in the leftmost word of p).
'
DECLARE SUB Lftj (p AS INTEGER, k AS INTEGER)
'Left-justify large integer p on from array element k.
'
DECLARE FUNCTION PrimCeil& ()
'Return the upper limit of primetable PrimFlgs.bin.
'
DECLARE SUB Work (f AS STRING)
'Return working directory in ENVIRON$("LARGEINT"),
'buffer f should be large enough to hold the path.
'
' **************************************************************************
'               BASIC prototypes for machine language support
' **************************************************************************
DECLARE FUNCTION EstQ% (SEG a AS INTEGER, c AS LONG)
'( n(j, ..)·MB^2 + n(j-1, ..)·MB + n(j-2, ..) ) \ c
'45-bit dividend \ 30-bit divisor.
'
DECLARE SUB Ffix ()
'Fix up the FWAIT bug in the Microsoft FPU emulation library.
'
DECLARE SUB SftL (c AS LONG, k AS INTEGER)
'ShL long integer c over k (mod 32) bits.
'
DECLARE SUB SftR (c AS LONG, k AS INTEGER)
'ShR long integer c over k (mod 32) bits.
'
' **************************************************************************

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 (m ≤ n) linear Diophantine system
using Hermite and Smith normal forms.

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
for the quadratic ring K = Z(√d). Also given are the continued fraction expansion
of the irrational and a decimal approximation of the real 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.

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).
If the cycle of reduced forms contains a non-trivial ambiguous ideal,
the discriminant can be factorized.

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 ‘N-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.

Pi.bas

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

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

Carl Friedrich Gauss (1777-1855)

Download

View the ReadMe.txt

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.

Learn about Number Theory at Eric Weisstein's great
World of Mathematics or the online Wikipedia encyclopedia.

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