@Carl Love 

You have invested a lot of work and have transformed the problem into a larger project. Congratulations!
I have some remarks/questions.

- If you use the compiler, why do you keep the evalhf part? evalhf is slower and works for smaller integers (in Maple 64 bits).
BTW, the compiler does not accept max with >2 arguments, you must change this.
- I think that the compiled part shoud be used whenever p<pmax, without try/catch. I also don't thing that recursion is useful (I have mentioned this in a previous comment).
- For
#p:=nextprime(8*10^7); n:=p*2-3;  k:=floor(3*p/2);  # answer=50000014
p := 80000023;n := 160000043;k := 120000034;

Binomial needs minutes instead of seconds and a lot of gc. Is this because of recursion? Does the compiler work only once?

In maths we must accept the definitions. Unfortunately it may happen that our prefered vesion is changed in another source. In Maple  the CDF is F(x) = P(X<=x);  I would prefer  "<", but that's it.
I stopped reading the worksheet when I saw that the definitions where not included or discussed.

You said that the quantile results are not correct for the discrete uniform distribution. But it conforms to the definition given in the help page (and widely accepted). So, où est le reproche?

@rlopez 

Why challenging? I'd say that the problem is simply incorrectly or/and incompletely formulated.

@Mariusz Iwaniuk 

Actually the type AndProp(integer,RealRange(0,0))  is  0  (0 is a type itself).
So, _Z9  should have been replaced with 0, but probably it's too much to ask for this.
 

@Carl Love 

For p<<k the Lucas method is of course better than BIN which I have posted just for fun. It was not supposed to be integrated with Lucas because here only the case p<k is needed (as in "methods"). From the same reason the recursive approach is not useful because for n,k < p Lucas does nothing.
For "small" k,n,p  the evalhf idea seems to be effective.
Let's compare:

BIN := proc(N, K::nonnegint, p::prime)  # combined with Carl's
local r:= 1, n:= N, k:= min(K, N-K), Fac,Bin;

Fac:=proc(n,m)
  local r:=1,k;
  for k from n to m do r := irem(r*k,p) od
end;

Bin:=proc(N,K)  # 0 <=K < p, 0 <= N
local k:=min(K,N-K);
`if`(k<0, 0, Fac(N-k+1,N)/Fac(1,k) mod p)
end;

while k > 0 and r > 0 do 
   r:= irem(r*Bin(irem(n,p,'n'), irem(k,p,'k')), p) 
od;
r
end:

p:=nextprime(8*10^7): n:=p*2-1:  k:=floor(p/2):
CodeTools:-Usage(Binomial(n,k) mod p);
memory used=8.35GiB, alloc change=0 bytes, cpu time=4.65m, real time=4.65m, gc time=2.55m
                            80000022
CodeTools:-Usage(BIN(n,k,p));
memory used=1.79GiB, alloc change=0 bytes, cpu time=65.55s, real time=65.54s, gc time=30.65s
                            80000022

@mmcdara 

If you do not have any information about F, the problem is too vague and a solution is practically impossible.
- if F is far for being linear, using the convexity is not reasonable. Probably a better approach will be to find the boundary.
- if F is not injective it will be useful to find a minimal and "nice" X i.e. a convenient selection of the multivalued map F^(-1). [BTW, even the existence of a continuous selection is also a hard problem].
Onece again, it depends heavily on a theoretical study of F and also on your future intention with A (or X).

@vv 

I have mentioned that binomial(N,K) mod p can be adapted for K>=p too. Here is the stand-alone procedure

BIN:=proc(N::nonnegint,K::nonnegint,p::prime)
local k:=min(K,N-K), i,j, r1:=1, r2:=1, e:=0;
if k<0 then return 0 fi;
for i from N-k+1 to N do 
    j:=i; while irem(j,p)=0 do j:=j/p;e:=e+1 od;
    r1 := irem(r1*j,p) od;
for i from 1 to k do 
    j:=i; while irem(j,p)=0 do j:=j/p;e:=e-1 od;
    r2 := irem(r2*j,p) od;
if e>0 then 0 else  r1/r2 mod p fi
end:

It is rather fast and it does not use Lucas' theorem.
However, for p << k, using this theorem can increase the speed, provided that   "smaller" binomials in the product are computed mod p.
So, the best method seems to be the combination presented in the comment "methods".
 

@Carl Love 

The two methods should be deeper combined. E.g. for
p:=nextprime(8*10^5): n:=p*2-1: k:=floor(p/2);
the direct method is still much faster.

Example of combinaton

Fac:=proc(n,m,p)
local r:=1,k;
for k from n to m do r := irem(r*k,p) od
end:

Bin:=proc(N::nonnegint,K::nonnegint,p::prime)  
local k:=min(K,N-K);
`if`(k<0, 0, Fac(N-k+1,N,p)/Fac(1,k,p) mod p)
end:

Bin1 := proc(N::nonnegint, K::nonnegint, p::prime)  # combined with Carl's
local r:= 1, n:= N, k:= min(K, N-K);
   while k > 0 and r > 0 do 
      r:= irem(r*Bin(irem(n,p,'n'), irem(k,p,'k'),p), p) 
   od;
   r
end proc:

N.B. I prefer to write the procedures non-optimized and not integrated in `mod/...` (for better understanding the algorithms).

@Rohith 

It's elementary to scale the variables. If x is in a..b with step h, simply take x[i] = a + i*h, i=0,1,...
This way you can use arrays or Tom's code.

@Carl Love 

Unfortunately it seems that  in shadebetween  one must remember to reverse the order of the ranges compared with int.

@Mariusz Iwaniuk 

You have d1 instead of d3 for the 3rd sphere.

@brian bovril 

You have found a bug.

restart;
e:=5:
seq(e,e=[sin(u*Pi),2]);# assuming u::integer;

works. But after deleting ;#  it gives

Error, (in assuming) when calling '`one of {seq, sin}`'. Received: 'illegal use of an object as a name'

@mehdi jafari 

R gives the parametrization of the surface in spherical coordinates, R = R(u,v). Recall that the Jacobian is r^2*sin(v).
So, the volume equals
int( r^2*sin(v),  r=0..R, u=0..2*Pi, v=0..Pi);
(or equivalently, the double integral above).
 

@brian bovril 

It seems that you have installed a maple init file which does something wrong.
The code works in Maple 2015, 2017 and 2018.

@Carl Love 

Thank you Carl. I think that the situation is clear now.