@vv The two methods (VV's and mine) should be combined into a single procedure, so here it is. I think that p > k is sufficient for your method, right?

`mod/Mul`:= proc(E, J::(name=range), p::posint)
description "Just like `mul`, but in modular arithmetic. The modulus needn't be prime.";
local r:= 1, k, e:= subs(lhs(J)= k, E), s:= true;
   for k from lhs(rhs(J)) to rhs(rhs(J)) while r > 0 do
      r:= r*eval(e,2);
      if r < 0 then r:= -r; s:= not s fi;
      if r >= p then r:= irem(r,p) fi
   od;
   `if`(s, r, -r) mod p
end proc:

`mod/Binomial`:= proc(N::nonnegint, K::nonnegint, p::prime)
description "Computes binomial(N,K) mod p, p::prime, by Lucas's theorem";
option
   author= "Carl Love <carl.j.love@gmail.com> 31-Oct-2018",
   reference= "https://en.wikipedia.org/wiki/Lucas%27s_theorem"
; 
local r:= `if`(K > N, 0, 1), n:= N, k:= min(K, N-K), i;
   if r = 0 then return 0 fi;
   #Straight-forward method for large p: 
   if p > k then return Mul(n-i, i= 0..k-1)/Mul(i, i= 1..k) mod p fi; 
   while k > 0 and r > 0 do #Lucas's method for small p:
      r:= r*irem(binomial(irem(n,p,'n'), irem(k,p,'k')), p) 
   od;
   r mod p
end proc:

Your example:

(n,k,p):= (8*10^5, 5*10^5, nextprime(8*10^5)):
CodeTools:-Usage(Binomial(n,k) mod p);
memory used=84.86MiB, alloc change=6.01MiB, cpu time=781.00ms, 
real time=780.00ms, gc time=31.25ms
                             494383

@acer I agree with what you did: It's a good reason to edit a Question. And thank you for being a conscientious Moderator. There should still be a way to do it without changing the position, or to reset the position once it's done.

It's likely possible, but you need to clean it up. I mean clarify. There's no p in your equations. What is G? Will you be able to specify p and q numerically at the time that you solve it?

This is definitely not handled by the command solve, which operates strictly over the complex numbers. (Even trying to restrict solve to real numbers is often problematic.) But there are several integer solvers, such as chrem, msolve, and isolve. Look up their help by typing 

?chrem

etc., at a command prompt.

@Rohith Since you only deal with polynomials, there may be a symbolic solution available for you. See ?solve,parametric. Using your posted example:

f:=2*a*b-3*a-2*b-c:  # 15 <= f <= 30

Pick any two variables to get piecewise solutions broken down by the third. For example, do

solve({f > 30}, {a,b}, parametric)

This'll tell you (among a great many other things) that c=33 is a critical value upon which the branching of solutions happens. Apply this recursively:

solve({f > 30, c < 33}, {a}, parametric)

We thus learn that b=3/2 is critical. We are building a tree of all valid combination of a, b, c, f. This is hard symbolic work, but it's a start. In particular, deconstructing any piecewise programmatically is annoying; deconstructing those spit out by solve(..., ..., parametric) is brutal^[*1].

The parametric option to solve is only available for sets of polynomial inequalities and/or equations.

^[*1]It would help a little if there were an option akin to the ifactor / ifactors dichotomy: one version being prefrerred for visual display and one for programmatic uses.
 

@Joe Riel The call pobj:-eval(..., pobj) can be (and I think should be) shortened to pobj:-eval(...) by giving up eval's static status and changing it to 

eval:= e-> :-eval(e, [eqs])

I don't mind needing to type pobj:-eval(...) instead of eval(pobj, ...) (or some such). I am quite skeptical about calling module procedures without their module prefix anyway because I think that it makes code difficult to read. The only time that I find it acceptable is for binary and unary operators.

I think that we're more than halfway towards having a generic container object for managing parameters!

@acer Yeah, I think that your check is a good idea, one random point being sufficient. I too am hesitant about trig integrals over more than one period or polar integrals over more than one revolution, even when doing them by hand.

@acer Why did you feel the need for this particular problem to verify the results numerically? Was there some step in the process that seemed not entirely trustworthy to you?

@Joe Riel In my working copy, I actually had 

params:= Record(a::algebraic= 'a', b::algebraic= 'b', c::algebraic= 'c', e::algebraic= 'e')

specifically to deal with the issue of the value of an unassigned Record member including its type specification. (The fact that it does include that seems like a kernel bug to me.) I forgot to include that change in my posted code.

Also, it was my intention that passing equations to the ModuleApply whose lhs didn't match any member of params would be an error, although I can certainly see why some other programmers might not want it that way.

Thanks for the object code. Yes, the dual nature of eval is quite a nuisance, and not just in this case. It behaves badly in ways that only a builtin can. I just tell people "There are two builtin commands named eval. They are completely different. They are both documented at ?eval. To add to the confusion, one of them is sometimes called 'two-argument eval' even though the other can be used with two arguments." Sigh. It's quite a shame that the newer eval wasn't named evalat, but it'd be impossible to fix now.
 

@tomleslie I think that it's the duty of a responsible Answer Giver to warn that using cat like this creates global names, regardless of the status of the originals. But, also, a responsible OP should read the help pages of the commands recommended in an Answer.

@brian bovril To me, it was clear from the start that you wanted to index the names, not create new symbols, but I see how your meaning could've possibly been misinterpreted.

@Joe Riel Easier to find is the library function index. `?[]`(N, J) is equivalent to index(N, J[]) (even if J is the empty list).

@Joe Riel Yes, and the 43 capitalized semi-inert functions described at ?mod, plus &^, are implemented in Maple as `mod/...`.

Did you see my Reply to your Answer of my thismodule Question?

This is called the inverse image or preimage problem (the two terms mean the same thing). It can be extremely difficult to get even an approximate solution when there are continuous variables and the domain has more dimensions than the codomain, as in your case. But it's a very important problem with applications in 3D imaging, for example, How does one construct a 3D model of an object from a series of its 2D projections?. (Most of our brains manage this computation in real time---we don't quite know how.) The example that you've given is significantly simpler with a 1D codomain and real-interval range/projection.

I notice that your f is quadratic. Is it always so? It'll be easier to get a solution for that case.

@Carl Love If one changes the module as I suggest above, then it becomes useful for it to export a procedure for evaluation of expressions wrt the record params. This is quite simple: Declare a local (preferably immediately after the Record's declaration):

Params:= [exports(params)],

and declare an export Eval (or other appropriate name):

Eval:= e-> eval(e, Params=~ map2(index, params, Params)).

Although the above is easy to code, it's mind-boggling to me why evaluation wrt a module or record is not part of the regular eval command. It seems like such a fundamental operation. The use command is not an adequate substitute because it only works at the top level (i.e., on a command line).

@Carl Love It should be noted that when a single subsop command is used with multiple indices (as in this case), and some of the substitutions cause the target object's operand numbers (aka indices) to change (as in this case, where NULL is being substituted into a list), then the indices are all specified as their original values (i.e., the values that pertained at that start of the command). Since this is the more-convenient way to do it, and the way the you were doing it anyway, I didn't mention it earlier. However, some Maple commands do similar operations recursively or have modes for both recursive and non-recursive operation (such as subs and eval).