GAP Sqrt() does not work on cyclotomics - bug in Unitaryrepresenation

[dimpase]

As pointed out in comments to https://mathoverflow.net/a/387421/11100

[ttiiddee]

You mentioned in mo that the norm equation can be solved by PARI-GP. Would you have a plan to publish the code or share the method? Is it possible to take the sqrt using PARI-GP? (I want to do the calculation, but I know nothing about PARI-GP and number theory.)

[ttiiddee]

@dimpase

1. I think I find a way to solve the equation $f^2=g$ over cyclotomics, see the PARI-GP code at the end. Basically, we first try to solve the equation in the cyclotomic field $K_n$ defined by $n$, where $g$ is in $K_n$. If it has a solution, then we are done. If there is no solution, then $T=x^2-g$ is an irreducible polynomial over $K_n$. So we can make a field extension $L/K_n$ defined by $T$. Then we note that $L/\mathbb{Q}$ is an abelian extension, iff $f^2=g$ has a solution in cyclotomics, according to the Kronecker-Weber theorem. So we check if $L/\mathbb{Q}$ is an abelian extension, if it isn't, then there is no solutions, if it is, we just need to find the cyclotomic field $K_m$ of which $L$ is a subfield. Finally, we solve the equation over $K_m$.

2. In the problem of unitarization of a group representation, we apply LDL decomposition to the symmetric bilinear form, and then we take the square root of $D$. I have some concerns that, it could happen that the square root of $D$ has no solution, but $D=U^{\dagger}U$ has solution in cyclotomics, where $U$ is a unitary. Do you have some arguments to rule out that kind of cases?

3. Could you teach me how to solve the norm equations? I have no idea about how to enlarge the cyclotomic field.

/* g is a polynomial with variable x as an element in cyclotomic field defined by primitive root of unity of x^n = 1. It returns the [rst, m, f], if rst = 1, f^2 = g such that f is in the cyclotomic field defined by m. if rst = 0, the m and f are meaningless. */

sqrtInUniversalCyclotomics(g, n) = {
   my(pol, K, rst, T, f, Lpolabs, rat, m, Kext, fext);

   /* Define the cyclotomic field K of n. */
   pol = polcyclo(n);
   /* Substitute x with y */
   pol = subst(pol, x, y);
   K = nfinit(pol);

   g = subst(g, x, y);
   T = x^2 - g;

   /* Solve the square root */
   rst = nfeltissquare(K, g, &f);
   f = nfbasistoalg(K, f);
   f = lift(f);
   f = subst(f, y ,x);
   if(rst,
       return ([1, n, f]); 
   ,
   /* If we fail to solve f^2 = g, it means we don't have a solution in K. We need a field extension L/K. We first check if the extension L/K is abelian, if it is, then the equation has a solution. */
   /* The minimal polynomial of the extension L/Q. */
   Lpolabs = rnfinit(K, T).polabs;
   rat = bnfinit(y); /* rational number */
   rnfisabelian(rat, Lpolabs),
       /* Next we need to find the cyclotomic field of m of which L is a subfield. */
       m = rnfconductor(rat, Lpolabs)[1][1][1 ,1];

       /* We define the cyclotomic field of m. */
       extpol = polcyclo(m);
       extpol = subst(extpol, x, y);
       Kext = nfinit(extpol);

       /* Solve the square root */
       g = subst(g, x, x^(m/n));
       rst = nfeltissquare(Kext, g, &fext);
       fext = nfbasistoalg(Kext, fext);
       fext = lift(fext);
       fext = subst(fext, y ,x);
       return ([1, m, fext]); 
   ,       
       return ([0, n, f]);
   );
};
/* Test */
n = 10;
g = 7;
[rst, m, f] = sqrtInUniversalCyclotomics(g, n);
[rst, m, f, Mod(f^2, polcyclo(m)) == Mod(g, polcyclo(m))]

[dimpase]

On August 5, 2025 6:46:15 AM PDT, ttiiddee ***@***.***> wrote: ttiiddee left a comment (gap-packages/RepnDecomp#22) @dimpase 1. I think I find a way to solve the equation $f^2=g$ over cyclotomics, see the PARI-GP code at the end. Basically, we first try to solve the equation in the cyclotomic field $K_n$ defined by $n$, where $g$ is in `K_n`. If it has a solution, then we are done. If there is no solution, then $T=x^2-g$ is an irreducible polynomial over $K_n$. So we can make a field extension $L/K_n$ defined by $T$. Then we note that $L/K_n$ is an abelian extension, iff $f^2=g$ has a solution in cyclotomics, according to the Kronecker-Weber theorem. So we check if $L/K_n$ is an abelian extension, if it isn't, then there is no solutions, if it is, we just need to find the cyclotomic field $K_m$ of which $L$ is a subfield. Finally, we solve the equation over $K_m$. 2. In the problem of unitarization of a group representation, we apply LDL decomposition to the symmetric bilinear form, and then we take the square root of $D$. I have some concerns that, it could happen that the square root of $D$ has no solution, but $D=U^{dagger}U$ has solution in cyclotomics, where $U$ is a unitary. Do you have some arguments to rule out that kind of cases? 3. Could you teach me how to solve the norm equations? I have no idea about how to enlarge the cyclotomic field. `/* g is a polynomial with variable x as an element in cyclotomic field defined by primitive root of unity of x^n = 1. It returns the [rst, m, f], if rst = 1, f^2 = g such that f is in the cyclotomic field defined by m. if rst = 0, the m and f are meaningless. */ sqrtInUniversalCyclotomics(g, n) = { my(pol, K, rst, T, f, Lpolabs, rat, m, Kext, fext); /* Define the cyclotomic field K of n. */ pol = polcyclo(n); /* Substitute x with y */ pol = subst(pol, x, y); K = nfinit(pol); g = subst(g, x, y); T = x^2 - g; /* Solve the square root */ rst = nfeltissquare(K, g, &f); f = nfbasistoalg(K, f); f = lift(f); f = subst(f, y ,x); if(rst, return ([1, n, f]); , /* If we fail to solve f^2 = g, it means we don't have a solution in K. We need a field extension L/K. We first check if the extension L/K is abelian, if it is, then the equation has a solution. */ /* The minimal polynomial of the extension L/Q. */ Lpolabs = rnfinit(K, T).polabs; rat = bnfinit(y); /* rational number */ rnfisabelian(rat, Lpolabs), /* Next we need to find the cyclotomic field of m of which L is a subfield. */ m = rnfconductor(rat, Lpolabs)[1][1][1 ,1]; /* We define the cyclotomic field of m. */ extpol = polcyclo(m); extpol = subst(extpol, x, y); Kext = nfinit(extpol); /* Solve the square root */ g = subst(g, x, x^(m/n)); rst = nfeltissquare(Kext, g, &fext); fext = nfbasistoalg(Kext, fext); fext = lift(fext); fext = subst(fext, y ,x); return ([1, m, fext]); , return ([0, n, f]); ); }; /* Test */ n = 10; g = 7; [rst, m, f] = sqrtInUniversalCyclotomics(g, n); [rst, m, f, Mod(f^2, polcyclo(m)) == Mod(g, polcyclo(m))]`

1) on vacation now, sorry for the slow reply. We have the code for solving norm equation in Pari-GP (as well as GAP code to use it to extend the package with the capacity to make representations real) it's in a private GitHub repo of my former student, I have to ask him to allow to make it public. 2) generally speaking, it's preferable to solve the norm equation rather than take square roots, as there are cases where the latter needs an extension of the field, while the former former does not. We work with LDL* decomposition in the unitary case, and so we can make good use of D_k=x_k x_k* decompositions. 3) see 1).

[ttiiddee]

@dimpase Thanks for your reply! Hope I don’t waste too much of your vacation time.
I think my point in 2) is that I’m worried about the norm equation $D_k=x_k x_k^*$ (or square root) don’t have a solution in cyclotomics. But the full matrix equation $D=UU^{\dagger}$ have a solution in cyclotomics. And you mentioned norm equation doesn’t need a field extension, does that mean the norm equation always has a solution? That’s a little bit suprised.

[dimpase]

the norm equation arising in my paper always has a solution, without an extension of the field.

in general, the norm equation (with RHS a totally real cyclotomic number) always has a solution in a cyclotomic extension, and a bound on the extension degree is known).

[ttiiddee]

@dimpase
So for the realization problem, since we are considering irreducible representation, and according to the theorem 1 you proved in the paper, it must have a solution.
And for the unitarization of representation, we still need a cyclotomic extension? Unless we also consider irreducible representation (this is not really a condition) and there is a similar theorem to ensure the solution exists.
(I hope there exists an algorithm to solve the norm equation of totally positive cyclotomic numbers that needs a cyclotomic extension. Because I'm actually working on representation of graded group and trying to find the Wedderburn block, and I also want to "unitarize" all the representation since I'm from physics. There are a bunch of different theorems we need to prove if there is no such algorithm to solve the norm equation...)

[dimpase]

an algorithm to solve the norm equation would just repeatedly extend the field until a solution is found.

[ttiiddee]

@dimpase Thanks for your help!
I finished the GP code and tested, it seems the algorithm inside rnfisnorm is very sensitive, I tried $\zeta_{n}$, $n=37$, if I use a simple $d$, then it solves very quickly, but if I use another element, then it meets stackoverflow. The same thing would happen if I change $n$.
I still don't know how to use pari in gap, it will still be very helpful, if the repo could be public, thanks again!

[dimpase]

it will be public, I'll let you know

[dimpase]

@ttiiddee please see https://github.com/OldAnchovyTopping/realrep for GAP/Pari code etc
(the final year project of the student who did the programming etc work is also there in report)
Let us know if you have questions about it.

[ttiiddee]

@dimpase Thanks a lot, it works fine.

By the way, I find a bug in the function galoissubcyclo of Pari. It is used to get the real cyclotomic polynomial, but when the real cyclotomic field is $\mathbb{Q}$, it will output x - a where a is not $\zeta_n + \zeta_n^{-1}$. The bug only appears when $n=3, 4, 6$, a should be -1, 0, 1 but in Pari, a are both 1. So if we replace the variable back to $\zeta_n + \zeta_n^{-1}$, it will go wrong.

This is actually a bug of Pari, because the manual says galoissubcyclo returns the Gaussian period.

Another thing is, maybe it would be better to use x^2-y*x+1 other than x^2+1, the former is actually $(x-\zeta_n)(x- \zeta_n^{-1})$. This may make the cyclotomic degree a little bit smaller in some cases.

[dimpase]

By the way, I find a bug in the function galoissubcyclo of Pari

can you provide concrete gp input to show this? I am not a Pari/GP user, and it's always a pain to figure out the correct format from their hard to grok docs.

[dimpase]

Another thing is, maybe it would be better to use x^2-y*x+1 other than x^2+1,

where exactly?

[ttiiddee]

@dimpase Here is the details:

1. galoissubcyclo(n, H) is a function to get the minimal polynomial of a subfield of $L_n=\mathbb{Q}(\zeta_n)$ that is fixed by the subgroup H of the galois group $Gal(L_n/\mathbb{Q})=\mathbb{Z}_n^{\times}$. If we let H = -1, it just means the subgroup of $\mathbb{Z}_n^{\times}$ generated by -1 which is the subgroup of complex conjugation. So galoissubcyclo(n, -1) is the real cyclotomic field $\mathbb{Q}(\zeta_n+\zeta_n^{-1})$.

Here is the bug, when $n=1,2,3,4,6$ (and it only appears in such cases.), we have

? galoissubcyclo(1,-1)
%66 = x - 1
? galoissubcyclo(2,-1)
%67 = x - 1
? galoissubcyclo(3,-1)
%61 = x - 1
? galoissubcyclo(4,-1)
%62 = x - 1
? galoissubcyclo(6,-1)
%63 = x - 1

And in the manual, it said

This function output is somewhat canonical, as it returns the minimal polynomial of a Gaussian period $Tr_{\mathbb{Q}(\zeta_n)/L}(\zeta_n).$

So we expect they should have (n=1 and 6 happens to be correct.)

? galoissubcyclo(1,-1)
%66 = x - 1
? galoissubcyclo(2,-1)
%67 = x + 1
? galoissubcyclo(3,-1)
%61 = x + 1
? galoissubcyclo(4,-1)
%62 = x
? galoissubcyclo(6,-1)
%63 = x - 1

So it is a bug. It seems the reason is that Pari make some special process to the case when $\mathbb{Q}(\zeta_n+\zeta_n^{-1})$ is $\mathbb{Q}$.

Note that $n=1,2$ are special, where $\mathbb{Q}(\zeta_n)/\mathbb{Q}(\zeta_n+\zeta_n^{-1})$ is $\mathbb{Q}/\mathbb{Q}$, the gaussian period is just $\zeta_n$, and other $n$ are degree two extensions, so the gaussian period is $\zeta_n+\zeta_n^{-1}$.

So the first two are still not the real cyclotomic polynomial, since the gaussian period is not $\zeta_n+\zeta_n^{-1}$. When $n=1$, the real cyclotomic polynomial is x - 2, when $n=2$, it is x + 2.

2. In your student's code

\\ Inputs: (to be added from GAP...)
\\ cyclotomic_degree = ??? ; desired_norm = ???

\\ Create real and complex cyclotomic fields:
\\ reals = bnfinit(galoissubcyclo(cyclotomic_degree,-1,0,y));
\\ ^^^^^ Supplied from GAP, because we may need it to define the norm.
complexes = rnfinit(reals, x^2 + 1);

\\ Make an element of a given norm:
norms = rnfisnorminit(reals.pol, complexes.pol, 1);
norm_root = rnfisnorm(norms, desired_norm);
\\ return:
cyclotomic_x = liftall(norm_root[1]);
print(cyclotomic_x);

rnfisnorminit is used to initialize the calculation, the input is K and L, where K is the minimal polynimial of $\mathbb{Q}(\zeta_n+\zeta_n^{-1})$, that is galoissubcyclo(n, -1), and L is the minimal polynomial of the simple extension $\mathbb{Q}(\zeta_n)/\mathbb{Q}(\zeta_n+\zeta_n^{-1})$.

The canonical way to do this is to use $(x-\zeta_n)(x-\zeta_n^{-1})=x^2-(\zeta_n+\zeta_n^{-1})x+1$ to define the extension. In pari, the base field K is represented by a polynomial over y, so our polynomial is x^2-y*x+1.

And in your student's code, he used x^2 + 1 to add the imaginary unit $i$. I think this is also okay yesterday, but now I find some problems. When $n$ is not a multiple of 4, the extension by adding $i$ may not be $\mathbb{Q}(\zeta_n)/\mathbb{Q}(\zeta_n+\zeta_n^{-1})$, but $\mathbb{Q}(\zeta_n+\zeta_n^{-1})(i)/\mathbb{Q}(\zeta_n+\zeta_n^{-1})$. So the norm equation is solved in a $\mathbb{Q}(\zeta_n+\zeta_n^{-1})(i)$, and it is not considered in the code. The following shows an example that goes wrong, it is because there is no solution in the field $\mathbb{Q}(\zeta_n+\zeta_n^{-1})(i)$.

gap> f:=E(7)+2*E(7)^2;
E(7)+2*E(7)^2
gap> fbar:=ComplexConjugate(f);
2*E(7)^5+E(7)^6
gap> g:=f*fbar;
-3*E(7)-5*E(7)^2-5*E(7)^3-5*E(7)^4-5*E(7)^5-3*E(7)^6
gap> sol:=PariNorms(7,g);
1
gap> sol*ComplexConjugate(sol)=g;
false

For the same reason above, we also need to treat $n=1,2$ as special cases, since $\mathbb{Q}(\zeta_n)/\mathbb{Q}(\zeta_n+\zeta_n^{-1})=\mathbb{Q}/\mathbb{Q}$, the number theoretical norm is also not the norm what we want.

[dimpase]

Here is the bug, when n = 1 , 2 , 3 , 4 , 6 (and it only appears in such cases.), we have

do you mean to say that it's a documentation bug ? (the function does return correct polynomials, but the documentation isn't correct on this).