Algorithmic Number Theory for Cryptography and Cryptanalysis: Primality, Factoring and Discrete Logarithms - Introduction to Cryptography with Maple - page 388

i:=0;
while i<k+addrels do
x := RandomTools:-MersenneTwister:-GenerateInteger(range = 1 .. p-1);
rels := [op(rels), ICTrialDivision(p, g, x, fb)];
i := nops(rels)
end do;
M := LinearAlgebra:-Modular:-Mod(p-1, rels, integer);
if verbose then
s := interface(rtablesize = k + addrels + 1);
printf("Solving the linear system modulo %d of augmented matrix:\n", p-1);
print(M)
end if;
LinearAlgebra:-Modular:-LinearSolve(p-1, M, 1);
fblogs := LinearAlgebra:-Modular:-Copy(p-1, M, 1 .. k, k+1);
if verbose then
printf("Discrete logs of the factor base primes:\n");
print(LinearAlgebra:-Modular:-Copy(p-1, fblogs, 'transpose'));
interface(rtablesize = s)
end if;
[fb, fblogs]
end proc:

Next we give the main function of the algorithm, namely the function
IndexCalculus below. This function has three required input parameters that
serve to specify the prime p , the primitive root g , and the element a

∈ Z p whose
discrete logarithm is to be computed. The optional input parameter, add is used to
specify the number of relations to be computed, which will be equal to the size of
the factor base plus the value of add (this value will be subsequently passed to the
function IndexCalculusPrecomp ). The remaining optional input parameters
are the keyword parameters verbose and pre . The purpose of the first is to tell
the function whether to print some additional information or not (with default value
false ). Finally, pre is used to feed the function with the result of the precomputa-
tion. The default value calls IndexCalculusPrecomp but, if the result of calling
this function is already known, it can be passed directly to pre without having to
redo the precomputation. The output is log g a

∈ Z p − 1 .

> IndexCalculus := proc(p::posint, g::posint, a::posint, add:=10,
{verbose:=false, pre:=IndexCalculusPrecomp(p, g, addrels=add)})
uses LinearAlgebra:-Modular;
local fb, fblogs, k, rels, s, r, m;
RandomTools:-MersenneTwister:-SetState();
fb := pre[1];
fblogs := pre[2];
k := nops(fb);
if fblogs[k] = 0 then
error "Insufficient number of relations"
end if;
rels := [];
while nops(rels) = 0 do
r := RandomTools:-MersenneTwister:-GenerateInteger(range = 2 .. p-2);
rels := [op(rels), ICTrialDivision(p, g, r, fb, true, a)]
end do;
if verbose then
printf("(%dˆ%d)*%d mod %d has exponent vector over the factor base:\n",g,r,a,p);
s := interface(rtablesize =k+1);
print(op(rels));
interface(rtablesize = s)
end if;
m := Multiply(p-1, Mod(p-1, op(rels), integer), fblogs);