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

@Carl Love 

Why are you mentioning only undefined symbols? E.g. Beta is defined and Beta('''f''', A)(x)  acts the same.
The problem is that evalapply is builtin and its action on most functions is not documented. Probably it matters whether the procedure is associated to an operator (such as `[]` and `+`); that's why I have mentioned `&+`.
Note that it's possible to define `evalapply/Beta`  but `evalapply/[]` is ignored.

@Carl Love 

You say "this behavior is a natural and logical result of the overall design".
However, the procedures `[]`, `+`   behave differently than  `&+`, `or`  wrt evaluation of  '''f'''
 

['''f''', A](x);
`[]`('''f''', A)(x);
`+` ('''f''', A)(x);
`&+`('''f''', A)(x);
`or`('''f''', A)(x);
    [('f')(x), A(x)]
    [('f')(x), A(x)]
    ('f')(x)+A(x)
    (`&+`(''f'', A))(x)
    f(x) or A(x)

@Christian Wolinski 

I don't quite understand this type of comments. It would be neater to see the command line followed by its result (given by Maple) and your expected result (+ maybe some explanations). Otherwise it's just some obfuscated code, and the reader is supposed to decrypt it.

On the other side, do you want to program something in Maple and the evaluation of  [...](...)  makes this impossible? In this case please present an example.