https://hal.archives-ouvertes.fr/hal-02398953
Tip revision: c33910a29d53f4e137c225b21a8d59e43327cbf9 authored by Software Heritage on 08 December 2019, 12:26:32 UTC
hal: Deposit 351 in collection hal
hal: Deposit 351 in collection hal
Tip revision: c33910a
CMfield.cpp
#include "CMfield.h"
#include <cmath>
#include <iostream>
#include <string>
using namespace std;
using namespace giac;
map<string,gen> * sigma_values = 0;
map<string,gen> * T_values = 0;
CMfield::CMfield() {
verbose = false;
}
CMfield::CMfield(string s) {
verbose = false;
ten = gen("10",&ct);
t1 = gen("t1",&ct);
t2 = gen("t2",&ct);
vars1d.push_back(s__IDNT_e);
vars1d.push_back(t1);
vars2d.push_back(s__IDNT_e);
vars2d.push_back(t1);
vars2d.push_back(t2);
nd = 40;
ndmax = 200;
ndrecord = 30;
pf = 10;
zeronumber = eval(gen("0",&ct),&ct);
onenumber = eval(gen("1",&ct),&ct);
twonumber = eval(gen("2",&ct),&ct);
tag = s;
minpol = eval(gen("x-1",&ct),&ct);
genasrootof = onenumber;
if (tag.at(0)=='m') {
// Mostow group
unsigned k=1;
bool done = false;
string p_str = "";
while (!done && k<tag.size()) {
char te = tag.at(k);
done = te=='-';
if (!done) {
p_str = p_str + te;
}
k++;
}
if (!done) {
cerr << "Trouble reading tag" << endl;
exit(1);
}
done = false;
string tnum_str = "";
while (!done && k<tag.size()) {
char te = tag.at(k);
done = te=='%';
if (!done) {
tnum_str = tnum_str + te;
}
k++;
}
if (!done) {
cerr << "Trouble reading tag" << endl;
exit(1);
}
done = false;
string tden_str = "";
while (k<tag.size()) {
tden_str = tden_str + tag.at(k);
k++;
}
int pval = atoi(p_str.c_str());
int tnum = atoi(tnum_str.c_str());
int tden = atoi(tden_str.c_str());
if (verbose) {
cout << "p=" << pval << endl;
cout << "t=" << tnum << "/" << tden << endl;
}
gen arg_t3 = _simplify(((9*pval+2)*onenumber - 2*pval*tnum*onenumber/tden)/(4*pval*onenumber),&ct);
cout << "arg(t3)/2pi=" << arg_t3 << endl;
int arg_t3_num = _numer(arg_t3,&ct).val;
int arg_t3_den = _denom(arg_t3,&ct).val;
arg_t3_num = arg_t3_num % arg_t3_den;
cout << "arg(t3)/2pi num: " << arg_t3_num << endl;
cout << "arg(t3)/2pi den: " << arg_t3_den << endl;
gen t3;
if (arg_t3_num==0) {
t3 = onenumber;
}
else {
if (arg_t3_den==2) {
// num must be one
t3 = -onenumber;
}
else {
t3 = normal(_pow(makesequence(rootof(_cyclotomic(arg_t3_den,&ct),&ct),arg_t3_num),&ct),&ct);
}
}
// if (verbose) {
cout << "tau^3 as rootof: " << t3 << endl;
// }
if (!operator_equal(simplify(t3-onenumber,&ct),zeronumber,&ct) && !operator_equal(simplify(t3+onenumber,&ct),zeronumber,&ct)) {
extend(t3);
}
gen zeta = rootof(_cyclotomic(pval,&ct),&ct);
// if (verbose) {
cout << "zeta=" << zeta << endl;
// }
if (!operator_equal(simplify(zeta-onenumber,&ct),zeronumber,&ct) && !operator_equal(simplify(zeta+onenumber,&ct),zeronumber,&ct)) {
extend(zeta);
}
// At this stage should have a (large) number field that works for the whole group
// what follows is used only to compute the actual trace field.
// DM weights lcd of (1/2-1/p, 1/4+3/2p-t/2, 1/4+3/2p+t/2)
gen mu1 = onenumber/2 - onenumber/pval;
gen mu4 = onenumber/4 + 3*onenumber/(2*pval) - tnum*onenumber/(2*tden);
gen mu5 = onenumber/4 + 3*onenumber/(2*pval) + tnum*onenumber/(2*tden);
gen d1 = _denom(mu1,&ct);
gen d4 = _denom(mu4,&ct);
gen d5 = _denom(mu5,&ct);
vecteur dl;
dl.push_back(d1);
dl.push_back(d4);
dl.push_back(d5);
gen d = _lcm(dl,&ct);
// Should give optimal trace field
cout << "Computed taugen, d=" << d << endl;
taugen = rootof(_cyclotomic(d.val,&ct),&ct);
/* gen bla = tnum*onenumber;
bla = bla/(3*tden);
bla = bla + onenumber/pval;
bla = bla - onenumber/2;
int nu = _numer(bla,&ct).val;
int de = _denom(bla,&ct).val;
cout << "tau=exp(pi*i*" << nu << "/" << de << ")" << endl;*/
/* if (de!=1) {
// taugen = rootof(_cyclotomic(2*de,&ct),&ct);
if (nu>0) {
taugen = normal(pow(taugen,nu,&ct),&ct);
}
else {
if (nu<0) {
taugen = normal(onenumber/pow(taugen,-nu,&ct),&ct);
}
else {
taugen = onenumber;
}
}
}
else {
taugen = onenumber;
}
if (verbose) {
cout << "taugen=" << taugen << endl;
}*/
}
else {
if (tag.at(0)=='t') {
// thompson group
unsigned k=1;
bool done = false;
string p_str = "";
while (!done && k<tag.size()) {
char te = tag.at(k);
//cout << "te=" << te << endl;
done = (te=='H' || te=='S' || te=='E');
if (!done) {
p_str = p_str + te;
}
k++;
}
if (!done) {
cerr << "Trouble reading tag" << endl;
exit(1);
}
k--;
while (k<tag.size()) {
tautag = tautag + tag.at(k);
k++;
}
int pval = atoi(p_str.c_str());
if (verbose) {
cout << "p=" << pval << endl;
cout << "T=" << tautag << endl;
}
map<string,gen> m = * T_data();
map<string,gen>::iterator it = m.find(tautag), itend = m.end();
if (it==itend) {
cerr << "CM field not in data base" << endl;
exit(1);
}
if (pval>2) {
gen mu = rootof(_cyclotomic(pval,&ct),&ct);
extend(mu);
}
T = eval(it->second,1,&ct);
cout << "T values: " << T << endl;
cout << endl;
extend(T[0]);
extend(T[1]);
extend(T[2]);
}
else {
// should be sporadic...
unsigned k=0;
bool done = false;
string p_str = "";
while (!done && k<tag.size()) {
char te = tag.at(k);
done = te=='s';
if (!done) {
p_str = p_str + te;
}
k++;
}
if (!done) {
cerr << "Trouble reading tag" << endl;
exit(1);
}
tautag = "s";
while (k<tag.size()) {
tautag = tautag + tag.at(k);
k++;
}
int pval = atoi(p_str.c_str());
if (verbose) {
cout << "p=" << pval << endl;
cout << "tau=" << tautag << endl;
}
// Used to encode sporadic values
map<string,gen> m = * sigma_data();
map<string,gen>::iterator it = m.find(tautag), itend = m.end();
if (it==itend) {
cerr << "CM field not in data base" << endl;
exit(1);
}
taugen = eval(it->second,1,&ct);
cout << " this is also " << taugen << endl;
cout << endl;
extend(taugen);
if (pval>2) {
gen zeta = rootof(_cyclotomic(pval,&ct),&ct);
// gen zeta = eval(gen("exp(2*pi*i/"+print_INT_(pval)+")",&ct),1,&ct);
extend(zeta);
}
}
}
if (verbose) {
cout << "Done constructing number field, will now study real subfield" << endl;
}
findRealStructure();
vecteur e1o, e2o, e3o;
e1o.push_back(onenumber);
e1o.push_back(zeronumber);
e1o.push_back(zeronumber);
e2o.push_back(zeronumber);
e2o.push_back(onenumber);
e2o.push_back(zeronumber);
e3o.push_back(zeronumber);
e3o.push_back(zeronumber);
e3o.push_back(onenumber);
e1 = gen(e1o);
e2 = gen(e2o);
e3 = gen(e3o);
vecteur zv;
zv.push_back(zeronumber);
zv.push_back(zeronumber);
zv.push_back(zeronumber);
zerovector = gen(zv);
vecteur coid(9);
coid[0] = 1; coid[1] = 0; coid[2] = 0;
coid[3] = 0; coid[4] = 1; coid[5] = 0;
coid[6] = 0; coid[7] = 0; coid[8] = 1;
Id = createMatrix(coid);
vecteur ro1;
ro1.push_back(onenumber);
ro1.push_back(zeronumber);
vecteur ro2;
ro2.push_back(zeronumber);
ro2.push_back(onenumber);
vecteur ros;
ros.push_back(ro1);
ros.push_back(ro2);
Id2 = gen(ros);
// Not clear that what follows is needed...
gen gar = _evalf(genasrootof,&ct);
gen cxrts = _complexroot(makesequence(minpol,prec),&ct);
int countcheck = 0;
for (unsigned k=0; k<(*cxrts._VECTptr).size(); k++) {
if (is_greater(1e-6,_abs(cxrts[k][0]-gar,&ct),&ct)) {
// this does not always work, sometimes it is not an interval!!
if (operator_equal(cxrts[k][0],gen("i",&ct),&ct)) {
xvi = convert_interval(cxrts[k][0],nd,&ct);
}
else {
xvi = cxrts[k][0];
}
countcheck++;
}
}
if (countcheck==1) {
if (verbose) {
cout << "Found a unique interval containing evalf(genasrootof), will use " << xvi << endl;
}
}
else {
if (countcheck==0) {
cerr << "No interval contains evalf(genasrootof), this is probably due to roundoff errors in evalf " << endl;
exit(1);
}
else {
cerr << "Found several intervals containing evalf(genasrootof), I'd rather not use rootof." << endl;
exit(1);
}
}
}
CMfield::CMfield(gen pol, string s) {
verbose = false;
if (verbose) {
cout << "Creating CM field with tag " << s << endl;
}
ten = gen("10",&ct);
t1 = gen("t1",&ct);
t2 = gen("t2",&ct);
vars1d.push_back(s__IDNT_e);
vars1d.push_back(t1);
vars2d.push_back(s__IDNT_e);
vars2d.push_back(t1);
vars2d.push_back(t2);
nd = 40;
ndmax = 200;
ndrecord = 30;
pf = 10;
zeronumber = gen("0",&ct);
onenumber = gen("1",&ct);
twonumber = gen("2",&ct);
vecteur e1o, e2o, e3o;
e1o.push_back(onenumber);
e1o.push_back(zeronumber);
e1o.push_back(zeronumber);
e2o.push_back(zeronumber);
e2o.push_back(onenumber);
e2o.push_back(zeronumber);
e3o.push_back(zeronumber);
e3o.push_back(zeronumber);
e3o.push_back(onenumber);
e1 = gen(e1o);
e2 = gen(e2o);
e3 = gen(e3o);
vecteur zv;
zv.push_back(zeronumber);
zv.push_back(zeronumber);
zv.push_back(zeronumber);
zerovector = gen(zv);
vecteur coid(9);
coid[0] = 1; coid[1] = 0; coid[2] = 0;
coid[3] = 0; coid[4] = 1; coid[5] = 0;
coid[6] = 0; coid[7] = 0; coid[8] = 1;
Id = createMatrix(coid);
vecteur ro1;
ro1.push_back(onenumber);
ro1.push_back(zeronumber);
vecteur ro2;
ro2.push_back(zeronumber);
ro2.push_back(onenumber);
vecteur ros;
ros.push_back(ro1);
ros.push_back(ro2);
Id2 = gen(ros);
tag = s;
minpol = pol;
if (verbose) {
cout << "minpol: " << minpol << endl;
}
dc = _degree(minpol,&ct);
// if minpol is x^2+1, sto will not work!
gen tepo = simplify(minpol - gen("x^2+1",&ct),&ct);
if (!operator_equal(tepo,zeronumber,&ct)) {
genasrootof = rootof(minpol,&ct);
// if (verbose) {
cout << "Sto j with " << genasrootof << endl;
// }
gen jsto("j",&ct);
sto(jsto,genasrootof,&ct);
}
else {
genasrootof = rootof(minpol,&ct);
}
/* genasrootof = rootof(minpol,&ct);
gen jsto("j",&ct);
sto(jsto,genasrootof,&ct);
// if (verbose) {
cout << "Sto j with " << genasrootof << endl;
cout << "gen as rootof: " << genasrootof << endl;
cout << "Test interval (for gen as rootof): " << convert_interval(genasrootof,nd,&ct) << endl;
// }*/
gen gar = _evalf(genasrootof,&ct);
if (verbose) {
cout << "gar(2)=" << gar << endl;
cout << " evalf (for gen as rootof): " << gar << endl;
}
gen cxrts = _complexroot(makesequence(minpol,prec),&ct);
int countcheck = 0;
for (unsigned k=0; k<(*cxrts._VECTptr).size(); k++) {
if (is_greater(1e-6,_abs(cxrts[k][0]-gar,&ct),&ct)) {
// this does not always work, sometimes it is not an interval!!
if (operator_equal(cxrts[k][0],gen("i",&ct),&ct)) {
xvi = convert_interval(cxrts[k][0],nd,&ct);
}
else {
xvi = cxrts[k][0];
}
countcheck++;
}
}
if (countcheck==1) {
if (verbose) {
cout << "Found a unique interval containing evalf(genasrootof), will use " << xvi << endl;
}
}
else {
if (countcheck==0) {
cerr << "No interval contains evalf(genasrootof), this is probably due to roundoff errors in evalf " << endl;
exit(1);
}
else {
cerr << "Found several intervals containing evalf(genasrootof), I'd rather not use rootof." << endl;
exit(1);
}
}
if (!findInField(_subst(makesequence(minpol,x__IDNT_e,z__IDNT_e),&ct),conj(xvi,&ct),genconj)) {
cerr << "Trouble describing complex conjugation in field!" << endl;
exit(1);
}
if (verbose) {
cout << "Conjugate of generator: " << genconj << endl;
}
// cout << "gar: " << genasrootof << endl;
gen realgen1 = _simplify( (genasrootof + _subst(makesequence(genconj,x__IDNT_e,genasrootof),&ct)) ,&ct);
if (verbose) {
cout << "realgen1 = " << realgen1 << endl;
}
gen rmp1;
if (realgen1.type==0) {
rmp1 = s__IDNT_e-realgen1;
}
else {
rmp1 = expand(_poly2symb(makesequence(_pmin(eval(realgen1,1,&ct),&ct),s__IDNT_e),&ct),&ct);
}
if (!operator_equal(rmp1,s__IDNT_e,&ct) && operator_equal(_degree(makesequence(rmp1,s__IDNT_e),&ct),dc/2,&ct)) {
realgenasrootof = realgen1; // this is not right, will change it later
realminpol = rmp1;
if (verbose) {
cout << "Found a generator for real subfield: " << realgenasrootof << endl;
cout << "realminpol: " << rmp1 << endl;
}
}
else {
if (verbose) {
cout << "At first attempt, found degree " << _degree(makesequence(rmp1,s__IDNT_e),&ct) << endl;
}
gen realgen2 = _simplify((genasrootof*_subst(makesequence(genconj,x__IDNT_e,genasrootof),&ct)),&ct);
if (verbose) {
cout << "realgen2 = " << realgen2 << endl;
}
gen rmp2;
if (realgen2.type==0) {
rmp2 = s__IDNT_e-realgen1;
}
else {
rmp2 = expand(_poly2symb(makesequence(_pmin(eval(realgen2,1,&ct),&ct),s__IDNT_e),&ct),&ct);
}
if (!operator_equal(rmp2,s__IDNT_e,&ct) && operator_equal(_degree(makesequence(rmp2,s__IDNT_e),&ct),dc/2,&ct)) {
realgenasrootof = realgen2; // this is not right, will change it later
realminpol = rmp2;
if (verbose) {
cout << "Found a generator for real subfield: " << realgenasrootof << endl;
cout << "realminpol: " << rmp2 << endl;
}
}
else {
cout << "At second attempt, found degree " << _degree(makesequence(rmp2,s__IDNT_e),&ct) << endl;
cerr << "Trouble finding generator for real subfield" << endl;
exit(1);
}
}
gen rgar;
if (_degree(makesequence(realminpol,s__IDNT_e),&ct).val>1) {
rgar = rootof(realminpol,&ct);
gen jrsto("jr",&ct);
sto(jrsto,rgar,&ct);
}
else {
if (verbose) {
cout << "Setting gar to " << realgenasrootof << endl;
}
rgar = realgenasrootof;
}
if (verbose) {
cout << "real gen as rootof: " << realgenasrootof << endl;
}
// then need to compute realgen from realgenasrootof
gen yoo;
if (realgenasrootof.type==0) {
realgen = onenumber;
}
else {
if (contains(realgenasrootof,gen("j",&ct))) {
if (verbose) {
cout << "lidnt: " << lidnt(realgenasrootof) << endl;
}
int deg = _degree(makesequence(realgenasrootof,gen("j",&ct)),&ct).val;
vecteur coeffs;
for (int uu=deg; uu>=0; uu--) {
coeffs.push_back(_coeff(makesequence(realgenasrootof,gen("j",&ct),uu),&ct));
}
realgen = giac::expand(_poly2symb(makesequence(coeffs,x__IDNT_e),&ct),&ct);
}
else {
if (realgenasrootof.is_symb_of_sommet(at_prod)) {
yoo = _simplify(realgenasrootof[1][1]*realgenasrootof[2],&ct);
}
else {
if (realgenasrootof.is_symb_of_sommet(at_rootof)) {
yoo = _simplify(realgenasrootof[1],&ct);
}
else {
cout << "Problem reading factors!" << endl;
cout << "yoo=" << yoo << endl;
exit(1);
}
}
realgen = giac::expand(_poly2symb(makesequence(yoo,x__IDNT_e),&ct),&ct);
}
}
dr = _degree(makesequence(realminpol,s__IDNT_e),&ct);
for (int k=0; is_greater(dr,gen(k),&ct); k++) {
realminpolcoeffs.push_back(_coeff(makesequence(realminpol,s__IDNT_e,dr-k),&ct));
}
if (!2*dr==dc) {
cout << "This does not look like a CM field, degrees do not match!" << endl;
exit(1);
}
if (verbose) {
cout << "real min pol: " << realminpol << endl;
}
unsigned n = dr.val;
vars = vecteur(n);
for (unsigned j=0; j<n; ++j) {
vars[j] = gen("a"+print_INT_(n-j-1),&ct);
}
genericreal = _subst(makesequence(_poly2symb(makesequence(vars,s__IDNT_e),&ct),s__IDNT_e,realgen),&ct);
red(genericreal);
xv = _evalf(genasrootof,&ct);
if (verbose) {
cout << "xv:" << xv << endl;
}
sv = re(_evalf(realgenasrootof,&ct),&ct);
if (verbose) {
cout << "sv:" << sv << endl;
}
if (verbose) {
cout << "Computing num val of complex/real generators... " << endl;
}
computeXval();
}
bool CMfield::findRealStructure() {
xvi = convert_interval(genasrootof,nd,&ct);
if (verbose) {
cout << "xvi=" << xvi << endl;
cout << "conj(xvi)=" << conj(xvi,&ct) << endl;
}
if (!findInField(_subst(makesequence(minpol,x__IDNT_e,z__IDNT_e),&ct),conj(xvi,&ct),genconj)) {
cerr << "Trouble describing complex conjugation in field!" << endl;
exit(1);
}
if (verbose) {
cout << "Conjugate of generator: " << genconj << endl;
}
// cout << "gar: " << genasrootof << endl;
gen realgen1 = _simplify( (genasrootof + _subst(makesequence(genconj,x__IDNT_e,genasrootof),&ct)) ,&ct);
if (verbose) {
cout << "realgen1 = " << realgen1 << endl;
}
gen rmp1;
if (realgen1.type==0) {
rmp1 = s__IDNT_e-realgen1;
}
else {
rmp1 = expand(_poly2symb(makesequence(_pmin(eval(realgen1,1,&ct),&ct),s__IDNT_e),&ct),&ct);
}
if (!operator_equal(rmp1,s__IDNT_e,&ct) && operator_equal(_degree(makesequence(rmp1,s__IDNT_e),&ct),dc/2,&ct)) {
realgenasrootof = realgen1; // this is not right, will change it later
realminpol = rmp1;
realgen = _simplify(x__IDNT_e+genconj,&ct);
red(realgen);
if (verbose) {
cout << "Found a generator for real subfield: " << realgenasrootof << endl;
cout << "realminpol: " << rmp1 << endl;
cout << "realgen: " << realgen << endl;
cout << "realgen as rootof: " << realgenasrootof << endl;
}
}
else {
if (verbose) {
cout << "At first attempt, found degree " << _degree(makesequence(rmp1,s__IDNT_e),&ct) << endl;
}
gen realgen2 = _simplify((genasrootof*_subst(makesequence(genconj,x__IDNT_e,genasrootof),&ct)),&ct);
if (verbose) {
cout << "realgen2 = " << realgen2 << endl;
}
gen rmp2;
if (realgen2.type==0) {
rmp2 = s__IDNT_e-realgen2;
}
else {
rmp2 = giac::expand(_poly2symb(makesequence(_pmin(eval(realgen2,1,&ct),&ct),s__IDNT_e),&ct),&ct);
}
// cout << "rmp2 = " << rmp2 << endl;
if (!operator_equal(rmp2,s__IDNT_e,&ct) &&operator_equal(_degree(makesequence(rmp2,s__IDNT_e),&ct),dc/2,&ct)) {
realgenasrootof = realgen2; // this is not right, will change it later
realminpol = rmp2;
realgen = _simplify(x__IDNT_e*genconj,&ct);
red(realgen);
if (verbose) {
cout << "Found a generator for real subfield: " << realgenasrootof << endl;
cout << "realminpol: " << rmp2 << endl;
cout << "realgen: " << realgen << endl;
cout << "realgen as rootof: " << realgenasrootof << endl;
}
}
else {
cout << "At second attempt, found degree " << _degree(makesequence(rmp2,s__IDNT_e),&ct) << endl;
gen realgen3 = _simplify(((genasrootof+1)*(_subst(makesequence(genconj,x__IDNT_e,genasrootof),&ct)+1)),&ct);
if (verbose) {
cout << "realgen3 = " << realgen3 << endl;
}
gen rmp3;
if (realgen3.type==0) {
rmp3 = s__IDNT_e-realgen3;
}
else {
rmp3 = giac::expand(_poly2symb(makesequence(_pmin(eval(realgen3,1,&ct),&ct),s__IDNT_e),&ct),&ct);
}
// cout << "rmp2 = " << rmp2 << endl;
if (!operator_equal(rmp3,s__IDNT_e,&ct) && operator_equal(_degree(makesequence(rmp3,s__IDNT_e),&ct),dc/2,&ct)) {
realgenasrootof = realgen3; // this is not right, will change it later
realminpol = rmp3;
realgen = _simplify((x__IDNT_e+1)*(genconj+1),&ct);
red(realgen);
if (verbose) {
cout << "Found a generator for real subfield: " << realgenasrootof << endl;
cout << "realminpol: " << rmp3 << endl;
cout << "realgen: " << realgen << endl;
cout << "realgen as rootof: " << realgenasrootof << endl;
}
}
else {
cout << "At third attempt, found degree " << _degree(makesequence(rmp3,s__IDNT_e),&ct) << endl;
cerr << "Trouble finding generator for real subfield" << endl;
exit(1);
}
}
}
dr = _degree(makesequence(realminpol,s__IDNT_e),&ct);
for (int k=0; is_greater(dr,gen(k),&ct); k++) {
realminpolcoeffs.push_back(_coeff(makesequence(realminpol,s__IDNT_e,dr-k),&ct));
}
if (!2*dr==dc) {
cout << "This does not look like a CM field, degrees do not match!" << endl;
exit(1);
}
unsigned n = dr.val;
vars = vecteur(n);
for (unsigned j=0; j<n; ++j) {
vars[j] = gen("a"+print_INT_(n-j-1),&ct);
}
genericreal = _subst(makesequence(_poly2symb(makesequence(vars,s__IDNT_e),&ct),s__IDNT_e,realgen),&ct);
red(genericreal);
xv = _evalf(genasrootof,&ct);
if (verbose) {
cout << "xv:" << xv << endl;
}
sv = re(_evalf(realgenasrootof,&ct),&ct);
if (verbose) {
cout << "sv:" << sv << endl;
cout << "Computing num val of complex/real generators... " << endl;
}
computeXval();
}
bool CMfield::myCheckEqual(gen a0, gen a1, bool &res) {
gen pol0 = _pmin(a0,&ct);
if (verbose) {
cout << "pol0=" << pol0 << endl;
}
gen pol1 = _pmin(a1,&ct);
if (verbose) {
cout << "pol1=" << pol1 << endl;
}
if (operator_equal(pol0,pol1,&ct)) {
if (verbose) {
cout << "Same minimal polynomial.." << endl;
cout << "a0=" << a0 << " will convert to interval" << endl;
}
gen val0 = convert_interval(a0,nd,&ct);
if (verbose) {
cout << "val0=" << val0 << endl;
cout << "a0=" << a1 << " will convert to interval" << endl;
}
gen val1 = convert_interval(a1,nd,&ct);
if (verbose) {
cout << "val1=" << val1 << endl;
}
gen val0_re = re(val0,&ct);
if (operator_equal(val0_re,zeronumber,&ct)) {
val0_re = convert_interval(zeronumber,nd,&ct);
}
gen val0_im = im(val0,&ct);
if (operator_equal(val0_im,zeronumber,&ct)) {
val0_im = convert_interval(zeronumber,nd,&ct);
}
gen val1_re = re(val1,&ct);
if (operator_equal(val1_re,zeronumber,&ct)) {
val1_re = convert_interval(zeronumber,nd,&ct);
}
gen val1_im = im(val1,&ct);
if (operator_equal(val1_im,zeronumber,&ct)) {
val1_im = convert_interval(zeronumber,nd,&ct);
}
gen sols = _complexroot(makesequence(pol0,prec),&ct);
int count0 = 0;
int count1 = 0;
for (int k=0; k<(*sols._VECTptr).size(); k++) {
gen root = sols[k][0];
if (_intersect(makesequence(val0_re,re(root,&ct)),&ct).type==3) {
if (_intersect(makesequence(val0_im,im(root,&ct)),&ct).type==3) {
count0++;
}
}
if (_intersect(makesequence(val1_re,re(root,&ct)),&ct).type==3) {
if (_intersect(makesequence(val1_im,im(root,&ct)),&ct).type==3) {
count1++;
}
}
}
if (count0==0 || count1==0) {
res = false;
return true;
}
else {
if (count0>1 || count1>1) {
return false;
}
else {
res = true;
return true;
}
}
}
else {
res = false;
return true;
}
}
bool CMfield::extend(gen a) {
//
// Modifies minpol if needed, in order for a to be expressible in the number field
// if this works, b will contain a polynomial expression for a in terms of genasrootof
//
// the procedure returns true if an extension was needed
gen q = giac::expand(_poly2symb(makesequence(_pmin(normal(a,&ct),&ct),z__IDNT_e),&ct),&ct);
/* cout << "q=" << q << endl;
cout << "genasrootof=" << genasrootof << endl;*/
gen fa;
if (operator_equal(simplify(genasrootof-gen("i",&ct),&ct),zeronumber,&ct)) {
cout << "rootof was replaced by i, I will bypass usual factorization over number field" << endl;
vecteur fakefa;
fakefa.push_back(q);
fakefa.push_back(onenumber);
fa = gen(fakefa);
}
else {
if (operator_equal(simplify(genasrootof-onenumber,&ct),zeronumber,&ct)) {
if (verbose) {
cout << "Initial field seems to be Q, I will bypass usual factorization over number field" << endl;
}
vecteur fakefa;
fakefa.push_back(q);
fakefa.push_back(onenumber);
fa = gen(fakefa);
}
else {
fa = _factors(makesequence(q,genasrootof),&ct);
}
}
if (verbose) {
cout << "" << endl;
}
// cout << "fa=" << fa << endl;
int l = 0;
while (2*l<(*fa._VECTptr).size()) {
// cout << "l=" << l << endl;
int deg = _degree(makesequence(fa[2*l],z__IDNT_e),&ct).val;
if (deg==1) {
gen a0 = normal(-_coeff(makesequence(fa[2*l],z__IDNT_e,0),&ct),&ct);
gen pol0 = _pmin(a0,&ct);
if (verbose) {
cout << "a0=" << a0 << ", pmin :" << pol0 << endl;
}
gen a1 = normal(a*_coeff(makesequence(fa[2*l],z__IDNT_e,1),&ct),&ct);
gen pol1 = _pmin(a1,&ct);
if (verbose) {
cout << "a1=" << a1 << ", pmin :" << pol1 << endl;
}
if (verbose) {
cout << "My check? " << endl;
}
bool test;
bool worked = myCheckEqual(a0,a1,test);
if (worked) {
if (test) {
return false;
}
}
if (verbose) {
cout << "Done with mycheck" << endl;
}
// gen te = _simplify(a0-a1,&ct);
gen te = normal(eval(a0-a1,&ct),&ct);
// cout << "te=" << te << endl;
// cout << "pmin(te)=" << _pmin(te,&ct) << endl;
if (operator_equal(te,zeronumber,&ct)) {
return false;
}
else {
if (is_greater(1e-5,abs(_evalf(te,&ct)),&ct)) {
cerr << "Supposedly non-zero: " << _evalf(te,&ct) << endl;
exit(1); // may want to remove this eventually...
}
}
}
else {
if (deg>1) {
// cout << "deg=" << deg << endl;
gen te = normal(normal(_subst(makesequence(fa[2*l],z__IDNT_e,a),&ct),&ct),&ct);
/* cout << "a=" << a << endl;
cout << "te=" << te << endl;
cout << "evalf(te)=" << _evalf(te,&ct) << endl;*/
if (operator_equal(te,zeronumber,&ct)) {
// need to extend, will use RUR
if (verbose) {
cout << "Splitting this factor: " << fa[2*l] << endl;
}
vecteur coe;
for (int uu=deg; uu>=0; uu--) {
gen teuu = _coeff(makesequence(fa[2*l],z__IDNT_e,uu),&ct);
if (contains(teuu,gen("j",&ct))) {
coe.push_back(_subst(makesequence(teuu,gen("j",&ct),x__IDNT_e),&ct));
}
else {
if (teuu.type==8) {
if (teuu.is_symb_of_sommet(at_prod)) {
gen yoo = _poly2symb(makesequence(_simplify(teuu[1][1]*teuu[2],&ct),x__IDNT_e),&ct);
coe.push_back(yoo);
}
else {
if (teuu.is_symb_of_sommet(at_rootof)) {
gen yoo = _poly2symb(makesequence(_simplify(teuu[1],&ct),x__IDNT_e),&ct);
coe.push_back(yoo);
}
else {
cout << "Problem reading factors!" << endl;
}
}
}
else {
coe.push_back(teuu);
}
}
}
vecteur eqs;
eqs.push_back(minpol);
eqs.push_back(giac::expand(_poly2symb(makesequence(coe,a__IDNT_e),&ct),&ct));
vecteur vars;
vars.push_back(x__IDNT_e);
vars.push_back(a__IDNT_e);
if (verbose) {
cout << "eqs: " << eqs << endl;
cout << "vars: " << vars << endl;
}
gen gb = _gbasis(makesequence(eqs,vars,change_subtype(_RUR_REVLEX,_INT_GROEBNER)),&ct);
if (gb[0]==-4) {
gen varprov = lidnt(gb[2])[0];
if (verbose) {
cout << "varprov=" << varprov << endl;
cout << "gb[2]=" << gb[2] << endl;
}
minpol = _subst(makesequence(gb[2],varprov,x__IDNT_e),&ct);
if (operator_equal(_coeff(makesequence(minpol,x__IDNT_e,_degree(makesequence(minpol,x__IDNT_e),&ct).val),&ct),-onenumber,&ct)) {
// cout << "Changing sign, polynomial is not monic" << endl;
minpol = -minpol;
}
// still need to locate the right root, find complex conjugate, real generator
if (verbose) {
cout << "New minimal polynomial: " << minpol << endl;
}
// if minpol is x^2+1, sto will not work!
gen tepo = simplify(minpol - gen("x^2+1",&ct),&ct);
if (!operator_equal(tepo,zeronumber,&ct)) {
gen fafa = _factors(minpol,&ct);
int nfafa = (*fafa._VECTptr).size();
bool foundfactor = false;
int m = 0;
while (!foundfactor && 2*m<nfafa) {
if (_degree(makesequence(fafa[2*m],x__IDNT_e),&ct).val>0) {
gen pa = rootof(fafa[2*m],&ct);
gen subfa = _factors(makesequence(fafa[2*m],pa),&ct);
int nsubfa = (*subfa._VECTptr).size();
int sn=0;
while (!foundfactor && 2*sn<nsubfa) {
gen ro = normal(-_coeff(makesequence(subfa[2*sn],x__IDNT_e,0),&ct)/_coeff(makesequence(subfa[2*sn],x__IDNT_e,1),&ct),&ct);
gen xro = normal(genasrootof-_subst(makesequence(gb[4],varprov,ro),&ct)/_subst(makesequence(gb[3],varprov,ro),&ct),&ct);
gen aro = normal(a-_subst(makesequence(gb[5],varprov,ro),&ct)/_subst(makesequence(gb[3],varprov,ro),&ct),&ct);
foundfactor = (operator_equal(xro,zeronumber,&ct) && operator_equal(aro,zeronumber,&ct));
sn++;
}
if (foundfactor) {
minpol = fafa[2*m];
genasrootof = rootof(minpol,&ct);
dc = _degree(fafa[2*m],&ct);
}
}
m++;
}
if (!foundfactor) {
cout << "Trouble extending" << endl;
}
/* if (verbose) {
cout << "Sto j with " << genasrootof << endl;
}
gen jsto("j",&ct);
sto(jsto,genasrootof,&ct);*/
gen tepo = simplify(minpol - gen("x^2+1",&ct),&ct);
if (!operator_equal(tepo,zeronumber,&ct)) {
genasrootof = rootof(minpol,&ct);
if (verbose) {
cout << "Sto j with " << genasrootof << endl;
}
gen jsto("j",&ct);
sto(jsto,genasrootof,&ct);
}
else {
genasrootof = rootof(minpol,&ct);
}
return true;
}
else {
genasrootof = rootof(minpol,&ct);
dc = _degree(minpol,&ct);
return true;
}
}
else {
cout << "gb=" << gb << endl;
cerr << "Trouble computing primitive element" << endl;
exit(1);
}
}
else {
cout << a << " is not a root of this factor, ignoring " << fa[2*l] << endl;
}
}
/* else {
cout << "deg=" << deg << ", is this zero? if so don't want to do anything" << endl;
}*/
// else deg==0, we don't do anything
}
l++;
}
// if we get here there is a problem!
cerr << "Trouble extending!" << endl;
exit(1);
}
void CMfield::increasePrecision() {
nd = nd+pf;
if (nd>ndrecord) {
ndrecord = nd;
}
prec = _inverse(_pow(makesequence(ten,(int)nd),&ct),&ct);
computeXval();
}
void CMfield::lowerPrecision() {
prec = _inverse(_pow(makesequence(ten,(int)nd),&ct),&ct);
computeXval(); // not clear whether we need to do this..
}
bool CMfield::safeSift(gen &gb, vecteur &res) {
/*
* I am assuming gb is a RUR Groebner basis, with parameter named s
* (but this does not mean values of realminpol!)
*/
cout << "safeSift" << endl;
unsigned savend = nd;
bool doneSifting = false;
while (!doneSifting && nd<ndmax) {
doneSifting = siftAttempt(gb,res);
if (!doneSifting) {
increasePrecision();
}
}
nd = savend;
lowerPrecision();
// cout << "done safesift" << endl;
// cout << "done safesift, res= " << res << endl;
return doneSifting;
}
bool CMfield::siftAttempt(gen &gb, vecteur &res) {
gen sols = _realroot(makesequence(gb[2],prec),&ct);
// cout << "sols=" << sols << endl;
gen var = lidnt(gb[2]).front();
res.clear();
unsigned nbsols = (*sols._VECTptr).size();
for (unsigned k=0; k<nbsols; k++) {
gen den = _horner(makesequence(gb[3],sols[k][0],var),&ct);
gen svk = convert_interval(_horner(makesequence(gb[4],sols[k][0],var),&ct)/den,nd,&ct);
// cout << "svi=" << svi << endl;
// cout << "svk=" << svk << endl;
gen J = _intersect(makesequence(svi,svk),&ct);
// cout << "J=" << J << endl;
if (J.type==3) {
// To be safe, we will also check that the interval does not contain any other root!
bool isambiguous = false;
unsigned l=0;
unsigned lmax = (*othersvals._VECTptr).size();
while (!isambiguous && l<lmax) {
gen Jl = _intersect(makesequence(othersvals[l],svk),&ct);
isambiguous = (Jl.type==3);
l++;
}
if (isambiguous) {
return false;
}
else {
gen rt = convert_interval(sols[k][0],nd,&ct);
/* gen rt;
if (sols[k][0].type==8) {
cout << "Converting real root " << sols[k][0] << " to interval" << endl;
rt = convert_interval(sols[k][0],nd,&ct);
cout << " result " << rt << endl;
}
else {
cout << " NOT converting real root " << sols[k][0] << " to interval" << endl;
rt = sols[k][0];
}*/
res.push_back(rt);
}
}
}
return true;
}
bool CMfield::identifyQuadraticNumber(gen &te, vecteur &res) {
if (te.is_symb_of_sommet(at_prod)) {
if (te[1].is_symb_of_sommet(at_plus)) {
if (te[1][1].is_symb_of_sommet(at_neg)) {
if (te[1][1][1].is_symb_of_sommet(at_prod)) {
if (te[1][1][1][2].is_symb_of_sommet(at_pow)) {
res.push_back(te[1][1][1][2][1]);
res.push_back(_simplify(-te[1][1][1][1]*te[2],&ct));
res.push_back(_simplify(te[1][2]*te[2],&ct));
return true;
}
else {
cout << "Don't know what to do (1)..." << endl;
return false;
}
}
else {
// cout << "no *" << endl;
if (te[1][1][1].is_symb_of_sommet(at_pow)) {
res.push_back(te[1][1][1][1]);
res.push_back(_simplify(-te[2],&ct));
res.push_back(_simplify(te[1][2]*te[2],&ct));
return true;
}
else {
cout << "Don't know what to do (2)..." << endl;
return false;
}
}
}
else {
if (te[1][1].is_symb_of_sommet(at_prod)) {
if (te[1][1][2].is_symb_of_sommet(at_pow)) {
res.push_back(te[1][1][2][1]);
res.push_back(_simplify(te[1][1][1]*te[2],&ct));
res.push_back(_simplify(te[1][2]*te[2],&ct));
return true;
}
}
else {
cout << "te[1][1][0]: " << te[1][1][0] << endl;
cout << "Don't know what to do (3)... " << endl;
return false;
}
}
}
else {
cout << "Don't know what to do (4)..." << endl;
return false;
}
}
else {
if (te.is_symb_of_sommet(at_neg)) {
// could have a minus sign in front, say -a*sqrt(?)
if (te[1].is_symb_of_sommet(at_prod)) {
if (te[1][2].is_symb_of_sommet(at_pow)) {
res.push_back(te[1][2][1]);
res.push_back(-te[1][1]);
res.push_back(zeronumber);
return true;
}
else {
return false;
cout << "Don't know what to do (5)..." << endl;
}
}
else {
return false;
cout << "Don't know what to do (6)..." << endl;
}
}
else {
// need to double check this one
if (te.is_symb_of_sommet(at_prod)) {
if (te[2].is_symb_of_sommet(at_pow)) {
res.push_back(te[2][1]);
res.push_back(te[1]);
res.push_back(zeronumber);
return true;
}
else {
cout << "Don't know what to do (7)..." << endl;
return false;
}
}
else {
if (te.is_symb_of_sommet(at_plus)) {
if (te[1].is_symb_of_sommet(at_neg)) {
if (te[1][1].is_symb_of_sommet(at_prod)) {
if (te[1][1][2].is_symb_of_sommet(at_pow)) {
res.push_back(te[1][1][2][1]);
res.push_back(-te[1][1][1]);
res.push_back(te[2]);
return true;
}
else {
cout << "Don't know what to do (8)..." << endl;
return false;
}
}
else {
if (te[1][1].is_symb_of_sommet(at_pow)) {
res.push_back(te[1][1][1]);
res.push_back(-1);
res.push_back(te[2]);
return true;
}
else {
cout << "Don't know what to do (9)..." << endl;
return false;
}
}
}
}
}
}
}
}
gen CMfield::getCoeffs(gen &te) {
if (te.is_symb_of_sommet(at_prod)) {
return _simplify(te[1][1]*te[2],&ct);
}
else {
return te[1];
}
}
gen CMfield::getCoeffs(gen te, gen prec) {
/*
* not sure what this will apply to, probably want the parameter of the rootof to be minpol?
*/
if (te.is_symb_of_sommet(at_prod)) {
gen subfa = te._SYMBptr->feuille;
gen yo = _simplify(te[1][1]*te[2],&ct);
int umax = dc.val-(*yo._VECTptr).size();
for (int u=0; u<umax; u++) {
(*yo._VECTptr).insert((*yo._VECTptr).begin(),zeronumber);
}
gen ls = _inverse(genasrootofpowers,&ct)*yo;
return _poly2symb(makesequence(ls,x__IDNT_e),&ct);
}
else {
if (te.is_symb_of_sommet(at_rootof)) {
gen yo = _simplify(te[1],&ct);
int umax = dc.val-(*yo._VECTptr).size();
for (int u=0; u<umax; u++) {
(*yo._VECTptr).insert((*yo._VECTptr).begin(),zeronumber);
}
gen ls = _inverse(genasrootofpowers,&ct)*yo;
return _poly2symb(makesequence(ls,x__IDNT_e),&ct);
}
else {
cout << "Problem reading factors!" << endl;
}
}
}
bool CMfield::findInField(gen g, gen approx, gen &res) {
unsigned savend = nd;
bool done = false;
while (!done && nd<ndmax) {
done = findInFieldAttempt(g,approx,res);
if (!done) {
nd = nd+pf;
if (nd>ndrecord) {
ndrecord = nd;
}
cout << "Increasing precision, nd=" << nd << endl;
prec = _inverse(_pow(makesequence(ten,(int)nd),&ct),&ct);
computeXval();
}
}
nd = savend;
prec = _inverse(_pow(makesequence(ten,(int)nd),&ct),&ct);
computeXval();
return done;
}
bool CMfield::findInFieldAttempt(gen g, gen approx, gen &res) {
vecteur variables = lidnt(g);
if (variables.size()!=1) {
cout << "Polynomial is not univariate?" << endl;
}
else {
gen fa = _factors(makesequence(g,genasrootof),&ct);
// gen fa = _factors(makesequence(eval(g,&ct),genasrootof),&ct);
int nf = (*fa._VECTptr).size()/2;
int countsols = 0;
for (int k=0; k<nf; k++) {
if (verbose) {
cout << "fa=" << fa[2*k] << endl;
cout << "variable:" << variables[0] << endl;
}
gen df = _degree(makesequence(fa[2*k],variables[0]),&ct);
if (df==1) {
gen yo;
gen te = normal(-_coeff(makesequence(fa[2*k],variables[0],0),&ct)/_coeff(makesequence(fa[2*k],variables[0],1),&ct),&ct);
if (operator_equal(te-genasrootof,zeronumber,&ct)) {
yo = x__IDNT_e;
}
else {
if (operator_equal(te+genasrootof,zeronumber,&ct)) {
yo = -x__IDNT_e;
}
else {
// gen pol_te = giac::expand(_poly2symb(makesequence(_pmin(eval(te,&ct),&ct),z__IDNT_e),&ct),&ct);
gen pol_te = giac::expand(_poly2symb(makesequence(_pmin(te,&ct),z__IDNT_e),&ct),&ct);
// cout << "pol_te=" << pol_te << endl;
if (contains(te,gen("j",&ct))) {
yo = _subst(makesequence(te,gen("j",&ct),x__IDNT_e),&ct);
}
else {
if (contains(te,gen("jr",&ct))) {
yo = _subst(makesequence(te,gen("jr",&ct),realgen),&ct);
}
else {
// cout << "te=" << te << endl;
if (te.is_symb_of_sommet(at_prod)) { // this is not necessarily correct (problem when -sign for instance)
// yo = normal(te[1][1]*te[2],&ct);
yo = _simplify(te[1][1]*te[2],&ct);
}
else {
if (te.is_symb_of_sommet(at_rootof)) {
yo = te[1];
}
else {
yo = te; // this is not great, but it's all I can do!
}
}
yo = giac::expand(_poly2symb(makesequence(yo,x__IDNT_e),&ct),&ct);
}
}
}
}
gen tei;
if (contains(yo,x__IDNT_e)) {
tei = _horner(makesequence(yo,xvi,x__IDNT_e),&ct);
}
else {
tei = convert_interval(yo,nd,&ct);
}
/* cout << "approx=" << approx << endl;
cout << "tei=" << tei << endl;*/
gen tei_r = re(convert_interval(tei,nd,&ct),&ct);
gen approx_r = re(convert_interval(approx,nd,&ct),&ct);
/* cout << "approx_r=" << approx_r << endl;
cout << "tei_r=" << tei_r << endl;*/
gen Jr = _intersect(makesequence(tei_r,approx_r),&ct);
if (Jr.type==3){
gen tei_i = convert_interval(im(convert_interval(tei,nd,&ct),&ct),nd,&ct);
gen approx_i = convert_interval(im(convert_interval(approx,nd,&ct),&ct),nd,&ct);
/* cout << "approx_i=" << approx_i << endl;
cout << "tei_i=" << tei_i << endl;*/
gen Ji = _intersect(makesequence(tei_i,approx_i),&ct);
if (Ji.type==3){
// cout << "intersect" << endl;
res = yo;
countsols++;
}
}
/* else {
cout << "not intersecting?" << endl;
}*/
}
else {
if (df>1) {
cout << "Polynomial doesn't split over nb field" << endl;
}
}
}
if (countsols!=1) {
cout << "Trouble, found " << countsols << " values! " << endl;
return false;
}
else {
return true;
}
}
}
void CMfield::findgenasrootof() {
gen xtest("xx",&ct);
gen rtx=rootof(minpol,&ct);
sto(xtest,rtx,&ct);
gen fa = _factors(makesequence(minpol,rtx),&ct);
int nf = (*fa._VECTptr).size()/2;
int k=0;
bool found = false;
while (k<nf && !found) {
gen df = _degree(makesequence(fa[2*k],x__IDNT_e),&ct);
if (df==1) {
gen te = _simplify(-_coeff(makesequence(fa[2*k],x__IDNT_e,0),&ct)/_coeff(makesequence(fa[2*k],x__IDNT_e,1),&ct),&ct);
found = 1e-8 > abs(_evalf(te,&ct)-xv);
if (found) {
genasrootof = te;
cout << "Gen as rootof: " << genasrootof << endl;
cout << " check: " << _evalf(genasrootof,&ct) << endl;
cout << " " << _evalf(xvi,&ct) << endl;
int m = dc.val;
vecteur gpv;
vecteur col0;
for (unsigned k=0; k<m-1; k++) {
col0.push_back(zeronumber);
}
col0.push_back(onenumber);
for (unsigned k=0; k<m-1; k++) {
if (te.is_symb_of_sommet(at_prod)) {
gen subfa = te._SYMBptr->feuille; // unused below!!
gen yo = _simplify(te[1][1]*te[2],&ct);
int umax = m-(*yo._VECTptr).size();
for (int u=0; u<umax; u++) {
(*yo._VECTptr).insert((*yo._VECTptr).begin(),zeronumber);
}
gpv.insert(gpv.begin(),yo);
}
else {
if (te.is_symb_of_sommet(at_rootof)) {
gen yo = _simplify(te[1],&ct);
int umax = m-(*yo._VECTptr).size();
for (int u=0; u<umax; u++) {
(*yo._VECTptr).insert((*yo._VECTptr).begin(),zeronumber);
}
gpv.insert(gpv.begin(),yo);
}
else {
cout << "Problem with rootof powers!" << endl;
}
}
te = _simplify(te*genasrootof,&ct);
}
gpv.push_back(col0);
genasrootofpowers = _tran(gen(gpv),&ct);
}
}
k++;
}
if (!found) {
cout << "Problem computing field generator as rootof!" << endl;
}
else {
gen stest("ss",&ct);
gen rt=rootof(realminpol,&ct);
sto(stest,rt,&ct);
fa = _factors(makesequence(realminpol,rt),&ct);
nf = (*fa._VECTptr).size()/2;
k=0;
found = false;
while (k<nf && !found) {
gen df = _degree(makesequence(fa[2*k],s__IDNT_e),&ct);
if (operator_equal(df,onenumber,&ct)) {
gen te = _simplify(-_coeff(makesequence(fa[2*k],s__IDNT_e,0),&ct)/_coeff(makesequence(fa[2*k],s__IDNT_e,1),&ct),&ct);
found = is_greater(1e-8,abs(_evalf(te,&ct)-svi),&ct);
if (found) {
realgenasrootof = te;
cout << "Real gen as rootof: " << realgenasrootof << endl;
cout << " check: " << _evalf(realgenasrootof,&ct) << endl;
cout << " " << _evalf(svi,&ct) << endl;
}
}
k++;
}
if (!found) {
cout << "Problem computing real subfield generator as rootof!" << endl;
}
}
}
gen CMfield::convertToReal(gen g) {
gen eqs = _symb2poly(makesequence(genericreal-g,x__IDNT_e),&ct);
gen sol = _solve(makesequence(eqs,vars),&ct);
if ((*sol._VECTptr).size()==1) {
return _simplify(_poly2symb(makesequence(sol[0],s__IDNT_e),&ct),&ct);
}
else {
cerr << "Error converting to real generator, is this number real??" << endl;
exit(1);
}
}
void CMfield::computeXval() {
xvi = convert_interval(genasrootof,nd,&ct); // this does not always produce the same value!!
// cout << "xvi=" << xvi << endl;
svi = re(convert_interval(_horner(makesequence(realgen,xvi,x__IDNT_e),&ct),nd,&ct),&ct);
// cout << "svi=" << svi << endl;
gen rtx = _complexroot(makesequence(minpol,prec),&ct);
int counter = 0;
int k = 0;
int kmax = (*rtx._VECTptr).size();
while (k<kmax) {
gen valk = convert_interval(rtx[k][0],nd,&ct);
// cout << "comparing " << valk << " and " << xvi << endl;
gen Jt = _intersect( makesequence(re(valk,&ct),re(xvi,&ct)), &ct );
if (Jt.type==3) {
Jt = _intersect(makesequence(im(valk,&ct),im(xvi,&ct)),&ct);
if (Jt.type==3) {
xvi = valk;
counter++;
}
}
k++;
}
if (counter!=1) {
cerr << "Trouble with generator value, xvi is not isolating (should increase precision), counter=" << counter << endl;
exit(1);
}
gen rts = _realroot(makesequence(realminpol,prec),&ct);
vecteur osv;
counter = 0;
k = 0;
kmax = (*rts._VECTptr).size();
while (k<kmax) {
gen te = convert_interval(rts[k][0],nd,&ct);
// cout << "te= " << te << endl;
gen J = _intersect(makesequence(svi,te),&ct); // rts[k][0] may not be an interval!
// cout << "J= " << J << endl;
if (J.type==3) {
svi = te;
// svi = rts[k][0];
// cout << "svi =" << svi << endl;
counter++;
}
else {
osv.push_back(te);
}
k++;
}
if (counter!=1) {
// cout << "Found " << counter << " values for s (so far)" << endl;
// cout << " (this should not be a problem)" << endl;
}
else {
othersvals = gen(osv);
}
}
gen CMfield::createMatrix(int m, int n, vecteur v) {
vecteur te;
for (int j=0; j<m; j++) {
vecteur tej;
for (int k=0; k<n; k++) {
tej.push_back(v[m*j+k]);
}
te.push_back(tej);
}
return gen(te);
}
gen CMfield::createMatrix(vecteur v) {
vecteur r1 = vecteur(3);
vecteur r2 = vecteur(3);
vecteur r3 = vecteur(3);
for (int j=0; j<3; j++) {
r1[j] = v[j];
}
for (int j=0; j<3; j++) {
r2[j] = v[3+j];
}
for (int j=0; j<3; j++) {
r3[j] = v[6+j];
}
vecteur res(3);
res[0] = gen(r1);
res[1] = gen(r2);
res[2] = gen(r3);
return gen(res);
}
gen CMfield::mypower(gen M, int k) {
gen res = Id;
for (int j=0; j<k; j++) {
res = res*M;
matred(res);
}
return res;
}
int CMfield::order(gen M, int k) {
bool done = isScalar(M);
gen te = M;
int j=1;
if (done) {
return j;
}
while (!done && j<k) {
j++;
te = te * M;
matred(te);
done = isScalar(te);
}
if (done) {
return j;
}
else {
return -1;
}
}
int CMfield::order(gen M) {
int k = order(M,1000);
if (k>-1) {
return k;
}
else {
// cout << "Order is > 1000" << endl;
return -1;
}
}
int CMfield::linearorder(gen M, int k) {
bool done = operator_equal(M,Id,&ct);
gen te = M;
int j=1;
if (done) {
return j;
}
while (!done && j<k) {
j++;
te = te * M;
matred(te);
done = operator_equal(te,Id,&ct);
}
if (done) {
return j;
}
else {
return -1;
}
}
int CMfield::linearorder(gen M) {
int k = linearorder(M,1000);
if (k>-1) {
return k;
}
else {
// cout << "Order is > 1000" << endl;
return -1;
}
}
int CMfield::linearorder2d(gen M, int k) {
// bool done = isSameMatrix(M,Id2);
// cout << "Computing linear order of " << M << endl;
bool done = operator_equal(M,Id2,&ct);
gen te = M;
int j=1;
while (!done && j<k) {
j++;
te = te * M;
matred(te);
done = operator_equal(te,Id2,&ct);
// cout << "te=" << te << endl;
// cout << "Id2=" << Id2 << endl;
// cout << "Equal? " << done << endl;
}
if (done) {
return j;
}
else {
return -1;
}
}
int CMfield::linearorder2d(gen M) {
// int k =
return linearorder2d(M,1000);
/* if (k>-1) {
return k;
}
else {
// cout << "Order is > 1000" << endl;
return -1;
}*/
}
int CMfield::braidOrder(gen M, gen N) {
int k = braidOrder(M,N,1000);
if (k>-1) {
return k;
}
else {
// cout << "Order is > 1000" << endl;
return -1;
}
}
int CMfield::braidOrder2d(gen M, gen N) {
int k = braidOrder2d(M,N,1000);
if (k>-1) {
return k;
}
else {
// cout << "Order is > 1000" << endl;
return -1;
}
}
bool CMfield::areDep(gen v1, gen v2) {
bool allzero = true;
int c = 0;
while (allzero && c<3) {
gen te = v1[c]*v2[(c+1)%3] - v2[c]*v1[(c+1)%3];
red(te);
allzero = operator_equal(te,zeronumber,&ct);
c++;
}
return allzero;
}
bool CMfield::isSameMatrix(gen M, gen N) {
gen v = _row(makesequence(M,0),&ct);
int m = (*v._VECTptr).size();
bool same = true;
int j = 0;
while (same and j<m) {
int k = 0;
while (same and k<m) {
gen te = M[j][k]-N[j][k];
red(te);
same = operator_equal(te,zeronumber,&ct);
k++;
}
j++;
}
return same;
}
bool CMfield::isMultiple(gen M, gen N) {
// starts by finding a nonzero element in M, then uses it as a "pivot"
bool fdnz = false;
int k=-1;
int l;
int d = (*M._VECTptr).size();
int dm1 = d-1;
while (!fdnz && k!=dm1) {
k++;
l=-1;
while (!fdnz && l!=dm1) {
l++;
fdnz = !operator_equal(M[k][l],zeronumber,&ct);
if (!fdnz && !operator_equal(N[k][l],zeronumber,&ct)) {
return false;
}
}
}
if (!fdnz) {
cerr << "Matrix is zero, this should not happen in a group!" << endl;
exit(1);
}
else {
// cout << "fdnz" << endl;
// cout << "k=" << k << ", l=" << l << endl;
gen a = M[k][l];
gen b = N[k][l];
bool ismult = true;
while (ismult && l<d) {
gen te = b*M[k][l]-a*N[k][l];
red(te);
ismult = operator_equal(te,zeronumber,&ct);
l++;
}
k++;
while (ismult && k<d) {
l=0;
while (ismult && l<d) {
// cout << "k=" << k << ", l=" << l << endl;
gen te = b*M[k][l]-a*N[k][l];
red(te);
ismult = operator_equal(te,zeronumber,&ct);
l++;
}
k++;
}
return ismult;
}
}
/*
bool CMfield::isMultiple(gen M, gen N) {
// not sure what this is supposed to be doing!
// it certainly does not check if matrices are multiples of each other...
bool res = true;
int c = 0;
while (res && c<3) {
res = areDep(M[c],N[c]);
c++;
}
return res;
}*/
int CMfield::braidOrder(gen M, gen N, int k) {
gen A = M;
gen B = N;
if (isMultiple(A,B)) {
return 1;
}
else {
int o = 2;
while (o<k) {
A = A*N;
B = B*M;
matred(A);
matred(B);
if (isMultiple(A,B)) {
return o;
}
else {
o++;
A = A*M;
B = B*N;
matred(A);
matred(B);
if (isMultiple(A,B)) {
return o;
}
else {
o++;
}
}
}
}
return -1;
}
int CMfield::braidOrder2d(gen M, gen N, int k) {
gen A = M;
gen B = N;
if (isSameMatrix(A,B)) {
return 1;
}
else {
int o = 2;
while (o<k) {
A = A*N;
B = B*M;
matred(A);
matred(B);
if (isSameMatrix(A,B)) {
return o;
}
else {
o++;
A = A*M;
B = B*N;
matred(A);
matred(B);
if (isSameMatrix(A,B)) {
return o;
}
else {
o++;
}
}
}
}
return -1;
}
bool CMfield::testDependent(gen V, gen W) {
bool res = false;
int k=0;
while (!res && k<3) {
int kp1 = (k+1)%3;
gen te = V[k]*W[kp1] - V[kp1]*W[k];
red(te);
res = operator_equal(te,zeronumber,&ct);
}
}
bool CMfield::isScalar(gen M) {
bool test = true;
int n = (*M._VECTptr).size();
int j = 0;
gen diag = M[0][0];
while (j<n && test) {
int k = 0;
while (k<n && test) {
if (j!=k) {
test = operator_equal(M[j][k],zeronumber,&ct);
}
else {
if (j>0) {
gen te = (M[j][j]-diag);
red(te);
test = operator_equal(te,zeronumber,&ct);
}
}
k++;
}
j++;
}
return test;
/* if (M[0][0]==M[1][1] && M[1][1]==M[2][2]) {
gen ze = gen("0",&ct);
return (M[0][1]==ze && M[0][2]==ze && M[1][2]==ze && M[1][0]==ze && M[2][0]==ze && M[2][1]==ze);
}
else {
return false;
}*/
}
void CMfield::red(gen &g) {
g = _simplify(g,&ct);
if (contains(g,s__IDNT_e)) {
if (contains(g,x__IDNT_e)) {
gen gx = _simplify(_subst(makesequence(g,s__IDNT_e,realgen),&ct),&ct);
g = _simplify(_rem(makesequence(gx,minpol,x__IDNT_e),&ct),&ct);
}
else {
g = _simplify(_rem(makesequence(g,realminpol,s__IDNT_e),&ct),&ct);
}
}
else {
g= _simplify(_rem(makesequence(g,minpol,x__IDNT_e),&ct),&ct);
}
}
void CMfield::vectred(gen &v) {
int n = (*v._VECTptr).size();
for (int j=0; j<n; j++) {
red((*v._VECTptr)[j]);
}
}
void CMfield::matred(gen &M) {
int m = (*M._VECTptr).size();
for (int j=0; j<m; j++) {
vectred((*M._VECTptr)[j]);
}
}
gen CMfield::myconj(gen g) {
gen gco = _subst(makesequence(g,x__IDNT_e,genconj),&ct);
red(gco);
return gco;
}
gen CMfield::vectconj(gen v) {
int n = (*v._VECTptr).size();
vecteur res(n);
for (int j=0; j<n; j++) {
res[j] = myconj(v[j]);
}
return gen(res);
}
gen CMfield::matconj(gen M) {
int m = (*M._VECTptr).size();
vecteur res(m);
for (int j=0; j<m; j++) {
res[j] = vectconj(M[j]);
}
return gen(res);
}
gen CMfield::det(gen M) {
gen te = _det(M,&ct);
red(te);
return te;
}
gen CMfield::trace(gen M) {
gen te = M[0][0]+M[1][1]+M[2][2];
red(te);
return te;
}
gen CMfield::discr(gen z) {
gen zc = myconj(z);
gen z2 = z*z;
gen z2c = myconj(z2);
gen z3 = z2*z;
gen z3c = myconj(z3);
gen te = z2*z2c - 4*(z3+z3c) + 18*z*zc - 27;
red(te);
return convertToReal(te);
}
gen CMfield::discrAdjusted(gen M) {
gen tee = det(M);
red(tee);
gen de = inverse(tee);
gen z = trace(M);
gen zc = myconj(z);
gen z2 = z*z;
gen z2c = myconj(z2);
gen z3 = z2*z*de;
gen z3c = myconj(z3);
gen te = z2*z2c - 4*(z3+z3c) + 18*z*zc - 27;
red(te);
return convertToReal(te);
}
bool CMfield::isReal(gen g) {
gen gt = g-myconj(g);
red(gt);
return gt==0;
}
gen CMfield::realPart(gen g) {
gen gt = (g+myconj(g))/2;
red(gt);
return convertToReal(gt);
}
gen CMfield::inverse(gen g) {
if (contains(g,s__IDNT_e)) {
if (contains(g,x__IDNT_e)) {
gen gx = _subst(makesequence(g,s__IDNT_e,realgen),&ct);
gen e =_egcd(makesequence(gx,minpol,x__IDNT_e),&ct);
gen res = e[0]/e[2];
red(res);
return res;
}
else {
gen e =_egcd(makesequence(g,realminpol,s__IDNT_e),&ct);
gen res = e[0]/e[2];
red(res);
return res;
}
}
else {
gen e =_egcd(makesequence(g,minpol,x__IDNT_e),&ct);
gen res = e[0]/e[2];
red(res);
return res;
}
}
gen CMfield::evali(gen g) {
if (contains(g,s__IDNT_e)) {
if (contains(g,x__IDNT_e)) {
gen gx = _subst(makesequence(g,s__IDNT_e,realgen),&ct);
return _subst(makesequence(gx,x__IDNT_e,xvi),&ct);
}
else {
return _subst(makesequence(g,s__IDNT_e,svi),&ct);
}
}
else {
// cout << "no s" << endl;
if (contains(g,x__IDNT_e)) {
return _subst(makesequence(g,x__IDNT_e,xvi),&ct);
}
else {
// cout << "no x either" << endl;
// this will work only for numbers?!
// cout << "g=" << g << endl;
if (operator_equal(g,zeronumber,&ct)) {
return convert_interval(onenumber,nd,&ct)-onenumber;
}
else {
// cout << "no s no x" << endl;
// cout << convert_interval(g,nd,&ct) << endl;
return convert_interval(g,nd,&ct);
}
}
}
}
gen CMfield::vectevali(gen v) {
vecteur te;
int n = (*v._VECTptr).size();
for (int j=0; j<n; j++) {
te.push_back(evali((*v._VECTptr)[j]));
}
gen tee = gen(te);
return tee;
}
gen CMfield::matevali(gen M) {
vecteur te;
int m = (*M._VECTptr).size();
for (int j=0; j<m; j++) {
te.push_back(vectevali((*M._VECTptr)[j]));
}
gen tee = gen(te);
return tee;
}
gen CMfield::evaln(gen g) {
if (contains(g,s__IDNT_e)) {
if (contains(g,x__IDNT_e)) {
gen gx = _subst(makesequence(g,s__IDNT_e,realgen),&ct);
return _subst(makesequence(gx,x__IDNT_e,xv),&ct);
}
else {
return _subst(makesequence(g,s__IDNT_e,sv),&ct);
}
}
else {
if (contains(g,x__IDNT_e)) {
return _subst(makesequence(g,x__IDNT_e,xv),&ct);
}
else {
if (operator_equal(g,zeronumber,&ct)) {
return g;
}
else {
return _evalf(makesequence(g,(int)nd),&ct);
}
}
}
}
gen CMfield::vectevaln(gen v) {
vecteur te;
int n = (*v._VECTptr).size();
for (int j=0; j<n; j++) {
te.push_back(evaln((*v._VECTptr)[j]));
}
gen tee = gen(te);
return tee;
}
gen CMfield::matevaln(gen M) {
vecteur te;
int m = (*M._VECTptr).size();
for (int j=0; j<m; j++) {
te.push_back(vectevaln((*M._VECTptr)[j]));
}
gen tee = gen(te);
return tee;
}
bool CMfield::checkSign(gen g, int &res){
// cout << "Check sign" << endl;
int savend = nd;
bool done = false;
while (!done && nd<ndmax) {
done = checkSignAttempt(g,res);
if (!done) {
increasePrecision();
}
}
nd = savend;
lowerPrecision();
return done;
}
gen CMfield::transpose(gen M) {
vecteur coo;
for (int j=0; j<3; j++) {
for (int k=0; k<3; k++) {
coo.push_back(M[k][j]);
}
}
return createMatrix(coo);
}
bool CMfield::checkSignAttempt(gen g, int &sign) {
// We will assume g is a polynomial in s, and it is reduced...
// (this may not be a good idea...)
if (operator_equal(g,zeronumber,&ct)) {
sign = 0;
return true;
}
else {
gen va = convert_interval(_subst(makesequence(g,s__IDNT_e,svi),&ct),nd,&ct);
if (va.type==8) {
if (is_greater(-prec,va,&ct)) {
sign = -1;
return true;
}
else {
if (is_greater(va,prec,&ct)) {
sign = 1;
return true;
}
else {
if (va==zeronumber) {
sign = 0;
return true;
}
else {
return false;
}
}
}
}
else {
if (is_greater(zeronumber,_right(va,&ct),&ct)) {
sign = -1;
return true;
}
else {
if (is_greater(_left(va,&ct),zeronumber,&ct)) {
sign = 1;
return true;
}
else {
return false;
}
}
}
}
}
bool CMfield::checkSignPol(gen par, gen mpol, gen q, int &res) {
if (operator_equal(zeronumber,q,&ct)) {
res = 0;
return true;
}
int savend = nd;
bool worked = false;
while (!worked && nd<ndmax) {
worked = checkSignPolAttempt(par,mpol,q,res);
if (!worked) {
increasePrecision();
}
}
nd = savend;
lowerPrecision();
return worked;
}
bool CMfield::checkSignPolAttempt(gen par, gen mpol, gen q, int &res) {
// u is an expression in s = gen of totally real subfield
// par is an interval, containing a single root of mpol
// we want to evaluate the sign of q at that root
// cout << "Precision: " << prec << endl;
gen sols = _realroot(makesequence(mpol,prec),&ct);
int nbsols = (*sols._VECTptr).size();
int count = 0;
int ksave;
for (int k=0; k<nbsols; k++) {
// beware sols[k][0] may not be an interval...
gen J = _intersect(makesequence(par,sols[k][0]),&ct);
if (J.type==3) {
count++;
ksave = k;
// cout << "ksave = " << ksave << endl;
}
}
if (count!=1) {
cerr << "Interval is not isolating, please double check!" << endl;
exit(1);
}
else {
gen te = _horner(makesequence(q,sols[ksave][0],s__IDNT_e),&ct);
if (is_greater(zeronumber,_right(te,&ct),&ct)) {
res = -1;
return true;
}
else {
if (is_greater(_left(te,&ct),zeronumber,&ct)) {
res = 1;
return true;
}
else {
// need to check if it is zero!
gen e = _egcd(makesequence(mpol,q,s__IDNT_e),&ct);
gen cr = _realroot(makesequence(e[2],prec),&ct);
for (int k=0; k<(*cr._VECTptr).size(); k++) {
gen J = _intersect(makesequence(sols[ksave][0],cr[k][0]),&ct);
if (J.type==3) {
res = 0;
return true;
}
}
// if we get here, then the expression is non-zero but we don't know its sign...
return false;
}
}
}
}
gen CMfield::matinverse(gen M) {
gen de = det(M);
red(de);
if (operator_equal(de,zeronumber,&ct)) {
cerr << "Matrix is not invertible" << endl;
exit(1);
}
else {
vecteur co;
co.push_back(M[1][1]*M[2][2] - M[1][2]*M[2][1]);
co.push_back(- M[0][1]*M[2][2] + M[0][2]*M[2][1]);
co.push_back(M[0][1]*M[1][2] - M[0][2]*M[1][1]);
co.push_back(- M[1][0]*M[2][2] + M[1][2]*M[2][0]);
co.push_back(M[0][0]*M[2][2] - M[0][2]*M[2][0]);
co.push_back(- M[0][0]*M[1][2] + M[0][2]*M[1][0]);
co.push_back(M[1][0]*M[2][1] - M[1][1]*M[2][0]);
co.push_back(- M[0][0]*M[2][1] + M[0][1]*M[2][0]);
co.push_back(M[0][0]*M[1][1] - M[0][1]*M[1][0]);
gen A = createMatrix(co);
A = A*inverse(de);
matred(A);
return A;
}
}
gen CMfield::matinverse2d(gen M) {
gen de = det(M);
red(de);
if (operator_equal(de,zeronumber,&ct)) {
cerr << "Matrix is not invertible" << endl;
exit(1);
}
else {
vecteur co;
co.push_back(M[1][1]);
co.push_back(-M[0][1]);
co.push_back(-M[1][0]);
co.push_back(M[0][0]);
gen A = createMatrix(2,2,co);
A = A*inverse(de);
matred(A);
return A;
}
}
bool CMfield::isInList(gen a, vecteur li) {
bool res = false;
int k=0;
while (!res && k<li.size()){
res = operator_equal(a,li[k],&ct);
k++;
}
return res;
}
gen CMfield::create_identity(int d) {
vecteur lico;
for (int j=0; j<d; j++) {
for (int k=0; k<d; k++) {
if (k==j) {
lico.push_back(onenumber);
}
else {
lico.push_back(zeronumber);
}
}
}
vecteur te;
for (int j=0; j<d; j++) {
vecteur tej;
for (int k=0; k<d; k++) {
tej.push_back(lico[d*j+k]);
}
te.push_back(tej);
}
return gen(te);
}
map<string,gen> * CMfield::sigma_data() {
if (!sigma_values) {
sigma_values = new map<string,gen>;
// SPORADIC VALUES
(*sigma_values)["s1"] = gen("-1+i*sqrt(2)",&ct);
(*sigma_values)["s2"] = gen("-1+i*(sqrt(5)+1)/2",&ct);
(*sigma_values)["s3"] = gen("-1+i*(sqrt(5)-1)/2",&ct);
(*sigma_values)["s4"] = gen("(-1+i*sqrt(7))/2",&ct);
(*sigma_values)["s5"] = gen("(sqrt(5)+i*sqrt(3))/2",&ct);
// (*sigma_values)["s5"] = gen("((1+i*sqrt(3))/2+(sqrt(5)-1)/2)/rootof(cyclotomic(18))",&ct);
(*sigma_values)["s6"] = gen("((1+i*sqrt(3))/2+(sqrt(5)+1)/2)/rootof(cyclotomic(18))",&ct);
(*sigma_values)["s7"] = gen("((1+i*sqrt(3))/2+2*cos(2*pi/7))/rootof(cyclotomic(18))",&ct);
(*sigma_values)["s8"] = gen("((1+i*sqrt(3))/2+2*cos(4*pi/7))/rootof(cyclotomic(18))",&ct);
(*sigma_values)["s9"] = gen("((1+i*sqrt(3))/2+2*cos(6*pi/7))/rootof(cyclotomic(18))",&ct);
(*sigma_values)["s10"] = gen("(1+sqrt(5))/2",&ct);
(*sigma_values)["s11"] = gen("(1-sqrt(5))/2",&ct);
(*sigma_values)["s1c"] = gen("-1-i*sqrt(2)",&ct);
(*sigma_values)["s2c"] = gen("-1-i*(sqrt(5)+1)/2",&ct);
(*sigma_values)["s3c"] = gen("-1-i*(sqrt(5)-1)/2",&ct);
(*sigma_values)["s4c"] = gen("(-1-i*sqrt(7))/2",&ct);
(*sigma_values)["s5c"] = gen("(sqrt(5)-i*sqrt(3))/2",&ct);
// (*sigma_values)["s5c"] = gen("rootof(cyclotomic(18))*((1-i*sqrt(3))/2+(sqrt(5)-1)/2)",&ct);
(*sigma_values)["s6c"] = gen("rootof(cyclotomic(18))*((1-i*sqrt(3))/2+(sqrt(5)+1)/2)",&ct);
(*sigma_values)["s7c"] = gen("rootof(cyclotomic(18))*((1-i*sqrt(3))/2+2*cos(2*pi/7))",&ct);
(*sigma_values)["s8c"] = gen("rootof(cyclotomic(18))*((1-i*sqrt(3))/2+2*cos(4*pi/7))",&ct);
(*sigma_values)["s9c"] = gen("rootof(cyclotomic(18))*((1-i*sqrt(3))/2+2*cos(6*pi/7))",&ct);
}
return sigma_values;
}
map<string,gen> * CMfield::T_data() {
if (!T_values) {
T_values = new map<string,gen>;
// SPORADIC VALUES
(*T_values)["S1"] = gen(makesequence(gen("(1+i*sqrt(7))/2",&ct),gen("1",&ct),gen("1",&ct)));
(*T_values)["S2"] = gen(makesequence(gen("1+(-1+i*sqrt(3))*(1+sqrt(5))/4",&ct),gen("1",&ct),gen("1",&ct)));
(*T_values)["S3"] = gen(makesequence(gen("0",&ct),gen("1",&ct),gen("1",&ct)));
(*T_values)["S4"] = gen(makesequence(gen("0",&ct),gen("1",&ct),gen("sqrt(2)",&ct)));
(*T_values)["S5"] = gen(makesequence(gen("0",&ct),gen("1",&ct),gen("(1+sqrt(5))/2",&ct)));
(*T_values)["E1"] = gen(makesequence(gen("i*sqrt(2)",&ct),gen("1",&ct),gen("1",&ct)));
(*T_values)["E2"] = gen(makesequence(gen("sqrt(2)",&ct),gen("(1+i*sqrt(3))/2",&ct),gen("sqrt(2)",&ct)));
(*T_values)["E3"] = gen(makesequence(gen("0",&ct),gen("1",&ct),gen("sqrt(3)",&ct)));
(*T_values)["H1"] = gen(makesequence(gen("(-1+i*sqrt(7))/2",&ct),gen("rootof([1,0,0,0,0,0],[1,1,1,1,1,1,1])",&ct),gen("rootof([1,0,0,0,0,0],[1,1,1,1,1,1,1])",&ct)));
(*T_values)["H2"] = gen(makesequence(gen("-1-rootof(cyclotomic(5))^4",&ct),gen("rootof(cyclotomic(5))^2",&ct),gen("rootof(cyclotomic(5))^2",&ct)));
(*T_values)["S1c"] = gen(makesequence(gen("(1-i*sqrt(7))/2",&ct),gen("1",&ct),gen("1",&ct)));
(*T_values)["S2c"] = gen(makesequence(gen("1+(-1-i*sqrt(3))*(1+sqrt(5))/4",&ct),gen("1",&ct),gen("1",&ct)));
(*T_values)["E1c"] = gen(makesequence(gen("-i*sqrt(2)",&ct),gen("1",&ct),gen("1",&ct)));
(*T_values)["E2c"] = gen(makesequence(gen("sqrt(2)",&ct),gen("(1-i*sqrt(3))/2",&ct),gen("sqrt(2)",&ct)));
(*T_values)["H1c"] = gen(makesequence(gen("(-1-i*sqrt(7))/2",&ct),gen("rootof([1,0,0],[1,1,1,1,1,1,1])",&ct),gen("rootof([1,0,0],[1,1,1,1,1,1,1])",&ct)));
(*T_values)["H2c"] = gen(makesequence(gen("-1-1/rootof(cyclotomic(5))^4",&ct),gen("1/rootof(cyclotomic(5))^2",&ct),gen("1/rootof(cyclotomic(5))^2",&ct)));
}
return T_values;
}
