Because of its abstractness and generality, my previous modular Mul command wasn't efficient for its most-common use (or most-common use relevant to this thread): the product of a range of integers (such as a factorial). So, I've now implemented Mul(m::integer..n::integer) mod p::posint as syntax additional to the Mul(e::algebraic, k::name= m::integer..n::integer) mod p::posint that I'd already implemented. This brings my code into closer alignment with VV's Fac procedure.
You don't need to use or understand any Mul syntax to use this code for its primary purpose: computing binomial(n,k) modulo a prime. It's all "under the hood". You simply use the syntax
Binomial(n,k) mod p.
The code is now (code below) in the form of a module. Simply entering the code into a Maple session or reading it in from a file will install the `mod/...` procedures. If the module is stored in a library, then simply invoking Binomial_mod_p() will install them.
Here is the code. For convenience, it also exists in the Code Edit Region of the attached worksheet.
Binomial_mod_p:= module()
description "Some tools to speed up the computation of binomial(n,k) mod p, p::prime";
option
date= "9-Nov-2018",
authors= (
"VV from MaplePrimes <...>",
"Carl Love <carl.j.love@gmail.com> (Corresponding)"
),
reference= "https://en.wikipedia.org/wiki/Lucas%27s_theorem"
;
################################################################################################################
# This is a module of tools whose primary purpose is to speed up the computation of binomial(n,k) mod p, for #
# prime p. This code interacts with `mod` through the `mod/...` global hooks (an ancient pre-module #
# interface). So, one simply invokes Binomial(n,k) mod p, where the Binomial is an inert function (see ?mod). #
# #
# The speed of the computation, and the choice of algorithm, are based on the relative sizes of n, k, and p. #
# all cases where p <= binomial(n,k)*k!, this code makes some significant attempt to use modular arithmetic to #
# improve upon the speed of the naive command binomial(n,k) mod p. #
# #
# This module has no exports. It may "installed" any of these ways: #
# 1. Simply copy-and-paste this code into an execution group and execute; #
# 2. Use the `read` command to import it from a file; #
# 3. If it's stored in a library, execute the module's name: Binomial_mod_p(). #
# Any of these techniques will assign the necessary global `mod/...` procedures. #
# #
# The module has a few userinfo statements, all keyed to infolevel[Binomial]:= 1. #
################################################################################################################
local
#The methods implemented as `mod/...` procedures by this module are Mul and Binomial. If methods are added,
#then their names should local procedures in this module named X such that the corresponding `mod/...`
#procedure will be named `mod/X` (and thus the usage will be X(...) mod p). The names should be added to
#the list immediately below. Everything else will then happen automatically.
ModMethods:= [Mul, Binomial],
#The ModuleApply only exists so that invoking Binomial_mod_p() will force the ModuleLoad to
#run if the module is stored in a library.
ModuleApply:= proc($)
userinfo(1, Binomial, "The `mod/...` procedures have been protected and installed.");
return
end proc,
ModuleUnload:= proc()
userinfo(1, Binomial, "The module Binomial_mod_p is being garbage collected.");
return
end proc,
#Some constants specifying the largest "convenient" integer for various modes of computation:
Mhf, Mgmp, Mcomp,
#...and a proc to initialize them:
Init_Ms:= proc()
#In hfloat mode, we try to keep the mantissa half full so that products don't overflow it:
Mhf:= isqrt(Scale2(1, trunc(evalhf(DBL_MANT_DIG))-1));
#In GMP mode, we try to keep products to 1-word integers.
Mgmp:= kernelopts(maximmediate);
#Likewise for compiled mode, but the wordsize is slightly different.
Mcomp:= Scale2(1, kernelopts(wordsize)-1) - 1;
return
end proc,
###########################################################
# The forms of Mul currently implemented are
# - Mul(m..n):
# a range of integers;
# - Mul(e, k= m..n):
# an expression e (likely depending on name k), with k iterated over a range of integers.
Mul:= proc()
option remember;
description "like mul, but in modular arithmetic";
if nargs=2 then
if args[1]::range then `Mul/range`(args)
else `Mul/container`(args)
fi
elif nargs=3 then `Mul/2arg`(args)
else error "Expecting 1 or 2 arguments, but received", nargs-1
fi
end proc,
`Mul/range/compiled`:=
proc(m::integer, n::integer, p::posint)::integer;
option autocompile;
local k::integer, r::nonnegint:= 1;
for k from m to n do
r:= r*k mod p;
if r = 0 then break fi
od;
r
end proc,
`Mul/range`:= proc(R::range(integer), p::posint)
description
"Computes the modular product of a range of integers (just like 1-argument mul, but modular)."
" The modulus needn't be prime. Attempts compiled mode for word-sized p."
;
local r, k, m:= lhs(R), n:= rhs(R), M, MR:= max(1, CopySign(m,1), CopySign(n,1));
if p < Mcomp then #Use compiled code.
try return `Mul/range/compiled`(m, n, p) mod p
catch: userinfo(1, Binomial, "Compiled code failed. Will use endogenous arithmetic instead.")
end try
fi;
r:= 1; M:= max(p, iquo(Mgmp, max(1, MR)));
for k from `if`(m>=0, m, -n) to `if`(m>=0, n, -m) while r <> 0 do
r:= r*k;
if r >= M then r:= irem(r, p) fi
od;
`if`(m>0 or (k-m)::even, r, -r) mod p
end proc,
`Mul/2arg/evalhf`:= proc(e, m, n, p)
local k, r:= 1, s:= 1, M:= _M, M2:= M^2;
for k from m to n while r <> 0 do
r:= r*e(k);
if r > M2 then error fi; #too large for exact representation
if r < 0 then r:= -r; s:= -s fi;
if r >= M then r:= irem(r,p) fi
od;
s*r
end proc,
`Mul/2arg`:= proc(E, J::(name=range), p::posint)
description
"Just like two-argument `mul`, but in modular arithmetic. "
"The modulus needn't be prime. Attempts hfloat mode when p < 2^26."
;
local
r, k, s, m:= lhs(rhs(J)), n:= rhs(rhs(J)), M,
#This is an iterated subs--not an ordinary multiple subs--used like `unapply`:
e:= subs(_E= E, lhs(J)= _K, _K-> _E)
;
if p < Mhf then
#It's worth trying to compute exactly in hfloats because it's much faster.
try return trunc(evalhf(`Mul/2arg/evalhf`(e, m, n, p, M))) mod p
catch: userinfo(1, Mul, "Evalhf failed. Will use endogenous arithmetic instead.")
end try
fi;
r:= 1; s:= true; M:= max(p, iquo(Mgmp, p));
for k from m to n while r <> 0 do
r:= r*e(k);
if r < 0 then r:= -r; s:= not s fi;
if r >= M then r:= irem(r,p) fi
od;
`if`(s, r, -r) mod p
end proc,
Binomial:= proc(N::nonnegint, K::nonnegint, p::prime)
option remember;
local r, n, k:= min(K, N-K), i, num, M:= max(p, iquo(Mgmp, p));
if p > k and k >= 0 then #Straight-forward method for large p
num:= Mul(N-k+1..N, p);
return `if`(num=0, 0, num/':-Mul'(1..k) mod p)
fi;
#Lucas's method for small p, with callbacks to method above for the
#sub-binomials:
r:= `if`(K > N, 0, 1); n:= N;
while k <> 0 and r <> 0 do
r:= r*thisproc(irem(n,p,'n'), irem(k,p,'k'), p);
if r >= M then r:= irem(r,p) fi
od;
r mod p
end proc,
ModuleLoad:= proc()
local P, N;
#Set upper bounds of computation modes:
Init_Ms();
#Insert them into relevant procs:
`Mul/2arg/evalhf`:= subs(_M= Mhf, eval(`Mul/2arg/evalhf`));
#Compile compilable procedure:
Compiler:-Compile(`Mul/range/compiled`, 'warnings');
#Protect global `mod` names and assign `mod/...` procs:
for P in ModMethods do
N:= cat(`mod/`, P);
if not N::protected then assign(N= P) fi;
protect~([N, cat(``, P)])
od;
return
end proc,
`;`;
ModuleLoad()
end module:
Binomial_mod_p.mw
I need to add some compiled code to this, akin to what VV did. When I do, I'll just edit this Reply.
[Edit: I've now added the compiled code.]
[Edit 9-Nov-2018: The code and worksheet have been updated.]
