@2cUniverse I am almost certain that all of the bases that my program says are of a given order (or period) are indeed of that order. The error is that it doesn't necessarily produce all bases of that order. 

Yes, there was a fundamental theoretical error in my algorithm. I don't know if there's an efficient way to fix it, but I remain hopeful. Hopefully someone with a deeper knowledge of Number Theory, specifically of the non-cyclic multiplicative groups modulo n, will read this and advise. 

Here's a theoretical description of my error. This requires no knowledge of Maple, or even programming, to understand; just knowledge of Number Theory.
 

Definitions and Notations:

All lowercase variables represent nonnegative integers. All lowercase variables other than n represent nonnegative integers less than n.

Multiplicative group modulo n: Let Z[n]* denote the multiplicative group modulo n, in other words, the set of positive integers less than n that are coprime with n. This is also called the group of units modulo n.

Order: The order of g in Z[n]* is the smallest positive integer e such that g^e = 1 (mod n).

Totient (Euler's totient): The totient of n is the number of positive integers less than n and coprime with n. Thus, totient(n) = |Z[n]*|, and the order of any element is a divisor of the totient. When used as a prefix functional operator symbol, it is usually denoted by the lowercase Greek letter phi. The totient can be very easily calculated from the prime factorization of n. 

Carmichael's lambda function: The Carmichael function of n is the largest order of any element of Z[n]*. It is also called the reduced totient. It is a divisor of the totient, and the order of any element is a divisor of it. Its prefix functional operator symbol is usually lowercase Greek lambda. It can also be easily calculated from the prime factorization of n.

Primitive root: An element of Z[n]* of order totient(n) is called a primitive root of n. Primitive roots only exist for some n. When they do exist, the number of them is totient(totient(n)).
 

Cyclic vs Non-cyclic groups of units: 

The assumption that n is NOT of the form 1, 2, 4, p^k, 2*p^k (where p is an odd prime) is equivalent to all of the following:

• Z[n]* is not a cyclic group.
• Carmichael(n) < totient(n).
• There does not exist an element in the group whose order is totient(n).
• The group does not have a primitive root.

The specific false assumptions that I made:

Let L = Carmichael(n). Let G = { g in Z[n]* | order(g) = L }.

Falsehood 1: I had assumed that for every x in Z[n]* there existed g in G and integer e such that g^e = x (mod n). In other words, I assumed that the powers of the elements of G "generated" the whole group (hence the symbol G).

However, that does seem to be true for most x. When it's true, we have this proposition:

Proposition: order(x) = L / gcd(e, L).

The consequence of this false belief was that I was trying to use that process to find all x of that order.

Falsehood 2: I had assumed that |G| = totient(totient(n)). Although this is certainly true if Z[n]* has a primitive root, in general |G| <> totient(totient(n)). For example, for n = 8, G = {3,5,7}, so |G| = 3, but totient(totient(n)) = 2.

The consequence of this false belief was that my algorithm may stop trying to find new elements of G too soon.
 

Reduction of the size of G for more efficient computation:

This process, which is already implemented in my algorithm, can be used to replace G with a much smaller subset GR that has identical computational power: 

Define an equivalence relation R on G via g R h iff there exists an e such that g^e = h (mod n). Any such e is necessarily coprime with L, so the size of the equivalence classes is totient(L). Let GR be any system of representatives of the equivalence classes. So, the size of GR is |G|/totient(L). (If Z[n]* has primitive roots, then GR has a single element, which could be any of the primitive roots.)

The way that this is implemented in my algorithm is more efficient than the naive interpretation of the definition of R may imply: Immediately after any g in G is found, all members of its equivalence class are cached and removed from further consideration.

Conjecture: For all postive integers n,

totient(totient(n)) / totient(Carmichael(n)) = totient(n) / Carmichael(n).

This seems to be true even if |G| <> totient(totient(n)).

Questions; Request for assistance:

1. Pseudo-primitive roots:
So, how can I enlarge GR so that it can function akin to a "set of primitive roots"? In other words, how can it be made into some sort of "basis" or minimal set of generators which can generate all of Z[n]*? Would taking products of distinct elements of GR be sufficient enlargement? What if they were only pairwise products of distinct elements? Should I use pseudo-primitive roots of n (see ?NumberTheory, PseudoPrimitiveRoot)? A potential problem with that with respect to computational efficiency is that I'd like to avoid computing pseudo-primitive roots that are equivalent under relation R (defined in the previous section) to ones already computed.

2. Counting elements of a given order:
If I could count all the elements of a given order before they are generated, the speed of the algorithm could be greatly improved. I realize that there is no simple closed-form formula for this, but perhaps there is a recurrence relation (perhaps akin to partitions of an integer). Computationally, a recurrence would be just as good as a closed-form.

@ianmccr Your original 1-D version was

CP2:= (X,Y)-> (local x,y; {seq(seq([x,y], y= Y), x= X)});
#             ^                                        ^ 

It was the parentheses that I've highlighted in red that were causing the problem. Parentheses enclosing the whole operator definition are not needed in this case. However, they are often needed when the operator definition is embedded within another expression. The correct parentheses placement is

CP2:= ((X,Y)-> local x,y; {seq(seq([x,y], y= Y), x= X)});
#     ^                                                ^

@2cUniverse There is a problem with the code above: It doesn't find some bases. I think that the problem is with my algorithm itself, not my Maple coding of it. In particular, it's missing some bases of order 2 and some of order Carmichael(n)/2. I'm still trying to solve it.

@2cUniverse Sorry, I mistakenly thought that you said that you were using Maple 2023. What version are you using? I modified my uses of $ so that they'll work for earlier Maple. Here is the worksheet link:

Carmichael.mw

@MalakMMK Thanks. You did that so quickly. You've made quite-sophisticated changes to the code. You know a lot more Maple than I originally thought. 

Any Maple computation that can be done in 2D Input can be done in 1D. (But the converse of that is very far from true.) Indeed, any 2D Input is internally converted to 1D as the input is being parsed.

Since you've made a bona fide effort to work with me on this, I'll show you how to do for loops and whatnot in 1D. It might take me a day to get to it.

@MalakMMK Okay, I read the paper. They're using Mathematica, not Maple. I don't think that the difference in the black-point count is significant. The only difference that might be significant is the comparison with the counts for the other algorithms.

It's possible to rewrite the code so that every arithmetic step is done in 64-bit hardware floating point "double precision", and then perhaps it would function exactly like the authors' version. Since (as far as I can see), the authors didn't provide actual code, that's hard to say for sure.

I'm not willing to continue working with code in 2D Input. If you want further advice on this from me, retype it in 1D. The first and last sections of your worksheet are already in 1D input anyway.

> 1D input is the stuff in red boldface monospaced characters, like this.

@MalakMMK I don't see where in the code N occurs. Is it the loop upper bound, which is 20 in your original worksheet? The difference is likely due to a different model of floating-point computation with different rounding errors. Were the authors using Maple?

@MalakMMK Your error is due to using 2D Input (the default input mode). It doesn't recognize the ,= or ++ operators (as well as a host of other modern niceties).

Perhaps going against my better judgement, I'll show you how to do this in 2D Input. (But, fair warning, I'm not about to go down the nearly bottomless rabbit hole of correcting errors caused by 2D Input.)

Make the line in the procedure
    :-BlackCt:= :-BlackCt + 1;  :-BlackZ(:-BlackCt):= x0 + I*y0;
and make the line outside the procedure
    BlackCt:= 0:  BlackZ:= Array(1..0, datatype= complex[8]):
The indexing of BlackZ inside the procedure must be done with parentheses ( ), as shown, not with the usual square brackets [ ].

If you run this on the code exactly as you had it, there are indeed 0 black points. But if you change the function F to z^4 - 1, then there will be 84.

Please use the green uparrow on the toolbar of the MaplePrimes editor to upload a link to the file that you're trying to upload to Maple. Or, if the file is publically accessible on the Internet, you may just give the URL.

@2cUniverse I ask again Why do you use square-free numbers? The techniques that we're discussing here will work for any positive integer (other than 1).

Use the green uparrow on the toolbar of the MaplePrimes editor to upload your worksheet. Even if it can't display your worksheet inline (which often happens), it will still put in a link such that your worksheet can be downloaded.

For the example number that you show (den = 11842585), we have Totient(den)/CarmichaelLambda(den) = 704, which is fairly high. I can probably speed it up some, but when that ratio is large, it becomes more difficult. It can also be speeded up significantly if you'd be content with some bases of the specified period rather than all bases.

Please define "divergence point". Is it any point every neighborhood of which intersects all basins?

@jonsimoniowa I confirm your report that there is a bug with respect to .graphml files in Maple 2023.1 (the current version). Futhermore, the bug is specifically in the Import. I inspected the result of your Export with a plain-text editor, and the file contained the correct weights. (The erroneous weight matrix is <0, 2, 2; 3, 0, 3; 4, 4, 0>.) Using ExportGraph and/or ImportGraph produce identical erroneous results.

The erroneous values aren't necessarily lower than the originals. It seems like the first n nonzero values encountered (reading left-to-right, then top-to-bottom) are duplicated as constant values for the nonzero entries in each of the n rows.

@2cUniverse You wrote:

• I have checked numbers up to 10^8 with your code. With MultiplicativeOrder its possible up t 10^6.

  I use this for visualizing small periods (< 40)  of a period-base combination (see attached file) .

  My code produces a lineplot which has a lot of polygones. Then I fill the Polygones random with colors.

I don't see any attached file in your Reply. If you attach the file, I can advise further on speeding up the code. For example, if you want to compute the multiplicative order of all members of the multiplicative group modulo n, there are ways to do that that are much much faster than the sum of the times to compute them individually.

@2cUniverse I think that this is what you're looking for: Carmichael's lambda function returns the largest divisor of the totient that can be a period, and all other periods will be divisors of this number. So

NumberTheory:-CarmichaelLambda(78);
                               12

NumberTheory:-CarmichaelLambda(123);
                               40

Like the totient, this function can be easily computed from the prime factorization of its argument. Details are given in this Wikipedia article https://en.wikipedia.org/wiki/Carmichael_function

CarmichaelLambda(n) is the largest possible value of MultiplicativeOrder(a,n), and any value of MultiplicativeOrder(a,n) is a divisor of CarmichaelLambda(n).