@acer No doubt your hypergeometric expression is equivalent to my regularized incomplete beta function in integral form. I lost my kernel connection about a dozen times (mostly on smaller p) while producing that post, so I guess that the fault is in the numeric integration rather than NextZero. Did you have any such trouble?

@vv I am extremely skeptical that your direct modular method for p <= k (posted immediately above) could possibly be faster than using Lucas's combined with the direct modular method for the smaller binomials. Please post an example that shows it.

@jga You wrote:

• To be honest I can't follow the discussion that's going on in here.

There are two algorithms being discussed here: Lucas's and what VV may have inadvertently dubbed "the direct method". Both of these are implemented in modular arithmetic. Both are much much faster than the naive direct method binomial(n,k) mod p. The only issue is when is the optimal time to automatically switch from one method to the other. You don't need to understand the inner workings of this in order to use the code; it's done automatically.

• I can tell you that in my algorithm there's no relation between N and K with p.

Even if they're completely independent, you must have some idea of their individual ranges or distributions.

• I'll try the integrated version and see how does it compare to the direct method in my case. 

You'd just be wasting your time if by "direct method" you mean binomial(n,k) mod p. Use the integrated version that I just posted. It's definitely faster.

• After this I still need to find the kernel (or equivalently, the rank) of them in prime characteristic.

The command LinearAlgebra:-Modular:-Rank is super fast for this, especially for p < 2^25.

@vv I agree that they should be recursively combined, and I've done it below. I've also implemented an automatic diversion to evalhf for the "direct method" when the computations can thus be done exactly, because it's much faster. (They can be done exactly when the magnitude of the result is less than 2^53.) Any deviation from exactness causes the code to revert to native GMP.

`mod/Mul`:= 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."
;
option author= "Carl Love <carl.j.love@gmail.com> 2-Nov-2018";
local 
   r, k, s, m:= lhs(rhs(J)), n:= rhs(rhs(J)),
   #This a is an iterated subs, not an ordinary multiple subs: 
   e:= subs(_E= E, lhs(J)= _K, _K-> _E),
   #In hfloat mode, we try to keep the mantissa half full
   #so that products don't overflow it: 
   M:= Scale2(1, iquo(trunc(evalhf(DBL_MANT_DIG)), 2)) #i.e., M = 2^26 
;
   if p < M then 
      #It's worth trying to compute exactly in hfloats because it's much faster.
      try
         return 
            trunc(evalhf(proc(e, m, n, p, M)
            local 
               k, r:= 1, s:= 1, 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
            (e, m, n, p, M)
            )) mod p
      catch:
         userinfo(1, Mul, "Evalhf failed. Will use GMP instead.")
      end try  
   fi;
   r:= 1;  s:= true;  M:= max(p, iquo(kernelopts(maximmediate), 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:

`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> 2-Nov-2018",
   reference= "https://en.wikipedia.org/wiki/Lucas%27s_theorem"
; 
local r, n, k:= min(K, N-K), i, num, M:= max(p, iquo(kernelopts(maximmediate), p));
   if p > k then #Straight-forward method for large p 
      num:= Mul(i, i= N-k+1..N) mod p;
      return `if`(num=0, 0, num/Mul(i, i= 1..k) mod p)
   fi;
   #Lucas's method for small p, with recursive 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*irem(thisproc(irem(n,p,'n'), irem(k,p,'k'), p), p);
      if r >= M then r:= irem(r,p) fi 
   od;
   r mod p
end proc:

Here's your most-recent example:

p:=nextprime(8*10^5): n:=p*2-1:  k:=floor(p/2):
CodeTools:-Usage(Binomial(n,k) mod p);
memory used=4.62KiB, alloc change=0 bytes, cpu time=312.00ms, real time=317.00ms, gc time=0ns
                             800010

@mmcdara According to the Wikipedia article "Negative binomial distribution", it is common (perhaps even standard) to use non-integer values for the "number of failures" parameter r. There it is said about the associated nomenclature:

• The Pascal distribution (after Blaise Pascal) and Polya distribution (for George Pólya) are special cases of the negative binomial distribution. A convention among engineers, climatologists, and others is to use "negative binomial" or "Pascal" for the case of an integer-valued stopping-time parameter r, and use "Polya" for the real-valued case.

So it seems that the use of non-integer r reflects how the distribution is used in actual practice and has nothing to do with internal representations of factorials as Gamma functions.

@Rouben Rostamian Ah, I didn't "see" that it was possible to eliminate varphi from the constraint. Indeed the constraint is linear wrt varphi. As I have always been suspicious of the implicitdiff command, someone should check whether the results are really equivalent. If they aren't, then I have far more confidence that your results are the correct ones.

@Joe Riel Are static procedures more time efficient on a per-use basis? Or are they just more efficient (timewise and memorywise) when the object is created? 

@acer Can you give an example of what you mean? (It is frustrating, yet understandable, that thismodule and static don't mix.)

@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.