@digerdiga I didn't say that it was only true for 0 < t < Pi; I just said that it was true for those t, which left open the situation for other t. I was just providing an example of a single thing that you could change in your command that would make the simplification valid. Indeed, the interval could be stretched to -Pi < t <= Pi.

@vv You can replace simplify(convert(..., exp)) with combine(...).

@brian bovril Yes, your iteration is correct. But I don't want you to get the wrong impression that the reason to use x=y on the first iteration is because it allows the simplification of ln(1) = 0. That's just a side benefit. The reason is that x+ln(x) is nearly the identity function on your interval of interest, so x=y is already close to the solution.

Also, three iterations are required to be certain of getting the three decimals places that you asked for over the domain that you asked for. You show two iterations. Two iterations will get three decimal places at the lower end of the domain, but not the higher end.

@Rohith The commands select, remove, and selectremove only search the top-level operands of an expression. This job requires indets. I'll give more details later.

[In the following, W(x) (without boldface) refers to the mathematical Lambert W function rather than to my procedure W. On MaplePrimes, I always use boldface for inline references to literal Maple code.]

As VV has astutely observed, my first "fold" W(x,x) makes x the starting point for the iteration, and this convenient choice eliminates the ln from the first iteration. That convenience was not the reason that I chose this.

Since W(x) is the inverse of x*exp(x), it's fairly close to ln(x). Both ln(x) and W(x) are steep for 0 < x < 1. Steepness causes difficulty for Newton's method. But the function G(x) = W(exp(x)) is nearly the identity function---no steepness. It's inverse is x+ln(x) (as can be verified with solve), and thus the solution x of x+ln(x)=ln(y) is W(y). So, it's this equation rather than the more usual x*exp(x) = y that I solve with Newton's and/or Householder's method. Because of the closeness to the identity function, x0 = y is a great choice for the initial guess.

x + ln(x) = ln(y)  =>  x+ln(x/y) = 0. I define g:= (x,y)-> x+ln(x/y). The Householder method of order 1 is the ordinary Newton's method. The Householder method of any order is given by this simple procedure:

Householder:= proc(F, ord::posint)
local x, y;
   unapply(x+ord*D[1$ord-1](1/F)(x,y)/D[1$ord](1/F)(x,y), (x,y))
end proc:    

where F is a two-parameter procedure such that F(x,y) = 0 is the equation that we want to solve for x. Doing Householder(g, 3) puts a reasonably sized expression on my screen. Then I meticulously hand-optimized it by removing common subexpressions, subfactoring, etc. Why did I choose order 3? Because it seemed to have a nice balance between the simplicity of the iterative function and the number of iterations required for the desired accuracy. I haven't  checked, but it may minimize the total operation count.
 

@acer It's fine by me. Duplicate answers are unavoidable. You know what they say about great minds. I was surprised that your initial Answer didn't mention eval.

If a function appears in the expression without any coefficient, would you like its coefficient listed as 1, or would you rather that it not be listed at all?

In your point 2, you have made an unsubstantiated claim of a serious error in Maple. You have had all day to post some evidence, but you have not. All claims of errors, even if unsubstantiated, cause damage to the reputation of the software and the company that makes it. As such, it is unethical to make those claims, just like it would be unethical to spread unsubstantiated claims that an associate of yours had a contagious disease.

So, I give you 12 more hours to post some evidence. If you won't I'll delete this thread, unless some other respected moderators can tell me some good reason not to.

@phbmc Yes, this is a 2D-Input problem that I've encountered before. It can't handle parameters whose type specifications refer to other parameters. So, if you make the typespec for L simply L::list, it'll work in 2D. You may, if you wish, put an error check in the body of the procedure:

     if nops(L) < 2*c then error "List not long enough." fi;

Maple, oddly enough, doesn't care if you pass unneeded arguments to a procedure. It'll just ignore them (unless the code specifically looks for them). So, you can just pass in m as a third argument (but I think that it's a weird thing to do). But, I don't see how it can be done without passing the list of items to be permuted. Does that make sense to you?

If you have a collection of related procedures, they can be (and probably should be) collected together into a module. That can easily depend on m, and I wouldn't even have any stylistic objections.

@zarara I know that it seems too easy to be true, but the limit command that VV shows does it all in this case, which is a super-easy problem. That is, if you're willing to accept the output of a computer program as "proof"...which I am. The big-O-ness is equivalent to that limit not being infinity, when that limit exists. 

By "snap" do you mean to click on, as with a mouse or other pointing device?

[This Reply is not related to binomial-mod-prime computations, and is intended only for expert-level module writers who write some numerical or semi-numerical code. It relates to a technique that I used in the module above.]

I just developed a coding technique---shown in the module above---that may be novel. It allows lexical variables (in this case module locals) to be used in evalhf'able and compilable procedures. The evalf'able and compilable procedures are defined containing a dummy variable _M. The values of _M are potentially machine dependent (at least theoretically), so they are computed in the ModuleLoad and subs into the procedures. The same thing could be done for any values that are not known until runtime.

I am saddened by the lack of enthusiasm for this Post. I believe that the code above is the only implementation in Maple of any exact^[*1] statistical test of independence in a contingency table (the same property that is tested by the inexact Statistics:-ChiSquareIndependenceTest). Some modern statisticians recommend that given that we have the necessary mathematical software, that exact tests always be used for small n (the sum of the entries in the table), say n < 1000. Furthermore the code above works for tables of any number of cells (i.e., the number of entries in the matrix) whereas some other programs that do Fisher's exact test are limited to two rows, two columns, or both.

^[*1]"Exact" in this context refers to the fact that the computations are done with the relevant discrete distributions, the multinomial (generalized hypergeometric), rather than with a continuous distribution such as Chi-squared, which is only asymptotically equivalent to the discrete distributions as n approaches infinity. "Exact" does not refer to the fact that the computations are done in exact rational-number arithmetic; however, the code above does that also. 

@vv I understand now why there was that factor-of-4 time difference between our codes in the last test that you showed: I was trying to do too much with just one variety of Mul, which handled arbitrary expressions akin to mul. I improved it by adding a separate Mul for ranges of integers, such as factorials, and I made it compiled for the appropriate cases. Please try it out.

I may have used the word "recursive" incorrectly earlier, leading to some confusion on your part. I simply meant that the procedure Binomial makes a single call back to itself to handle the smaller binomials. I didn't mean that the Lucas code was being called again.

Because of its abstractness and generality, my previous modular Mul command wasn't efficient for its most-common use (or most-common use relevant to this thread): the product of a range of integers (such as a factorial). So, I've now implemented Mul(m::integer..n::integer) mod p::posint as syntax additional to the Mul(e::algebraic, k::name= m::integer..n::integer) mod p::posint that I'd already implemented. This brings my code into closer alignment with VV's Fac procedure.

You don't need to use or understand any Mul syntax to use this code for its primary purpose: computing binomial(n,k) modulo a prime. It's all "under the hood".  You simply use the syntax

Binomial(n,k) mod p.

The code is now (code below) in the form of a module. Simply entering the code into a Maple session or reading it in from a file will install the `mod/...` procedures. If the module is stored in a library, then simply invoking Binomial_mod_p() will install them.

Here is the code. For convenience, it also exists in the Code Edit Region of the attached worksheet.

Binomial_mod_p:= module()
description "Some tools to speed up the computation of binomial(n,k) mod p, p::prime";
option
      date= "9-Nov-2018",
       authors= (
          "VV from MaplePrimes  <...>",
          "Carl Love <carl.j.love@gmail.com> (Corresponding)"
    ),     
    reference= "https://en.wikipedia.org/wiki/Lucas%27s_theorem" 
;
################################################################################################################
# This is a module of tools whose primary purpose is to speed up the computation of binomial(n,k) mod p, for   #
# prime p. This code interacts with `mod` through the `mod/...` global hooks (an ancient pre-module            #
# interface). So, one simply invokes Binomial(n,k) mod p, where the Binomial is an inert function (see ?mod).  #  
#                                                                                                              #
# The speed of the computation, and the choice of algorithm, are based on the relative sizes of n, k, and p.   #
# all cases where p <= binomial(n,k)*k!, this code makes some significant attempt to use modular arithmetic to #
# improve upon the speed of the naive command binomial(n,k) mod p.                                             #
#                                                                                                              #
# This module has no exports. It may "installed" any of these ways:                                            #
#     1. Simply copy-and-paste this code into an execution group and execute;                                   #
#    2. Use the `read` command to import it from a file;                                                       #
#    3. If it's stored in a library, execute the module's name: Binomial_mod_p().                              #
# Any of these techniques will assign the necessary global `mod/...` procedures.                               #
#                                                                                                              #
# The module has a few userinfo statements, all keyed to infolevel[Binomial]:= 1.                              #
################################################################################################################
local
    #The methods implemented as `mod/...` procedures by this module are Mul and Binomial. If methods are added,
    #then their names should local procedures in this module named X such that the corresponding `mod/...` 
    #procedure will be named `mod/X` (and thus the usage will be X(...) mod p). The names should be added to
    #the list immediately below. Everything else will then happen automatically. 
    ModMethods:= [Mul, Binomial],    

    #The ModuleApply only exists so that invoking Binomial_mod_p() will force the ModuleLoad to
    #run if the module is stored in a library.
    ModuleApply:= proc($)
        userinfo(1, Binomial, "The `mod/...` procedures have been protected and installed.");
        return
    end proc,

    ModuleUnload:= proc()
        userinfo(1, Binomial, "The module Binomial_mod_p is being garbage collected.");
        return
    end proc,
        
    #Some constants specifying the largest "convenient" integer for various modes of computation:
    Mhf, Mgmp, Mcomp,
    #...and a proc to initialize them:
    Init_Ms:= proc()
        #In hfloat mode, we try to keep the mantissa half full so that products don't overflow it: 
            Mhf:= isqrt(Scale2(1, trunc(evalhf(DBL_MANT_DIG))-1));  
        #In GMP mode, we try to keep products to 1-word integers.
         Mgmp:= kernelopts(maximmediate);
         #Likewise for compiled mode, but the wordsize is slightly different.
         Mcomp:= Scale2(1, kernelopts(wordsize)-1) - 1;
        return
    end proc,

     ###########################################################
    # The forms of Mul currently implemented are
    #    - Mul(m..n): 
    #        a range of integers;
    #    - Mul(e, k= m..n):
    #         an expression e (likely depending on name k), with k iterated over a range of integers. 
       Mul:= proc()
    option remember;
    description "like mul, but in modular arithmetic";
           if nargs=2 then 
               if args[1]::range then `Mul/range`(args)
            else `Mul/container`(args)
             fi
           elif nargs=3 then `Mul/2arg`(args)
           else error "Expecting 1 or 2 arguments, but received", nargs-1
           fi
    end proc,

    `Mul/range/compiled`:=  
    proc(m::integer, n::integer, p::posint)::integer;
    option autocompile;
    local k::integer, r::nonnegint:= 1;
        for k from m to n do
            r:= r*k mod p;
            if r = 0 then break fi
        od;
        r
    end proc,

    `Mul/range`:= proc(R::range(integer), p::posint)
    description
          "Computes the modular product of a range of integers (just like 1-argument mul, but modular)."
           " The modulus needn't be prime. Attempts compiled mode for word-sized p."
    ;
    local r, k, m:= lhs(R), n:= rhs(R), M, MR:= max(1, CopySign(m,1), CopySign(n,1)); 
          if p < Mcomp then #Use compiled code.
              try return `Mul/range/compiled`(m, n, p) mod p
              catch: userinfo(1, Binomial, "Compiled code failed. Will use endogenous arithmetic instead.")
              end try
        fi;
        r:= 1; M:= max(p, iquo(Mgmp, max(1, MR)));
        for k from `if`(m>=0, m, -n) to `if`(m>=0, n, -m) while r <> 0 do
            r:= r*k; 
            if r >= M then r:= irem(r, p) fi
        od;
        `if`(m>0 or (k-m)::even, r, -r) mod p
    end proc,

    `Mul/2arg/evalhf`:= proc(e, m, n, p)
     local k, r:= 1, s:= 1, M:= _M, 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,
            
    `Mul/2arg`:= 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."
    ;
    local 
           r, k, s, m:= lhs(rhs(J)), n:= rhs(rhs(J)), M,
           #This is an iterated subs--not an ordinary multiple subs--used like `unapply`: 
           e:= subs(_E= E, lhs(J)= _K, _K-> _E)
    ;
          if p < Mhf then 
              #It's worth trying to compute exactly in hfloats because it's much faster.
             try return trunc(evalhf(`Mul/2arg/evalhf`(e, m, n, p, M))) mod p
              catch: userinfo(1, Mul, "Evalhf failed. Will use endogenous arithmetic instead.")
             end try  
           fi;
           r:= 1;  s:= true;  M:= max(p, iquo(Mgmp, 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,

    Binomial:= proc(N::nonnegint, K::nonnegint, p::prime)
    option remember;
    local r, n, k:= min(K, N-K), i, num, M:= max(p, iquo(Mgmp, p));
        if p > k and k >= 0 then #Straight-forward method for large p 
             num:= Mul(N-k+1..N, p);
             return `if`(num=0, 0, num/':-Mul'(1..k) mod p)
        fi;
           #Lucas's method for small p, with 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*thisproc(irem(n,p,'n'), irem(k,p,'k'), p);
              if r >= M then r:= irem(r,p) fi 
           od;
           r mod p
    end proc,    
    
    ModuleLoad:= proc()
    local P, N;
        #Set upper bounds of computation modes:
        Init_Ms();
        #Insert them into relevant procs:
        `Mul/2arg/evalhf`:= subs(_M= Mhf, eval(`Mul/2arg/evalhf`));
        #Compile compilable procedure:
        Compiler:-Compile(`Mul/range/compiled`, 'warnings'); 
        #Protect global `mod` names and assign `mod/...` procs:
        for P in ModMethods do
            N:= cat(`mod/`, P);
            if not N::protected then assign(N= P) fi;
             protect~([N, cat(``, P)])
        od;
        return
    end proc,  
`;`;
    ModuleLoad()
end module:
   

Binomial_mod_p.mw

I need to add some compiled code to this, akin to what VV did. When I do, I'll just edit this Reply.

[Edit: I've now added the compiled code.]

[Edit 9-Nov-2018: The code and worksheet have been updated.]