It's not mentioned on its help page, but sequences of plots structures are also accepted by **plots:-display**. These may be followed by additional options.

Many Maple commands accept additional forms of input that they aren't documented to accept. I'm not quite sure if that's a good thing or a bad thing.

The code for the module named **action** is in a part of the worksheet called "Startup Code", as mentioned by Joe Riel. When you load the worksheet, you need to grant permission to Maple to execute that code. You'll get a popup dialog box:

- "Autoexecute Warning: This worksheet contains content that will execute automatically. Do you wish to proceed?"

You need to click "Yes".

There are always a few language additions in each new version of Maple. In the vast majority of cases, older code continues to work. (That is a design goal called "backward compatibility".) Once I grant permission, it seems to work. Let me know if you have any more problems.

Although you haven't said this explicitly, the way that you present the Question suggests that you expect to count those black points via some sort of image processing of the final plot. However, it's much easier to just save those points when they are found by your numeric algorithm. They are, of course, those points for which the algorithm returns **0**.

There are two ways to do this. The one you use will be determined by whatever after-processing you want to do. They are both very easy to implement, but the 1st is a little bit easier than the 2nd:

- Simply count those points.
- Save those points.

Both cases involve adding exactly one line of code to **Hafiz1Count**, immediately before the **return 0**, and exactly one line outside the procedure before the **plot3d** is called. To simply count them, the line in the procedure is

**:-BlackCt++; #or :-BlackCt:= :-BlackCt+1**

and the outside line is

**BlackCt:= 0;**

To save the points, the line inside the procedure is

**:-BlackZ,= x0+y0*I; # Yes, the operator is "comma equal"!**

the outside line is

**BlackZ:= Array(1..0, datatype= complex[8]);**

and the number that have been saved can be retrieved by

**numelems(BlackZ);**

Both of my suggestions for the inside line require that the procedure be entered in 1D input.

In both cases, the outside line, which initializes the counter or array, must be re-executed before redoing any plots that use **Hafiz1Count**.

You wrote:

- I would expect the way to do this would be
`ans_all[1]`

.

That would be true except for cases where **allvalues **only returns 1 solution, which is what is happening in this case. If you want an expression that'll work regardless of whether there's 1 or more than 1 solution, use

**[ans_all][1];**

Your inner procedure **computus** does not use any lexical variables. Thus, it is utterly trivial to de-nest it! With a few trivial* syntactic* (not *algorithmic*) modifications, it can be compiled.

The basic technique for creating dynamic nested loops is to apply the **foldl** (fold-left) command to **seq **or one of its companion looping commands such as **add**. This is what @sursumCorda did. As he showed, making the loops dynamic by **fold**ing doesn't make them any faster, nor should you expect it to. But I did find 2 major ways to make the code faster, speeding up yours and @sursumCorda 's by a factor of at least 40 (i.e., 40 times faster).

The first improvement was to eliminate your innermost loop. That loop is simply counting (by repeatedly adding **1**) a sequence of consecutive integers. If you know the integer bounds, the count can be done by simple subtraction of those bounds, plus 1. That's the significance of this expression, which is the innermost addend in my **add** loops:

**max(0, min(_v||n1 + 1, V[-1]) - Ceil(V[-1]/2))**

The second improvement is due to a Maple idiosyncracy: The command **trunc **is many times faster than **ceil** (or **floor** or the other commands of that ilk). By making a slight adjustment to **trunc**, it can replace **ceil**:

**ceil(x) = trunc(x) + `if`(trunc(x)=x or x < 0, 0, 1)**

With these changes, my final code is able to do the **n=6** case in under 9 seconds.

nequaln:= (n::And(posint, Not(1)))-> local i, %add, n1:= n-1, V:= (n-2)*24 -~ (0, seq['scan'= `+`](_v||i, i= 2..n1)),#seq of partial sums#Maple's library 'ceil' has a highly symbolic aspect that slows it significantly. (You #can confirm this by showstat(ceil).) But 'trunc' is kernel builtin and doesn't have #this problem. Thus, I wrote this version of 'ceil' that only uses 'trunc'.Ceil:= x-> local T:= trunc(x); `if`(T::integer, (thisproc(x):= T + `if`(T=x or x<0, 0, 1)), 'procname'(x)) ; (eval@subs)( _v||1= V[1] - n1, %add= add, foldl( %add, max(0, min(_v||n1 + 1, V[-1]) - Ceil(V[-1]/2)), seq(_v||i= Ceil(V[i-1]/(n-i+2)).._v||(i-1), i= n1..2, -1) ) ) : [seq](CodeTools:-Usage(nequaln(k)), k= 2..6);memory used=14.30KiB, alloc change=0 bytes, #Usage for k=2 cpu time=0ns, real time=2.00ms, gc time=0ns memory used=285.20KiB, alloc change=0 bytes, #Usage for k=3 cpu time=15.00ms, real time=14.00ms, gc time=0ns memory used=475.34KiB, alloc change=0 bytes, #Usage for k=4 cpu time=0ns, real time=9.00ms, gc time=0ns memory used=13.45MiB, alloc change=0 bytes, #Usage for k=5 cpu time=79.00ms, real time=80.00ms, gc time=0ns memory used=1.58GiB, alloc change=0 bytes, #Usage for k=6 cpu time=8.58s, real time=8.35s, gc time=437.50ms [0, 48, 816, 10642, 117788]

You're trying to use a variable named **delta[h]_range**. That's not allowed as an *unquoted* variable name. But, if you wish, you can use back-quotes (aka left-quotes) to make anything at all a variable name:

**`delta[h]_range`:= 0.55 .. 0.75;**

Click anywhere in the output field--no need to highlight it--and press Ctrl-Delete.

It can be done like this:

period:= 6: x:= 78: d:= period /~ NumberTheory:-PrimeFactors(period): select( b-> igcd(x,b)=1 and b&^period mod x = 1 and andseq(b&^i mod x <> 1, i= d), [$2..x-1] );[17, 23, 29, 35, 43, 49]

~~The asymptotic time complexity with respect to ~~There's no benefit to this if the period is square-free, as is **period** can be further reduced by using some kind of "factor tree" to compute the powers **b**** &^ i mod x**, where **i** runs over the proper divisors of **period**. **6**.

**Edit:** The idea mentioned in the previous paragraph is now implemented by the new 2nd line above:

**d:= period /~ NumberTheory:-PrimeFactors(period);**

You can put the evaluation of the Newton iteration in a **try**-**catch** statement to trap both singularities of **f** and points where its derivative is 0.

NM2:= proc( f, xold::complexcons, {precision::positive:= 10.^(1-Digits), max_iters::posint:= 9} ) local x0:= xold, x1:= x0+2*precision, Df:= D(f), k; for k to max_iters do try (x0,x1):= (x1, evalf(x1 - f(x1)/Df(x1))) catch: error "Singularity at %1", x1 end try; userinfo(1, procname, 'NoName', x1) until abs(x0-x1) < precision; if k > max_iters then error "Maximum iterations exceeded" fi; x1 end proc : p:= randpoly(x, degree= 9, dense);9 8 7 6 5 4 3 p := 72 x + 37 x - 23 x + 87 x + 44 x + 29 x + 98 x 2 - 23 x + 10 x - 61f:= unapply(p, x): infolevel[NM2]:= 1: r:= NM2(f, 1.3);1.121934407 .9609516325 .8274608051 .7423627912 .7136509949 .7110238945 .7110039767 .7110039755 .7110039755 r := 0.7110039755f(r);0.r:= NM2(f, 1.+I);.9145078329+.9127717630*I .8599212591+.8460365442*I .8403634004+.8059364156*I .8411389647+.7949972141*I .8416536274+.7947359908*I .8416536540+.7947375421*I .8416536539+.7947375422*I r := 0.8416536539 + 0.7947375422 If(r);-7 0. + 1.42 10 I

I think that the following is what you have in mind, although it's not quite correct to say that it's the surface *common to *the two cylinders. Rather, it's the portion of each cylinder that's inside the other.

plots:-display( plot3d(# y^2 + z^2 = 1:[[r, t, sqrt(1-r^2*sin(t)^2)], [r, t, -sqrt(1-r^2*sin(t)^2)]], r= 0..1, t= -Pi..Pi, coords= cylindrical ), plot3d(# x^2 + y^2 = 1:[1, t, z], t= -Pi..Pi, z= -sqrt(1-sin(t)^2)..sqrt(1-sin(t)^2), coords= cylindrical ), scaling= constrained, labels= [x,y,z] );