An equivalent to map(expand@expand, N^%H) is 

(expand~@@2)(N^*)

If you want some variables to represent real numbers, then you must use assume or assuming. If you want them all to represent reals, then why not just use the ordinary transpose?

@Laurenso I just realized that my procedure doesn't check whether the given point is actually on the sphere. If you require that check, I can easily add it.

@Laurenso The command to see the procedure is showstat(geom3d::onobjl).

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