Carl Love

## 26897 Reputation

11 years, 310 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

## Complete code for it...

You wrote:

• The bases that exit for one selected period (your proc) are pairwise conected to another base in the set.  I can see this in the plots.

So  it's possible to run only the half of the loop with some extra calculations for pairfinding. I think this may reduce the time to 60% of running the fulll loop.

Yes, both of those observations are essentially correct, but actually the time reduction is usually much greater. The "orbits" of a given base with respect to a given period are not just pairs but have size totient(period). So, for example, for period = 11, the bases come in groups of size totient(11) = 10. Knowing any one base in that orbit, the other 9 are trivial to generate. Here is a program for it:

```restart:
Periods:= module()
export
#prime factorization in convenient list form (see ?ifactors):
Ifactors:= (n::posint)-> (thisproc(n):= ifactors(n)[2]),

#prefactored form of the totient, making it more efficient to compute the
#Carmichael number at the same time:
TotSeq:= (n::posint)->
local p;
(thisproc(n):= [seq]((p[1]-1)*p[1]^(p[2]-1), p= Ifactors(n))),

#Euler's totient function:
Totient:= (n::posint)-> (thisproc(n):= mul(TotSeq(n))),

#Carmichael's lambda function:
Carmichael:= (n::And(posint, Not(1)))->
local T:= TotSeq(n);
(thisproc(n):= ilcm(`if`(n::even and T[1]>2, [T[1]/2, T[2..][]], T)[])),

#list of prime divisors of n:
PrimeDivisors:= (n::posint)-> (thisproc(n):= op~(1, Ifactors(n))),

#Carmichael number reduced by its prime divisors:
PreCarmichael:= (n::And(posint, Not(1)))->
local L:= Carmichael(n);
(thisproc(n):= iquo~(L, PrimeDivisors(L))),

(* Find any element of order Carmichael(n). The optional table Ex contains
elements to NOT check. The optional 3rd argument St is the element
Generator:= proc(
n::And(posint, Not({1,2})), Ex::table:= table('sparse'), St::posint:= 2
)
local E:= PreCarmichael(n), a, e;
for a from St to n-1 do
if Ex[a]=0 and igcd(a,n)=1 and andseq(a&^e mod n <> 1, e= E) then
return a
fi
od;
FAIL
end proc,

#Find the elements of order ord. The optional 3rd argument NG is the
#maximum number of generators to use.
Bases:= proc(
n::And(posint, Not({1,2})),
ord::And(posint, Not(1)),
NG::posint:= infinity
)
local
Ex:= table('sparse'), L:= Carmichael(n), LP:= iquo(L, ord),
E:= select(igcd=1, [\$1..ord-1], ord), EG:= select(igcd=1, [\$1..L-1], L),
G, e, g:= 1, r, ng:= 0, h, TT:= Totient(Totient(n))
;
if irem(L, ord) <> 0 then return {} fi;
E:= E[2..] -~ E[..-2];
EG:= EG[2..] -~ EG[..-2];
{
to NG while ng <> TT and (g:= Generator(n, Ex, g+1)) <> FAIL do
G:= e-> (thisproc(e):= g&^e mod n);
Ex[(h:= g)]:= 1; ng++;
for e in EG do
if Ex[(h:= h*G(e) mod n)]=0 then Ex[h]:= 1; ng++ fi
od;
r:= g&^LP mod n;
G:= e-> (thisproc(e):= r&^e mod n);
((h:= r), seq((h:= h*G(e) mod n), e= E))
od
}
end proc
;
end module
:
den:= 11842585:  #same number that you used
NumberTheory:-Divisors(Periods:-Carmichael(den));
{1, 2, 3, 4, 6, 7, 11, 12, 13, 14, 21, 22, 26, 28, 33, 39, 42,
44, 52, 66, 77, 78, 84, 91, 132, 143, 154, 156, 182, 231, 273,
286, 308, 364, 429, 462, 546, 572, 858, 924, 1001, 1092, 1716,
2002, 3003, 4004, 6006, 12012}

bases11:= CodeTools:-Usage(Periods:-Bases(den, 11));
memory used=287.81MiB, alloc change=0 bytes,
cpu time=4.28s, real time=3.79s, gc time=703.12ms

bases11 := {207496, 299716, 338141, 484156, 514896, 561006,
568691, 737761, 853036, 1029791, 1098956, 1237286, 1367931,
1452466, 1544686, 1590796, 1598481, 1752181, 1767551, 1844401,
2105691, 2128746, 2159486, 2267076, 2282446, 2359296, 2482256,
2620586, 2781971, 2797341, 3135481, 3143166, 3189276, 3312236,
3389086, 3573526, 3658061, 3673431, 4026941, 4188326, 4418876,
4457301, 4541836, 4603316, 4634056, 4733961, 4933771, 5202746,
5248856, 5448666, 5487091, 5663846, 5717641, 5733011, 5871341,
5963561, 6086521, 6148001, 6224851, 6247906, 6401606, 6601416,
6662896, 6747431, 6793541, 6901131, 7031776, 7208531, 7254641,
7277696, 7308436, 7431396, 7546671, 7631206, 7723426, 7777221,
8207581, 8292116, 8307486, 8338226, 8445816, 8660996, 8722476,
8799326, 8807011, 8822381, 8960711, 8976081, 9052931, 9091356,
9175891, 9268111, 9314221, 9321906, 9337276, 9368016, 9490976,
9606251, 9690786, 9752266, 9783006, 9852171, 9882911, 10082721,
10205681, 10397806, 10505396, 10636041, 10782056, 10812796,
10858906, 10912701, 11020291, 11035661, 11296951, 11373801,
11550556, 11627406, 11665831, 11811846}

```

## Last used for-loop index...

The last-used index of a for loop is still accessible after the loop ends. At that point, you can, for example, subtract 1 from it and use it. Like this:

for i to nops(testseq) while testseq[i] < 15 do od:
const:= `if`(i=1, FAIL, testseq[i-1]);

If all elements satisfy the while condition, then i will equal nops(testeq) + 1 at the end, and the last command will still work.

## try ... catch "time expired": ...; next ...

Yes, it can be done as in this example:

```for j to 99 do
p:= ithprime(j);
try
timelimit(1, ifactor(2^p-1))
catch "time expired":
printf("Time expired at j=%d; proceeding to next j.\n", j);
next  #use next j
end try
od:```

However, note that there are some commands that may use far more than their allowed amount of time before the timelimit command "catches" and stops them.

## Sequences also accepted...

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.

## You need to grant it permission...

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.

## Save from the numeric algorithm...

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:

1. Simply count those points.
2. 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

## In this case, there's only one...

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];

## Trivial to de-nest it...

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.

## Folding...

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 folding 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
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)(
foldl(
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;

## Ctrl-Delete...

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

## select...

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 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. There's no benefit to this if the period is square-free, as is 6

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

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

## try-catch statement...

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 - 61

f:= unapply(p, x):
infolevel[NM2]:= 1:
r:= NM2(f, 1.3);
1.121934407
.9609516325
.8274608051
.7423627912
.7136509949
.7110238945
.7110039767
.7110039755
.7110039755
r := 0.7110039755
f(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 I
f(r);
-7
0. + 1.42 10   I

```

## Cylindrical coordinates...

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]
);```