Motivation
Motivation
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
'Update : 06-06-2009
'Code : FreeBasic 0.20.0 adaptation
'
' ***************************************************************************************
REM #DEFINE VB
' : Remove this REMark or put "-d VB" on the fbc command line
' : to compile a Visual Basic-callable DLL.
'
#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 alias "LargeInit" (k AS SHORT, 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 maximal largeint length measured in array elements
'(='words'). The decimal number length is about 4.51 * largeint length.
'
DECLARE SUB Letf alias "Letf" (p AS SHORT, c AS INTEGER)
'Assign signed long integer c to large integer p.
'
DECLARE SUB Readst alias "Readst" (p AS SHORT, g AS ZSTRING PTR)
'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 alias "Rndf" (p AS SHORT, k AS SHORT)
'Assign a random positive value of bitlength k to large integer p.
'
DECLARE SUB Term alias "Term" ()
'Close the primetable and logfile.
'
' ***************************************************************************************
' Conversion functions
' ***************************************************************************************
DECLARE FUNCTION Decf alias "Decf" (p AS SHORT) AS INTEGER
'Convert largeint p to base MD (see below), return the bytelength of
'decimal(p), call before CnvSt. {t1..t2}
'
DECLARE SUB CnvSt alias "CnvSt" (g AS ZSTRING PTR)
'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 alias "Logf" (p AS SHORT) AS DOUBLE
'Return double precision natural logarithm of |largeint p|.
'
DECLARE SUB Ratdec alias "Ratdec" (g AS ZSTRING PTR, p AS SHORT, q AS SHORT)
'Convert the decimal expansion of rational p/q in [0, 1) to string. {t0..t2}
'
' ***************************************************************************************
' Obtaining output
' ***************************************************************************************
#IFNDEF VB
'To be replaced with Visual Basic-tailored printfunctions:
' --------------------------------------------------------------------------
DECLARE SUB Printf alias "Printf" (ByRef f AS STRING, ByRef g AS STRING, _
ByRef h AS STRING, sw AS SHORT)
'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 alias "Printn" (p AS SHORT, ByRef f AS STRING, ByRef h AS STRING, _
sw AS SHORT)
'Print largeint p with prefix f and suffix h using the Printf switches. {t1..t2}
'
DECLARE SUB Printr alias "Printr" (p AS SHORT, q AS SHORT, ByRef f AS STRING, k AS SHORT)
'Print the decimal expansion of p/q with prefix f and fraction length k. {t1..t2}
'
DECLARE SUB Prints alias "Prints" (ByRef f AS STRING, sw AS SHORT)
'Print string f. Set switch = 1 to attach a CrLf, switch = 2 for two.
' --------------------------------------------------------------------------
#ENDIF
DECLARE SUB Work alias "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 alias "Subt" (p AS SHORT, q AS SHORT)
'Set large integer p to p - q.
'
DECLARE SUB Add alias "Add" (p AS SHORT, q AS SHORT)
'Set large integer p to p + q.
'
DECLARE SUB Inc alias "Inc" (p AS SHORT, a AS SHORT)
'Increment largeint p by signed one-word value a.
'
DECLARE SUB Divd alias "Divd" (p AS SHORT, d AS SHORT, q AS SHORT)
'Divide p by d, set p to the remainder and q to the quotient. {t2}
'
DECLARE FUNCTION Divint alias "Divint" (p AS SHORT, a AS SHORT) AS SHORT
'Divide largeint p by |one-word a|,
'return |remainder| and set p to the quotient.
'
DECLARE SUB Ldiv alias "Ldiv" (p AS SHORT, d AS SHORT)
'Set p to quotient p / d. {t1..t2}
'
DECLARE SUB Mult alias "Mult" (p AS SHORT, q AS SHORT, r AS SHORT)
'Set large integer r to p * q.
'
DECLARE SUB Lmul alias "Lmul" (p AS SHORT, q AS SHORT)
'Set p to product p * q. {t2}
'
DECLARE SUB Squ alias "Squ" (p AS SHORT, q AS SHORT)
'Set q to the square of p, faster than p * p.
'
DECLARE SUB Isqrt alias "Isqrt" (p AS SHORT, q AS SHORT)
'Set q to the truncated integer part of square root(p). {t0..t2}
'
DECLARE SUB Chs alias "Chs" (p AS SHORT)
'Change the sign of largeint p.
'
DECLARE FUNCTION Absf alias "Absf" (p AS SHORT) AS SHORT
'Clear sign bit of largeint p and return the current value.
'
' ***************************************************************************************
' Copying and comparison
' ***************************************************************************************
DECLARE SUB Dup alias "Dup" (p AS SHORT, q AS SHORT)
'Copy largeint p to q.
'
DECLARE SUB Swp alias "Swp" (p AS SHORT, q AS SHORT)
'Exchange largeint's p and q.
'
DECLARE FUNCTION Cmp alias "Cmp" (p AS SHORT, q AS SHORT) AS SHORT
'Compare returns -1 if p < q, 0 if p = q, or 1 if p > q.
'
DECLARE FUNCTION Isf alias "Isf" (p AS SHORT, a AS SHORT) AS SHORT
'Boolean: check if largeint p equals one-word value a.
'
DECLARE FUNCTION Ismin1 alias "Ismin1" (p AS SHORT, m AS SHORT) AS SHORT
'Boolean: check if p = m - 1
'
' ***************************************************************************************
' Bit manipulation
' ***************************************************************************************
DECLARE SUB Boolf alias "Boolf" (p AS SHORT, q AS SHORT, 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 unaffected.
'
DECLARE SUB Lsft alias "Lsft" (p AS SHORT, k AS SHORT)
'Left-shift largeint p by k bits; k < 0 shifts |k| full words.
'
DECLARE SUB Rsft alias "Rsft" (p AS SHORT, k AS SHORT)
'Right-shift largeint p by k bits; k < 0 shifts |k| full words.
'
DECLARE FUNCTION Odd alias "Odd" (p AS SHORT) AS SHORT
'Remove trailing zero bits from largeint p, return (right) shift count.
'
DECLARE FUNCTION Bitl alias "Bitl" (p AS SHORT) AS INTEGER
'Return the bitlength of largeint p.
'
' ***************************************************************************************
' Modular arithmetic functions
' ***************************************************************************************
DECLARE SUB Moddiv alias "Moddiv" (p AS SHORT, m AS SHORT)
'Compute positive residue p modulo largeint m. {t1..t2}
'
DECLARE SUB Modbal alias "Modbal" (p AS SHORT, m AS SHORT)
'Balanced residue p modulo m: |p| <= m / 2. {t1..t2}
'
DECLARE FUNCTION Mp2 alias "Mp2" (p AS SHORT, a AS SHORT) AS SHORT
'Return largeint p modulo a, which must be a one-word power of 2.
'
DECLARE SUB Modmlt alias "Modmlt" (p AS SHORT, q AS SHORT, m AS SHORT)
'Set largeint p to p * q (modulo m). {t1..t2}
'
DECLARE SUB Modsqu alias "Modsqu" (p AS SHORT, m AS SHORT)
'Set largeint p to p ^ 2 (modulo m). {t1..t2}
'
DECLARE SUB Modpwr alias "Modpwr" (p AS SHORT, q AS SHORT, m AS SHORT)
'Set largeint p to base p ^ exponent q (modulo m). {t0..t2, t3 if q < 0}
'
' ***************************************************************************************
' Number theoretic functions
' ***************************************************************************************
DECLARE SUB Powr alias "Powr" (p AS SHORT, k AS SHORT)
'Raise largeint p to the power |one-word k|. {t0..t3}
'
DECLARE SUB Fctl alias "Fctl" (p AS SHORT, a AS SHORT)
'Set largeint p to factorial a(a-1)贩. {t1..t2}
'
DECLARE SUB Lcm alias "Lcm" (p AS SHORT, q AS SHORT, d AS SHORT)
'Set largeint d to the least common multiple (p, q). {t0..t2}
'
DECLARE SUB Gcd alias "Gcd" (p AS SHORT, q AS SHORT, d AS SHORT)
'Set largeint d to the greatest common divisor (p, q). {t0..t2}
'
DECLARE SUB Euclid alias "Euclid" (p AS SHORT, q AS SHORT, 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 alias "Bezout" (p AS SHORT, q AS SHORT, 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 FUNCTION Kronec alias "Kronec" (p AS SHORT, q AS SHORT) AS SHORT
'Return Kronecker's quadratic residuosity symbol (p/q) = 0, 1, or -1. {t0..t3}
'
DECLARE FUNCTION IsSqr alias "IsSqr" (p AS SHORT, q AS SHORT) AS SHORT
'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 alias "IsPPrm" (p AS SHORT, d AS SHORT, 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 alias "Nxtprm" (ByRef sw AS SHORT) AS INTEGER
'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 alias "PrimCeil" () AS INTEGER
'Return the upper limit of primetable PrimFlgs.bin.
'
DECLARE SUB Triald alias "Triald" (q AS SHORT, p AS SHORT, c AS INTEGER)
'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 alias "EnQ" (q AS SHORT, p AS SHORT, 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 alias "ExQ" (p AS SHORT, ByRef a AS SHORT, 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 alias "Getl" (p AS SHORT) AS SHORT
'Return the length of largeint p measured in words.
'
DECLARE FUNCTION Gets alias "Gets" (p AS SHORT) AS SHORT
'Return the sign bit of largeint p: -32768 if p < 0, else zero.
'
DECLARE FUNCTION Getw alias "Getw" (p AS SHORT, k AS SHORT) AS SHORT
'Return word number k of largeint p: a base MB digit (see below).
'
DECLARE SUB Setl alias "Setl" (p AS SHORT, a AS SHORT)
'Set the word length of largeint p to a.
'
DECLARE SUB Sets alias "Sets" (p AS SHORT, a AS SHORT)
'Set the sign bit of largeint p if a < 0, else clear.
'
DECLARE SUB Setw alias "Setw" (p AS SHORT, k AS SHORT, a AS SHORT)
'Set word number k of largeint p to |a| < MB.
'
' ***************************************************************************************
' Internal functions
' ***************************************************************************************
DECLARE SUB Errorh alias "Errorh" (ByRef f AS STRING)
'Print an error message.
#IFDEF __FB_WIN32__
Declare Function ErrMsg lib "User32" alias "MessageBoxA" (hwnd As Integer, _
lpText As Zstring Ptr, lpCaption As Zstring Ptr, wType As Integer) As Integer
#ENDIF
'
DECLARE SUB Lftj alias "Lftj" (p AS SHORT, k AS SHORT)
'Left-justify largeint p, beginning at word number k.
'
DECLARE FUNCTION Hp2 alias "Hp2" (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
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.
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.
References:
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.
