' **************************************************************************
'Subject: Basic library include file for polynomial arithmetic modulo p.
' (version exclusively for characteristic 2 in packet PolMods.zip)
'Author : Sjoerd.J.Schaper
'Update : 11-01-2007
'Code : PDS (QuickBasic 7.1), FreeBasic adaptation in packet PolMods.zip
' **************************************************************************
CONST nmx = 144, pmx = 27 ' max.degree, primes
CONST LB = 15, MB = &H8000& ' bigint base
DEFINT A-W: DEFLNG X-Z
TYPE Poly ' univariate polynomial:
n AS INTEGER ' degree and
c(nmx) AS LONG ' coefficients
END TYPE
TYPE Bigi ' big integer:
l AS INTEGER ' length,
s AS INTEGER ' sign and
d(nmx) AS INTEGER ' base 2^ 15 digits
END TYPE
DIM SHARED Modulus AS INTEGER ' field characteristic < 32768
'
' **************************************************************************
' Initialization and assignment
' **************************************************************************
DECLARE FUNCTION Setchar% (p)
'Set prime field characteristic, returns 0 if composite(p)
DECLARE SUB Letc (f AS Poly, c)
'Set f to constant c
DECLARE SUB Letm (f AS Poly, c)
'Set f to monomial X + c
DECLARE SUB Rndp (f AS Poly, n, m)
'f:= random monic polynomial of degree n
'with m non-zero coefficients
DECLARE SUB Leti (a AS Bigi, b AS LONG)
'Let bigint a = b
'
' **************************************************************************
' In- and output
' **************************************************************************
DECLARE FUNCTION Parsp% (f AS Poly, g AS STRING)
'Parse input polynomial
DECLARE SUB Prntp (f AS Poly, g AS STRING)
'Pretty print polynomial with prefix g$
DECLARE SUB Readp (f AS Poly, n)
'Input coefficient vector f.c[n..0]
DECLARE FUNCTION CRTfile% (p)
'Open an inputfile for LargeInt-module CRTtrans.bas
'set p = number of CRT-primes (<= pmx), p = 0 to close
DECLARE SUB CRTprnt (f() AS Poly, sw, g AS STRING)
'Print polynomial array f() to open file #Tfile,
'set sw to divide out content(f), g is a comment string
DECLARE SUB CRTchar (p)
'Fast field characteristic, without primality test
'
' **************************************************************************
' Polynomial arithmetic functions
' **************************************************************************
DECLARE SUB Addp (f AS Poly, g AS Poly)
'Polynomial f:= f + g
DECLARE SUB Divdp (f AS Poly, g AS Poly, q AS Poly)
'Division f = g * q + r; f:= r, q:= f / g
DECLARE SUB Multp (f AS Poly, g AS Poly)
'f:= f * g
DECLARE SUB Squp (f AS Poly)
'f:= f ^ 2
DECLARE SUB Subtp (f AS Poly, g AS Poly)
'f:= f - g
DECLARE SUB Shlp (f AS Poly, k)
'f:= f * X^ k
DECLARE SUB Shrp (f AS Poly, k)
'f:= f / X^ k
'
' **************************************************************************
' Modulo-polynomial functions
' **************************************************************************
DECLARE SUB Invp (f AS Poly, g AS Poly, d AS Poly)
'f:= f^ -1 Mod g, g:= g / gcd(f, g), d:= gcd(f, g)
DECLARE SUB Modmltp (f AS Poly, g AS Poly, h AS Poly)
'Polynomial f:= f * g Mod h
DECLARE SUB 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% (f AS Poly, c)
'Is f = constant c ?
DECLARE FUNCTION Ispm% (f AS Poly, c)
'Is f = monomial X + c ?
DECLARE FUNCTION Isirrd% (f AS Poly)
'Is f irreducible ?
'
' **************************************************************************
' Assorted functions
' **************************************************************************
DECLARE SUB 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 (f AS Poly, g AS Poly)
'g:= derivative f'
DECLARE FUNCTION Disc% (f AS Poly)
'Return discriminant(f)
DECLARE FUNCTION Eval% (f AS Poly, c AS LONG)
'Evaluate f(c) with Horner's rule
DECLARE SUB Gcdp (f AS Poly, g AS Poly, d AS Poly)
'd:= polynomial gcd(f, g)
DECLARE SUB Monic (f AS Poly)
'Make monic(f)
DECLARE FUNCTION Psdiv% (f AS Poly, g AS Poly, q AS Poly)
'Pseudo-division, returns the scalar
'L(g)^(f.n - g.n + 1) * f = g * q + r
DECLARE SUB Recip (f AS Poly)
'f:= reciprocal X^f.n * f(1/ X)
DECLARE FUNCTION Result% (f AS Poly, g AS Poly)
'Return resultant(f, g)
DECLARE SUB Subst (f AS Poly, g AS Poly)
'Horner substitution f:= f(g(X))
'
' **************************************************************************
' Modular arithmetic functions
' **************************************************************************
DECLARE FUNCTION Modd% (c AS LONG)
'Return positive(c) mod global modulus
DECLARE FUNCTION Balc% (c AS LONG)
'Return least absolute value(c) mod global modulus
DECLARE FUNCTION Inv% (c AS LONG)
'Return 1/ c mod global modulus
DECLARE FUNCTION Modpwr% (a AS LONG, k)
'Return a^ k mod global modulus
'
' **************************************************************************
' Building big integer exponents
' **************************************************************************
DECLARE SUB Addi (a AS Bigi, c)
'bigint a:= a + c
DECLARE FUNCTION Divi% (a AS Bigi, c)
'a:= a / c, returns remainder
DECLARE SUB Muli (a AS Bigi, c)
'a:= a * c
DECLARE SUB Pwri (p AS Bigi, a, k)
'p:= a ^ k
DECLARE SUB Squi (a AS Bigi)
'a:= a ^ 2
'
' **************************************************************************
' Internal functions
' **************************************************************************
DECLARE SUB Adjp (f AS Poly, n)
'Adjust degree n of polynomial f
DECLARE SUB Lftj (a AS Bigi, k)
'Left-justify bigint a, starting at array element k
DECLARE SUB Mulc (f AS Poly, c)
'f:= f * c mod global modulus
DECLARE SUB Negp (f AS Poly)
'f:= -f mod global modulus
DECLARE SUB PosLc (f AS Poly)
'Make positive leading coefficient
DECLARE FUNCTION Tail% (f AS Poly, g AS Poly)
'Zeroed poly tail for Addp and Subtp
DECLARE SUB Triald (Pf() AS LONG, fx, n AS LONG)
'Return primefactors of n in Pf(fx, k)
DECLARE SUB Work (f AS STRING)
'Set f to the working directory in ENVIRON$("LARGEINT"),
'the stringbuffer must 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 15 Kb
The above seven modules need transformation module CRTtrans.bas
to lift the results from the intersection of 28 small equivalence classes:
Plus GF(2^8) versions of the above three Elgamal cryptography modules.
References: