@sigl1982 It's quite refreshing to see such neatly written code posted here.

I suspect that there is some intentional randomization somewhere in the mod code that your program uses, perhaps in `mod/Factors`. Do this test with randomize: Make the first statement of your procedure randomize(1). Run the procedure several times with the same input. (Q1) Do you get identical results? Then change the first statement to randomize(2). (Q2) Do you now get different results?

• If Q1 and Q2 are both "yes", that's very very strong evidence that intentional randomization is being used somewhere.
• If Q1 is "no", then there's no point in asking Q2, and you have strong evidence that intentional randomization is not being used. Back to the drawing board.
• If Q1 is "yes" and Q2 is "no", then try changing the 2 to other numbers and re-asking Q2 until the answer is "yes". If you can't make the answer "yes" by changing the argument to randomize, then the test is inconclusive.

I recommend that you turn off the ASSERTion checking while doing this simply to reduce the volume of infomation that you need to look through.

Randomized algorithms are commonly used in finite-field polynomial arithmetic. No doubt you can find several such algorithms in the book Modern Computer Algebra by von zur Gaithen and Gerhard. (Gerhard works at Maplesoft, by the way.)

As usual, it'd help if you post your code.

In general, the order of sets of immutable objects (such as algebraic expressions, numbers, sets, lists) is fixed, and it doesn't change between runs, even runs on different computers. However, if your sets contain mutable objects like Vectors, Matrices, or tables, then the ordering may be based on address---I'll need to research that---and of course addresses change between runs. Mutable objects also possibly include modules (which includes Objects and Records); I'll need to research that also.

However, even if your set order is based on something ephemeral like addresses, there's no way that the order is being changed by set operations such as union, minus, and indexing. The order within any given session (i.e., between restarts) is fixed.

In command-line Maple there's a command-line flag that causes Maple to use different set orders. The purpose of this is to test whether one's code depends on set order.

I think that you'll get a better efficiency gain by switching to irem on the inside:

mods(mul(irem(coeff(g[i], x, 0), p), i= S), p);

You just need to wait for the result. On my computer it took about 5 minutes.

Just a small change is needed: taylor(y(t__n+h), h, 4). This is equivalent to making h the main variable (the differentiation variable) and expanding around h=0.

The keys to doing this are using the Units package, the function-definition operator ->, and the commands solve and evalf. A benefit of using Units is that if the units come out correct, that's a fairly good check that you've set up the problem correctly.

The additional key to doing it neatly is to use eval to substitute numeric values so that that step can be put off until the end of each computation.

Download airplane.mw

I've made an extensive investigation of this twice over the years. I haven't reinvestigated with Maple 2016, so this may have changed. What I found was:

• The primary order of layering of the objects when a PLOT structure is displayed is determined by the object type (CURVES, POLYGONS, POINTS, etc.).
• The order cannot be changed by changing the order that the objects appear in the PLOT structure.
• The objects most in the background are those of two topological dimensions such as area fillls.
• Next are objects of one topological dimension such as CURVES and the boundaries of POLYGONS.
• Next are objects of zero topological dimension: POINTS.
• Next is plotted text.
• The object most in the foreground are the axes.
• Within each class of objects, the order of layering can be changed by changing the order in the PLOT structure, which can in turn be controlled by the order that objects appear in a display command.

The plot that you want can be obtained through numeric differentiation with command fdiff:

fy:= y0-> fdiff((r-> rhs(r[-1]))@pds:-value(f(x,y), x=0), [1], [y0]):
plot(fy, 0..10, labels= [y, ``]);

You need to define an "epsilon cloud" around each integer so that if x is sufficiently close to an integer, then f(x) will be 1.

f:= x-> piecewise(abs((x-round(x))/x) < Float(6, 1-Digits), 1, 0);

The value Float(6, 1-Digits) = 6*10^(1-Digits) is often used as the target accuracy for floating-point computation.

Try any of

numtheory:-cfrac([1,1,2]);
seq(numtheory:-cfrac([1$k, 2]), k= 1..3);
Value(NumberTheory:-ContinuedFraction([1,1,2]));
seq(Value(NumberTheory:-ContinuedFraction([1$k, 2])), k= 1..3);

When using Maple's 2D Input, you can't have a space between a function name and the following left parenthesis. These extra spaces get translated as multiplication operators, which is why you see * after sphereplot and fieldplot above.

The trouble with your method is that you repeatedly add elements to existing lists or sets. Unfortunately, that operation is very inefficient in Maple. It is so natural to approach this problem recursively; however, I can't think of any way to efficiently do it recursively in Maple, even after thinking about it for hours.

Here is a way to do it using seq. This method is faster than Iterator. (In fairness, I must say that I learned this technique from a post by Joe Riel where he applied it to Cartesian products.) I also show a better way to do your indexing procedure TkSubSetList. And I also compare with VV's recursive method.

Download AllCombinations.mw

Indices, which are formed with square brackets, and different from subscripts, which are formed with double underscore. The two forms display pretty much the same in the GUI. However, a subscript is an integral part of the name to which it is attached, and hence it can't be substituted for an integer or anything else. A subscript is intended for display purposes, not programmatic purposes.

Additionally, you'll achieve better results with mul rather than product (lowercase p). Here's your worksheet, modified:

Download WhatHappensHere.mw

Here's a much-simplified version of your code:

funcs := [x, 3*x, x^2, 5*x];
for i to nops(funcs) do listM[i] := diff(funcs, [x $ (i-1)]) end do;
M := convert(convert(listM, list), matrix);
linalg:-det(M);

P.S. Here's a more-direct way:

M:= matrix([seq(diff(funcs, [x$(i-1)]), i= 1..nops(funcs))]);

You think that D[3](y)(L) refers to the 3rd derivative with respect to the first variable. That's incorrect: Rather, it refers to the 1st derivative with respect to the 3rd variable (regardless of the fact that y is a univariate function). There are several ways to represent the 3rd derivative with respect to the first (and only) variable:

(D@@3)(y)(L)
D[1,1,1](y)(L)
D[1$3](y)(L)

Likewise for the 2nd or any other order of derivative.