@Carl Love 

I am not able to provide such an example. A good candidate should be
isolve((y^2-x^3-1)*(x-2*y)*(y-2*x) = 0);
but Maple does not solve y^2-x^3-1 = 0 and however the flag is not set; it should be clear (for Maple and for us too) that not all solutions have been found e.g. x=2,y=3 [actually this is the only one which is missing].
The reason I asked the question is that I don't know why _SolutionsMayBeLost   has disappeared. It would be nice to have it back.

@Carl Love 

I don't think you are right about _SolutionsMayBeLost . solve has no reason to set it because the system is very simple.
And as I have said, the flag was set in earlier versions (and the answer was the same).

@boblas 

I have already shown this:

prod( wdeg([1,1,0,1,1,0,etc], [x,y,z,t,u,v,etc]), plex(x,y,z,t,u,v,etc) );

If you want to test it, please come with a concrete polynomial P (or a list of monomials) and provide also the sorted version of P.

@Carl Love 

If the starting point is changed from x to  `if`(x<2,x/2, ln(x/ln(x)))  then 2 iterations are enough for any x>0 with accuracy 4e-10 for evalhf.  And actually the accuracy is better; for evalf[30] ==> 1e-17.

W(W(`if`(x<2,x/2, ln(x/ln(x))),x),x):

But of course this is not important because an iteration of W is probably not much more expensive than ln(ln(x)).

Edit.  W(W(x/2,x),x)  has accuracy 1e-7 in the interval 0.3 .. 10

@Carl Love 

Very nice, vote up!
So, you took x as starting point for iterations. Of course there are better ones, but it has the advantage that the first iteration does not contain ln.

@brian bovril 

I did not derive it; it's in

D.A. Barry et al.,  Analytical approximations for real values of the Lambert W-function,
Math. Comput. Simul. 53(2000), 95-103.

It has not a great accuracy but it's simple and is valid for x>0.

limit(LambertW(x)-W, x=infinity) = 0.051...

@zarara 

We have just proved it.
Note that Maple is able to obtain the answer without knowing it:
MultiSeries:-asympt(n^2/2 + 2*n - 1, n, 0);
        O(n^2)

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