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.
' ************************************************************************* '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 ' ' ****************************************************************************************
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.
Mirror site (auf deutsch)
Computing Ackermann's function in Basic.
This site has an annex: Finite field arithmetic in BASIC.
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.