Use an embedded if statement rather than `if`, like this:

str:= "A";
x:= 10;
str:= sprintf(
    "%s it was%s 10", str, if x++=10 then "" else x:= 9; " not" fi
); 

The benefit of an embedded if statement is that any number of statements can appear in its then and else clauses. It's only the value of the final of these statements that becomes the return value of the embedded statement.

Your computation can be done by very simple arithmetic in compiled code. For the n=8 case, that avoids (2^8)^3 ~ 16 million calls to Bits:-Split, (2^8)^2 + (2^8)^3 calls to ListTools:-Occurences, and the creation and garbage collection of the same number of lists. My code below does the n=8 case in 17 seconds total on one processor. Further gains can be obtained by

• Writing the code directly in C so that you can use whole-word bit operators.
• This code can be parallelized with Grid or Threads (Threads is usually faster),
• Perhaps some speedup can be achieved by using 2-byte integers directly in C.

I made a great number of changes with the arithmetic. Please verify the accuracy of this by comparing against known-correct code using the same random vector M.

restart
:
BuildA_M:= proc(
    A::Array(datatype= integer[4]), M::Array(datatype= integer[4]),
    r::integer[4]
)::integer[4];
local 
    i0::integer[4], i::integer[4], j0::integer[4], j::integer[4],
    m::integer[4], t0::integer[4], t::integer[4],
    p::integer[4], S::integer[4]
;
    for i0 from 0 to r-1 do
        for j0 from 0 to r-1 do
            S:= 0;
            for t0 from 0 to r-1 do
                p:= 1; i:= i0; j:= j0; t:= t0; m:= M[t0+1];
                while i*t <> -j*m do            
                    if (i*t - j*m) mod 2 <> 0 then p:= 1-p fi;
                    i:= i/2; t:= t/2; m:= m/2; j:= j/2
                od;
                S:= S+p
            od;
            A[i0+1, j0+1]:= S
        od
    od;
    0
end proc
:
BuildA:= Compiler:-Compile(BuildA_M):

LArip:= proc(n::posint)
local 
    r:= 2^n, 
    M:= Array(1..r, combinat:-randperm(r) -~ 1, datatype= integer[4]),
    A:= Array((1..r)$2, 0, datatype= integer[4])
;
    BuildA(A, M, r);
    (A,M)
end proc 
:
CodeTools:-Usage(LArip(8));
memory used=263.30KiB, alloc change=0 bytes, 
cpu time=16.55s, real time=16.68s, gc time=0ns

​​​​​​

The Float(undefined) could be the result of dividing by very small numbers. If you increase the precision with Digits:= 25 the problem goes away. I just guessed that 25 might work, and it did. Some smaller value may work also.

Change the line where tab_err is created to this:

tab_err:= Vector(N-1); 

Note that your code had = not :=, so it wasn't creating anything.

Change all references of tab_err[nrow, 1] to tab_err[nrow]. (There are only two such references.)

Change the plotting command to

print(plot(<<seq(1..N-1)> | tab_err>, axis[2]= [mode= log])); 

In Maple, all plotting commands do nothing more than create a data structure that represents a plot; they do not display the plot (this even includes the command named display). Displaying the plot requires a print command. The reason that plots get displayed by top-level commands is (essentially) that you automatically get a "free" print at the end of every execution group.

Almost always, a plot of the absolute values of errors vs number of iterations should be done logarithmically.

Replace the entire for loop that prints tab_err with the single command

printf("Printing tab_err\n%6.4e\nDone printing tab_err\n\n", tab_err)

Then you will see for certain that tab_err is being printed.

Before executing the procedure, give the command

Digits:= 15; #or a larger integer

The default value of 10 digits is not enough to appreciate what's happening in your procedure.

I think that what you've done is a brilliant hack that gets around Histogram's unnecessary refusal to handle the legend option.

Much of your code is redundant. The entire thing can be replaced by

plots:-display(
    evalindets~(
        Statistics:-Histogram~(
             Statistics:-Sample~(Normal~([0,1], 1), 1000),
             color=~ [gold, blue]
        ),
        specfunc(COLOUR),
        1@@0, #multi-argument identity function
        LEGEND~(["N(0,1)", "N(1,1)"])
    ),
    transparency= 0.3
);

Unlike my previous Answer to you, this one does produce exactly the same results as yours.

As shown by Kitonum, a simple use of solve proves that g must be a square except *possibly* if b*d*f = 0. So let's consider that exception. 

eq:= sqrt(a+b*sqrt(c+d*sqrt(e +f*sqrt(g)))) = h;
eq||(b,d,f):= eval~(eq, [b,d,f]=~0)[];
                               (1/2)    
                       eqb := a      = h

                                      (1/2)    
                        /       (1/2)\         
                 eqd := \a + b c     /      = h

                                            (1/2)    
                 /                    (1/2)\         
                 |      /       (1/2)\     |         
          eqf := \a + b \c + d e     /     /      = h

None of those equations contains g. If you find solutions for those (and there are many), you could use any values of g and the other missing variables. It'd be a total waste of computational effort to loop through the variables that do not appear. It's also a waste to loop through h: If the left side is an integer, then it is h, and it's the only possible h.

The best way to do this is to build the integral into the ODE system. In the following code t(s) represents your integral, so diff(t(s), s) is the integrand and t(0) = tinit is the initial condition.

#I just made up these parameter values:
g:= 9.81:  h:= 0.1:  rinit:= 0.1: tinit:= 5:

dsysA:= {
    diff(r(s), s) = sqrt((g*r(s) + h)^2 - 1 + 1/r(s)),
    diff(t(s), s) = (g*r(s) + h)/(1 - 1/r(s)),
    r(0) = rinit, t(0) = tinit
}:
dsnA:= dsolve(dsysA, numeric, abserr= 1e-6):

dsnA(0.1);
              
  [s = 0.1, t(s) = 4.89144107647636, r(s) = 0.429604740501696]

You need to be able to distiguish the free variable eps (which could even be an expression) from the variable r, which must be a variable because it'll temporarily be bound by int before being returned free (free and bound are opposites). There is no "clean" way to do this without passing r to f. By "clean", I mean without subverting the normal evaluation rules and without globally designating r as the temporarily bound variable. So, I'd do this:

f:= (g, r::name)-> int(g(r), r):

Then you could achieve your two desired calls as 

f(h, r);
f(h@(x-> x+eps), r);

This works whether g or h are defined procedures, just names, or some concoction in between those two. The shift function x-> x+eps could be predefined:

sh:= eps-> x-> x+eps:
f(h@sh(eps), r);

I recommend that you never use the syntax f(...):= ... to define a procedure (sometimes called a "function", but that word is ambiguous when talking about Maple). That syntax works sometimes, but not always.
 

If you want some numbers to be ignored, make them names rather than strings. You can do this by replacing "2222" with `2222` and likewise for the other string.

Thus, I get

By the way, your entire for loop can be replaced by map[inplace](evalf@log[10], B).

I agree that that is a bug. I've seen it many times before. It seems to be related to large fractional exponents. They cause solve to create RootOf expressions of polynomial degree the numerator of the fraction. So, in this case, that's a polynomial of degree 131537.

A workaround is to convert the exponent to a symbol:

Eq:= 1 - 1/(1+203808*exp(-342569/506*t)*(1/131537))^(131537/203808) =
44983/56599;
Eq2:= evalindets(
   Eq, anything^fraction, e-> op(1,e)^convert(op(2,e), name)
);
sol:= {solve}(Eq2, t);
evalindets(sol, name, parse);
evalf(%);

Acer's Answer overrideoption is the best way to do this. I just want to show you a primitive alternative to your primitive method.

#Your way:
c  := op(select(has, op(1, b1), COLOUR));
b2 := eval(b1, c=COLOUR(RGB, 1, 0, 0)):

#Replace with:
b2:= evalindets(b1, specfunc(COLOUR), c-> COLOUR(RGB, 1, 0, 0));

I mention this because, IMO, evalindets is the most important command in all of computer algebra and everyone should learn it.

My code does not necessarily do exactly the same thing as yours. Mine replaces all colors.

The conversion of an expression sequence to a list and vice versa is a triviality both in terms of coding and efficiency. Kitonum has shown the coding syntax. Regarding the efficiency, both the time and memory needed for the conversions are very small constants; they do not depend on the length of the data.

It is often better to stick with lists because of the indexing problems that can occur when an expression sequence has 0 or 1 elements.

This should not have been a separate Question, but since Tom already put an Answer here, I won't delete it.

To me, the input diff(x-1 = 0, x) is not GI (garbage input) nor is 1=0 GO (garbage output). As I was trying to explain in the other thread, everything makes sense if you take into account the quantifications that are implicitly present.

Maple allows for the simultaneous differentiation of both sides of an equation, and this is often useful. However, you could argue that this shouldn't be allowed (because of this very problem). And thus you could argue that diff(a = b, x) should be considered GI and an error for any a and b. If this restriction were to be imposed, I wouldn't object. It's trivial to use overload to impose that restriction:

restart:
diff_orig:= eval(diff):
unprotect(diff):
diff:= overload([
    proc(e::`=`)
    option overload; 
        error "Differentiation of equations not allowed."
    end proc, 
    diff_orig
]):
protect(diff, diff_orig):

The output 1=0 could be viewed as an algebraic equation or as a boolean relation. If it's an equation, then it's an equation whose solution set is empty. That's not garbage. If it's a relation, then it evaluates to false; false is not garbage.

The input and output considered together constitute a proof-by-contradiction that there does not exist any differentiable function f such that f(x) = x - 1 and f(x) = 0 (for all x in whatever domain of differentiation you want to consider). That's knowledge (albeit trivial), and knowledge isn't garbage.

You could do it in the way that you suggest if you want, but all the standard differential operators are in the VectorCalculus package: Gradient, Jacobian, Hessian, and several others.

As a workaround, you can do this:

with(Physics):
Physics:-Setup(assumingusesAssume= false):
csgn(a - b[2]) assuming a > b[2];

The simplify and all the other assumptions are superfluous. You can put them if you want, but they make no difference.