Here's a naive implementation of a sampler that generates a uniform distribution.

Download SampleTriangle.mw

The specification of options for method= lsode follows an arcane protocol different than that of dsolve's other methods: Some of the options, including mxstep, must be passed in an array (not an Array). The size of the array must be (at least) 21 plus the number of ODEs after the system is converted to a first-order system. The position in the array of the options is documented at ?dsolve,numeric,lsode,advanced. Here's an example:

ode:= diff(y(x),x$4)=y(x):
ic:= y(0)=1, D(y)(0)=2, (D@@2)(y)(0)=3, (D@@3)(y)(0)=4:

#The size of the array must be 21 + #ODEs after conversion to a 1st-order system.
C:= convert(Array(1..25), array):    #The conversion to array is necessary.
C[2]:= 1e-5:   #Absolute error
C[10]:= 10000:  #mxstep: max # of steps allowed

Sol:= dsolve({ode,ic}, numeric, method= lsode, ctrl= C):

Your third question was

when I use lsode which way of numerical solution is applied (Euler,midpoint, rk3, rk4, rkf, heun, ... )?

Hmm, if you don't already know the answer to that, then why did you choose to use lsode---the most complex and difficult-to-use of dsolve's numerous methods? The short answer can be found at ?dsolve,numeric,lsode where it says "an implicit Adams method of evaluation [using] functional iteration (no Jacobian matrix is involved)."

The following procedure generalizes the Answers by Kitonum and One Man so that they'll work for any expression:

Display:= proc(e::algebraic)
     local r:= eval(e);
     print(nprintf("%a = %Q", e, subsindets(e, name, n-> nprintf("%Q", eval(n)))) = r);
     r;
end proc:

To use it, the desired expression must be entered with single quotes around every variable:

(a,b):= (3,4):
c:= Display('a'+'b'):

Now, the value of c has been assigned 7.

Of course, if you use this a lot, you may want to shorten the name Display to a single character.

The code is simple enough to analyze. The relevant code can be seen by

kernelopts(opaquemodules= false):
interface(verboseproc= 2):
eval(combinat:-permute1);

Note that I don't display the actual code here out of respect for copyright issues.

The code is similar to your procedure and only slightly more complicated. I believe that the significant difference is that the Maple code uses a remember table. That means it requires memory for a list of all the permutations, plus lists of all the permutations of all "tails" of the original input list.

Instead of break, put a readstat statement in your code, like this:

for k to 10 do
     if k=5 then readstat("Press ENTER or RETURN to continue:") end if
end do;

For more flexibility, you can start the debugger, like this:

Stop:= ()-> [][]:   #dummy procedure for stopat
stopat(Stop):   #Starts debugger when Stop is executed.

for k to 10 do
     if k=5 then Stop() end if
end do;

If you use sum, then L(i) is evaluated before i is given numeric values. This is called premature evaluation. There are many solutions. The simplest and most efficient is to use add instead of sum. With the special sequencing operators add, mul, and seq, the first argument is not evaluated until after the index variable is replaced by the appropriate values, and it is re-evaluated for each such value.

Premature evaluation is not a bug of sum. For the vast majority of Maple commands, including sum, arguments are fully evaluated before being passed. The evaluation is only condidered premature if it is done before the user wants it to be.

Whatever result you were trying to copy into your post didn't copy. But assuming that it was of the form Sol:= RealRange(a,b), then the b can be extracted as op(2,Sol). If b is of the form Open(c), then c can be extracted as op(op(2,Sol)) or op([2,1], Sol).

ListTools:-PartialSums(A);

It's not necessary that the message be less than N. If that were true, it would severly limit the size of messages which would be feasible to encrypt. Rather, you divide the message into chunks such that each chunk is less than N, and then you encrypt the chunks. It's convenient to use a fixed chunk size of an integer number of bytes. Using your generation technique, I will compute an appropriate chunk size:

(u,v):= 'nextprime(RandomTools:-Generate(integer(range= 10^100..10^(110-1))))' $ 2:
N:= u*v:
C:= iquo(ilog2(N), 8);

It follows that any chunk of 89 bytes or less can be represented as an integer less than N.

I can't test it without the code for RECTpred, RECTpredOR, RECTpredOL. But why aren't you using Threads:-Seq? That's much easier to use and manage than Threads:-Create.

The library commands abs and ceil are highly symbolic, and they're not thread safe (see ?index,threadsafe). They can't be used in code that uses Threads. If their arguments are strictly numeric, then you don't need the symbolic aspects, and you can write trivial thread-safe versions of them:

Abs:= x-> `if`(x < 0, -x, x);

Ceil:= proc(x)
local t:= trunc(x);
     `if`(x=t, x, `if`(x<0, t, t+1))
end proc:

Your code doesn't work because you are just reusing the same Matrix M each time. When you change the entries in a Matrix, it's still considered the same Matrix. This is different than the behavior of a list, which is considered immutable, meaning that if an entry changes, then it's considered an entirely new list. Try this:

restart:

CartProdSeq:= proc(L::seq(list), $)
local S,j;
    eval(subs(S= seq, foldl(S, [_i||(1..nargs)], seq(_i||j= L[j], j= nargs..1, -1))))
end proc:

PRIME:= 2:

for M in CartProdSeq([$0..PRIME-1] $ 16) do
     if Determinant(Matrix(M)) mod PRIME <> 0 then
          Mats[cat("", M[])]:= [][]
     end if
end do:

Mats:= [indices(Mats, nolist)]:

This code also handles your second request---it converts every Matrix to a compact string.

Edit: This has a bug. Corrected code is in my first Reply below.

Here are some bifurcation plots for this map. There is some quite interesting chaotic bifurcation for a small range of a.

Download bifurcation.mw

The "hairiness" of these plots is a defect of MaplePrimes. They don't appear like that when viewed in Maple.

Does this help?

a:= 1:  b:= 2:
Names:= <'a', 'b'>:
NewNames:= subsindets(Names, name, convert, `local`);

Names;

Now the names in NewNames will never evaluate.

The above won't work on indexed names such as a[0]. I'm working on code for that. (It's done: See my first Reply below.)

By the way, the purpose of the quotes `` isn't to prevent or delay evaluation; it's to construct names (specifically symbols) from strings that contain special characters or which would otherwise evaluate as code. So, `u0` is identical to u0 because it contains no special characters. The quotes are needed on `local` because local is a keyword which would evaluate as code.

Of course. Just add linestyle= spacedash to the verticalasymptoteoptions.

eq:= a*x+b*y=c*x+d*y:
`=`(selectremove(has, (lhs-rhs)(eq), x));