@sunit You can use the sign command in TransA, like this:

TransA:= proc(A)
local S:= sign(A), AA:= A/S;
    if AA = I*T[0]*(omega1-2*omega2) then S*(I*omega2*T[0]+I*si*T[2])
    elif AA = (3*I)*omega2*T[0] then S*(I*omega1*T[0]-I*si*T[2])
    elif AA = I*omegar*T[0] then S*(I*omega1*T[0]+I*si2*T[2])
    else A
    end if
end proc:

What kind of functions do you mean:

• "Standard" elementary functions like sin, cos, tan, ln, etc.?
• "Higher" mathematical functions like Gamma, erf, etc.?
• Functions that you define yourself?

@2cUniverse The "mod rule" that you're looking for to make your pairs is reciprocation. For example,

1/5 mod 78;
                               47
1/47 mod 78;
                               5

Reciprocals always have the same order. If the element has order 2, then it's its own reciprocal.

1/25 mod 78;
                               25

If an element x has order k, then its reciprocal is x^(k-1). For example, 17 has order 6, so

(1/17, 17^5) mod~ 78;
                             23, 23

If an element x has order k, and e is coprime with k, then x^e also has order k. For example, 7 has order 12. The coprimes of 12 are 1, 5, 7, 11, so

7^~(1,5,7,11) mod~ 78;
                         7, 37, 19, 67

@2cUniverse I think your intuition may be correct, and I was wrong about totient(totient(n)). It is well known that if there is a primitive root (a single element that generates the cyclic group), then the number of them is totient(totient(n)). But now I see that when there isn't a primitive root, there may be more than that number of elements of maximal order.

Consider Z[8]* = {1, 3, 5, 7}. The maximal order is 2, and there are 3 elements of that order ({3, 5, 7}). But totient(totient(8)) = 2, not 3.

So, I was wrong in my algorithm to stop looking for "generators" when the number found reached totient(totient(n)). But this alone cannot fully explain the error in my algorithm, for example when n = 78:

totient(78) = 24, Carmichael(78) = 12, totient(24) = 8; the subset of order 12 is G = {7, 11, 19, 37, 41, 59, 67, 71}, which has size 8. But the powers of these elements only generate 18 of the 24 elements of Z[78]*:
{seq}(seq(g^i mod 78, i= 1..12), g= G); nops(%);

@2cUniverse The totient(totient(n)) is the number of group elements of maximal order, regardless of whether that order is totient(n) or Carmichael(n).

However, I don't see what's intuitive about it.

@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.