@2cUniverse I think that you've found some relevant info. It'll take me a few days to digest it. Thanks.

Regarding your Reply before the one above: As I said before, that idea is already implemented in the presented code. Once any element x of order k is found, I immediately use it to generate another totient(k) - 1 elements of order k. What I'm missing is a fast way to generate the initial elements x needed to get all elements of order k. Starting with elements of order Carmichael(n) isn't enough.

So, do you want all pairs of perpendicular diameters with integer endpoints? 

It looks like you're trying to solve a system of 41 equations with only 8 unknowns.

@dharr I concur: vb has been changed from 0.5 to 5. Since vb is the numerator of the arguments of arccos, that change would have a major impact on the chance of there being a real solution (since arccos(x) is not real for real x such that abs(x) > 1).

Using the original vb = 0.5, I get this solution:

sol := {A = -1.261525434, B = 1.459322610, C = 0.03167646665,
    omega = 2.083929566, x1 = -1.965825318, x2 = -0.8500707662}

@dharr If NLPSolve detected complex numbers, then its error message would directly mention that. Since the error says "non-numeric", the input very likely still contains either symbolic non-decision variables or items with inappropriate indices (such as 2[1]).

I'm not saying that the input doesn't generate complex values; I don't know whether it does. I'm just saying that that can't be the main issue.

As you said, the OP needs to upload the worksheet.

@dharr The proximal cause of the error is that SparseDirectMKL is declared local to LinearAlgebra (for some unknown and likely unintentional reason), but the global SparseDirectMKL is being passed by the user's command. If this alone is the problem, it's easy to workaround; but there's another problem also:

kernelopts(opaquemodules= false):
:-SparseDirectMKL:= LinearAlgebra:-SparseDirectMKL:
kernelopts(opaquemodules= true):
sol:= LinearAlgebra:-LinearSolve(A, b, method= SparseDirectMKL);
Error, (in LinearAlgebra:-LinearSolve) data incompatible with method = LinearAlgebra:-SparseDirectMKL
 

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