' ****************************************************************************
'Subject: Basic library include file for polynomial arithmetic modulo p.
' (version exclusively for characteristic 2 in packet PolMods.zip)
'Author : Sjoerd.J.Schaper
'Date : 01-15-2008
'Code : FreeBasic 0.18.3 (QuickBasic 7.1 adaptation in packet PolMods.zip)
'
#INCLIB "polymodp"
' ****************************************************************************
CONST LB = 15, nmx = 256, pmx = 27 ' max.degree, primes
CONST MB = 32768, M1 = 32767 ' bigint base, mask
TYPE Poly FIELD = 1 ' univariate polynomial:
n AS SHORT ' degree and
c(nmx) AS INTEGER ' coefficients
END TYPE
TYPE Bigi FIELD = 1 ' big integer:
l AS SHORT ' length,
s AS SHORT ' sign and
d(nmx) AS SHORT ' base 2^ 15 digits
END TYPE
DIM SHARED Modulus AS SHORT ' field characteristic < 32768
'
' ****************************************************************************
' Initialization and assignment
' ****************************************************************************
DECLARE FUNCTION Setchar alias "Setchar" (ByVal p AS SHORT) AS SHORT
'Set prime field characteristic, returns 0 if composite(p)
DECLARE SUB Letc alias "Letc" (f AS Poly, ByVal c AS SHORT)
'Set f to constant c
DECLARE SUB Letm alias "Letm" (f AS Poly, ByVal c AS SHORT)
'Set f to monomial X + c
DECLARE SUB Rndp alias "Rndp" (f AS Poly, ByVal n AS SHORT)
'f:= random monic polynomial of degree n
DECLARE SUB Leti alias "Leti" (a AS Bigi, ByVal b AS INTEGER)
'Let bigint a = b
'
' ****************************************************************************
' In- and output
' ****************************************************************************
DECLARE FUNCTION Parsp alias "Parsp" (f AS Poly, ByVal g AS STRING) AS SHORT
'Parse input polynomial
DECLARE SUB Prntp alias "Prntp" (f AS Poly, ByVal g AS STRING)
'Pretty print polynomial with prefix g$
DECLARE SUB Readp alias "Readp" (f AS Poly, ByVal n AS SHORT)
'Input coefficient vector f.c[n..0]
DECLARE SUB CRTfile alias "CRTfile" (ByVal sw AS SHORT)
'sw = -1 opens, sw = 0 closes CRTtrans-inputfile
DECLARE SUB CRTchar alias "CRTchar" (ByVal p AS SHORT)
'Fast field characteristic, without primality test
DECLARE SUB CRTprnt alias "CRTprnt" (f() AS Poly, ByVal g AS STRING)
'Print polynomial array f() to open file #Tfile,
'creates an inputfile for LargeInt-module CRTtrans.bas
'
' ****************************************************************************
' Polynomial arithmetic functions
' ****************************************************************************
DECLARE SUB Addp alias "Addp" (f AS Poly, g AS Poly)
'Polynomial f:= f + g
DECLARE SUB Divdp alias "Divdp" (f AS Poly, g AS Poly, q AS Poly)
'Division f = g * q + r; f:= r, q:= f / g
DECLARE SUB Multp alias "Multp" (f AS Poly, g AS Poly)
'f:= f * g
DECLARE SUB Squp alias "Squp" (f AS Poly)
'f:= f ^ 2
DECLARE SUB Subtp alias "Subtp" (f AS Poly, g AS Poly)
'f:= f - g
DECLARE SUB Shlp alias "Shlp" (f AS Poly, ByVal k AS SHORT)
'f:= f * X^ k
DECLARE SUB Shrp alias "Shrp" (f AS Poly, ByVal k AS SHORT)
'f:= f / X^ k
'
' ****************************************************************************
' Modulo-polynomial functions
' ****************************************************************************
DECLARE SUB Invp alias "Invp" (f AS Poly, g AS Poly)
'f:= f^ -1 Mod g
DECLARE SUB Modmltp alias "Modmltp" (f AS Poly, g AS Poly, h AS Poly)
'Polynomial f:= f * g Mod h
DECLARE SUB Modpwrp alias "Modpwrp" (f AS Poly, k AS Bigi, g AS Poly)
'f:= f^ bigint(k) Mod g, set g = 0 to skip reduction
'
' ****************************************************************************
' Checks
' ****************************************************************************
DECLARE FUNCTION Ispc alias "Ispc" (f AS Poly, ByVal c AS SHORT) AS SHORT
'Is f = constant c ?
DECLARE FUNCTION Ispm alias "Ispm" (f AS Poly, ByVal c AS SHORT) AS SHORT
'Is f = monomial X + c ?
DECLARE FUNCTION Isirrd alias "Isirrd" (f AS Poly) AS SHORT
'Is f irreducible ?
'
' ****************************************************************************
' Assorted functions
' ****************************************************************************
DECLARE SUB Bezoup alias "Bezoup" (f AS Poly, g AS Poly, d AS Poly)
'Solution of the Diophantine equation fa + gb = gcd(f, g).
'Returns a in f, b in g, and gcd(f, g) in d
DECLARE SUB Deriv alias "Deriv" (f AS Poly, g AS Poly)
'g:= derivative f'
DECLARE FUNCTION Disc alias "Disc" (f AS Poly) AS SHORT
'Return discriminant(f)
DECLARE FUNCTION Eval alias "Eval" (f AS Poly, ByVal c AS INTEGER) AS SHORT
'Evaluate f(c) with Horner's rule
DECLARE SUB Gcdp alias "Gcdp" (f AS Poly, g AS Poly, d AS Poly)
'd:= polynomial gcd(f, g)
DECLARE SUB Monic alias "Monic" (f AS Poly)
'Make monic(f)
DECLARE FUNCTION Psdiv alias "Psdiv" (f AS Poly, g AS Poly, q AS Poly)_
AS SHORT ' Pseudo-division,
'returns the scalar L(g)^(f.n - g.n + 1) * f = g * q + r
DECLARE SUB Recip alias "Recip" (f AS Poly)
'f:= reciprocal X^f.n * f(1/ X)
DECLARE FUNCTION Result alias "Result" (f AS Poly, g AS Poly) AS SHORT
'Return resultant(f, g)
DECLARE SUB Subst alias "Subst" (f AS Poly, g AS Poly)
'Horner substitution f:= f(g(X))
'
' ****************************************************************************
' Modular arithmetic functions
' ****************************************************************************
DECLARE FUNCTION Modd alias "Modd" (ByVal c AS INTEGER) AS SHORT
'Return positive(c) mod global modulus
DECLARE FUNCTION Balc alias "Balc" (ByVal c AS INTEGER) AS SHORT
'Return least absolute value(c) mod global modulus
DECLARE FUNCTION Inv alias "Inv" (ByVal c AS INTEGER) AS SHORT
'Return 1/ c mod global modulus
DECLARE FUNCTION Modpwr alias "Modpwr" (ByVal a AS INTEGER, ByVal k AS SHORT)_
AS SHORT ' Return a^ k mod global modulus
'
' ****************************************************************************
' Building big integer exponents
' ****************************************************************************
DECLARE SUB Addi alias "Addi" (a AS Bigi, ByVal c AS SHORT)
'bigint a:= a + c
DECLARE FUNCTION Divi alias "Divi" (a AS Bigi, ByVal c AS SHORT) AS SHORT
'a:= a / c, returns remainder
DECLARE SUB Muli alias "Muli" (a AS Bigi, ByVal c AS SHORT)
'a:= a * c
DECLARE SUB Pwri alias "Pwri" (p AS Bigi, ByVal a AS SHORT, ByVal k AS SHORT)
'p:= a ^ k
DECLARE SUB Squi alias "Squi" (a AS Bigi)
'a:= a ^ 2
'
' ****************************************************************************
' Internal functions
' ****************************************************************************
DECLARE SUB Adjp alias "Adjp" (f AS Poly, ByVal n AS SHORT)
'Adjust degree(n) of polynomial f
DECLARE SUB Lftj alias "Lftj" (a AS Bigi, ByVal k AS SHORT)
'Left-justify bigint(a) on from array element(k)
DECLARE SUB Mulc alias "Mulc" (f AS Poly, ByVal c AS SHORT)
'f:= f * c mod global modulus
DECLARE SUB Triald alias "Triald" (Pf() AS INTEGER, ByRef fx AS SHORT,_
ByVal n AS INTEGER) ' Return primefactors(n) in Pf(fx, k)
DECLARE SUB Work alias "Work" (ByRef f AS STRING)
'Return working directory in ENVIRON$("LARGEINT"),
'buffer f should be large enough to hold the path
'
' ****************************************************************************
'28 largest primes < 2^ 15 for CRT-transform:
' ****************************************************************************
DATA 32749,32719,32717,32713,32707,32693,32687
DATA 32653,32647,32633,32621,32611,32609,32603
DATA 32587,32579,32573,32569,32563,32561,32537
DATA 32533,32531,32507,32503,32497,32491,32479
All this functionality in only 16 Kb
References: