default(parisizemax, 4000000000);
default(realprecision, 80);
default(recover, 0);

emitcoeffs(p, n, v) =
{
  for(k = 0, n,
    print1(polcoeff(p, k, v));
    if(k < n, print1(",")));
}

\\ Two different questions, kept apart on purpose.
\\
\\ Whether u is a square in the base decides whether adjoining its square root
\\ gives a new field at all. Testing that needs a variable ranked ABOVE the
\\ base's x, because nffactor wants the factored polynomial's main variable to
\\ outrank the field's; varhigher is the only way to get one.
\\
\\ Whether u generates the base decides whether charpoly(u, x^2) is the minimal
\\ polynomial of its square root. It is not the same question: charpoly is a
\\ power of the minimal polynomial whenever u sits in a proper subfield, so the
\\ substituted polynomial is reducible while the degree-24 field is perfectly
\\ good. u = -1 is the witness, and reading a reducible charpoly as "no field"
\\ silently discards exactly the imprimitive cases a kernel walk is hunting.
tv = varhigher("tv");
issquareinbase(nf0, M0, uu) = matsize(nffactor(nf0, tv^2 - Mod(lift(uu), M0)))[1] > 1;
generatesbase(uu) = polisirreducible(charpoly(uu, x));

{
M = Polrev([-4337,-14388,3142,34144,14133,-13068,-7174,1452,1044,-44,-56,0,1], x);
SPR = [17,41,4001];
n = 34;
G = vector(n);
G[1] = Polrev([-12047684744243/1010521154671,-19155536311258/1010521154671,33007530014598/1010521154671,970158434443/24646857431,-639111271345/59442420863,-16243652124447/1010521154671,86585457602/144360164953,2267371366854/1010521154671,49836642741/1010521154671,-123158160928/1010521154671,-2848479295/1010521154671,2282567358/1010521154671], x);
G[2] = Polrev([96515719785/144360164953,759103825690/144360164953,-249398958921/144360164953,-1875580523006/144360164953,221535487504/144360164953,893771924893/144360164953,-60682668408/144360164953,-153079870399/144360164953,5795728010/144360164953,10702072347/144360164953,-161207484/144360164953,-252417579/144360164953], x);
G[3] = Polrev([4097253336060/1010521154671,6428337088682/1010521154671,-16983362826770/1010521154671,-26666552348283/1010521154671,4995194169457/1010521154671,12011158500989/1010521154671,8963930242/144360164953,-102254502202/59442420863,-72749932554/1010521154671,2360297848/24646857431,2958967100/1010521154671,-1834190802/1010521154671], x);
G[4] = Polrev([-3021091832038/144360164953,-8613218001240/144360164953,6558467121606/144360164953,20533880665504/144360164953,-59691812402/8491774409,-8649158109568/144360164953,-377691350977/144360164953,1227591244300/144360164953,67439809479/144360164953,-67922537503/144360164953,-2339538370/144360164953,1285690719/144360164953], x);
G[5] = Polrev([1913757367295/1010521154671,6082920918230/1010521154671,-5435372097753/1010521154671,-14697465763651/1010521154671,1440320945061/1010521154671,6305856944901/1010521154671,9878638824/144360164953,-54714425900/59442420863,-29037735764/1010521154671,54531723532/1010521154671,1135463191/1010521154671,-1101099801/1010521154671], x);
G[6] = Polrev([43880801677618/1010521154671,83238393556365/1010521154671,-94806877318850/1010521154671,-202851916574322/1010521154671,24105374284427/1010521154671,86803966401352/1010521154671,255394050698/144360164953,-12411603455836/1010521154671,-573568136685/1010521154671,687834254605/1010521154671,21756433235/1010521154671,-13002838992/1010521154671], x);
G[7] = Polrev([99511842438240/1010521154671,290642024211977/1010521154671,-213293302278328/1010521154671,-693016559067411/1010521154671,24364204632428/1010521154671,295777211659516/1010521154671,2213719891050/144360164953,-42424475301949/1010521154671,-2539084124315/1010521154671,2369241282840/1010521154671,86437389469/1010521154671,-45242227010/1010521154671], x);
G[8] = Polrev([23613403870867/1010521154671,60159725795991/1010521154671,-90047496748035/1010521154671,-94831033276401/1010521154671,34189317193525/1010521154671,35927876709894/1010521154671,-463525969351/144360164953,-283236656288/59442420863,30929507809/1010521154671,251930010355/1010521154671,2829799290/1010521154671,-4486476280/1010521154671], x);
G[9] = Polrev([-6949511900173/1010521154671,-12448280055440/1010521154671,19376217049263/1010521154671,24598280919304/1010521154671,-453043281547/59442420863,-10548416285530/1010521154671,83846551078/144360164953,1503836848252/1010521154671,19478646381/1010521154671,-81991849052/1010521154671,-1592191473/1010521154671,1512304608/1010521154671], x);
G[10] = Polrev([21899464436435/144360164953,39355678010860/144360164953,-69392921881086/144360164953,-70249803390813/144360164953,25795653678533/144360164953,27808349209493/144360164953,-2285827488236/144360164953,-3810072366991/144360164953,-1154844563/144360164953,202762106062/144360164953,174294757/8491774409,-3671109714/144360164953], x);
G[11] = Polrev([-9313241592429/1010521154671,-26669404641088/1010521154671,24288747990999/1010521154671,63166045335544/1010521154671,-5277928444475/1010521154671,-26745003133018/1010521154671,-118100800418/144360164953,3815457456652/1010521154671,189928260507/1010521154671,-211399892092/1010521154671,-6911205657/1010521154671,3997988160/1010521154671], x);
G[12] = Polrev([-44638591385/8491774409,-2930974009876/144360164953,2663295506669/144360164953,7408099715549/144360164953,-92916062462/144360164953,-2913557226799/144360164953,-219080291543/144360164953,389143328893/144360164953,28506874382/144360164953,-1234011567/8491774409,-882736340/144360164953,395711827/144360164953], x);
G[13] = Polrev([22080795466718/1010521154671,42885274810785/1010521154671,-49919138129418/1010521154671,-102575597665551/1010521154671,11612387827732/1010521154671,43333236598132/1010521154671,155576488010/144360164953,-6150498945680/1010521154671,-299525257508/1010521154671,340071197728/1010521154671,652909108/59442420863,-6431601367/1010521154671], x);
G[14] = Polrev([4130255875904/1010521154671,18373734921943/1010521154671,-3814954801990/1010521154671,-37811835305327/1010521154671,-852291405345/1010521154671,15753459018104/1010521154671,106742102601/144360164953,-2190410999129/1010521154671,-98096567018/1010521154671,116942480129/1010521154671,3161408021/1010521154671,-2112707867/1010521154671], x);
G[15] = Polrev([-558544371388/144360164953,-12043346432/144360164953,2586082503815/144360164953,-1161802781909/144360164953,-1339133332788/144360164953,580173858222/144360164953,245170177752/144360164953,-90282497982/144360164953,-18221782801/144360164953,5481445097/144360164953,26480071/8491774409,-113956050/144360164953], x);
G[16] = Polrev([-4708600138551/1010521154671,-9840677086489/1010521154671,14835892772910/1010521154671,20979973064090/1010521154671,-5465023640534/1010521154671,-8505110875806/1010521154671,66148096914/144360164953,1166978009735/1010521154671,4385443598/1010521154671,-1505964490/24646857431,-805447864/1010521154671,1102048340/1010521154671], x);
G[17] = Polrev([96099604369170/1010521154671,209718796843960/1010521154671,-223207884487193/1010521154671,-518858241351927/1010521154671,49746770413671/1010521154671,12950374869207/59442420863,909129495572/144360164953,-31326646480518/1010521154671,-1551866820915/1010521154671,1733550669852/1010521154671,56789074501/1010521154671,-32781646320/1010521154671], x);
G[18] = Polrev([-122958447382291/1010521154671,-388336617489168/1010521154671,275370139541879/1010521154671,930110756203924/1010521154671,-679925005653/24646857431,-397535058055237/1010521154671,-3079266122405/144360164953,57063030045004/1010521154671,3455461305700/1010521154671,-3188417316092/1010521154671,-117036830314/1010521154671,1485621561/24646857431], x);
G[19] = Polrev([-279484502447/43935702377,-2939716092045/43935702377,-698416268420/43935702377,5656777728708/43935702377,1003039850636/43935702377,-2277262539668/43935702377,-1068632875/153086071,323658997456/43935702377,30828627898/43935702377,-18316218418/43935702377,-895591233/43935702377,357819906/43935702377], x);
G[20] = Polrev([-22207742446697/144360164953,-15973302625166/144360164953,50146921776044/144360164953,39126536521217/144360164953,-18042458804111/144360164953,-16269786593748/144360164953,1694058718434/144360164953,2258734455302/144360164953,-16494568675/144360164953,-120506100060/144360164953,-1453313389/144360164953,2177859576/144360164953], x);
G[21] = Polrev([25089798804264/1010521154671,71260535053854/1010521154671,-59007768100808/1010521154671,-167482660055057/1010521154671,675694987805/59442420863,72833885342528/1010521154671,421692445247/144360164953,-10535369292296/1010521154671,-583402316769/1010521154671,590038226726/1010521154671,20695987332/1010521154671,-11272787545/1010521154671], x);
G[22] = Polrev([300048407434/1010521154671,2642185183760/1010521154671,6720494080313/1010521154671,-2727520466155/1010521154671,-3099496177676/1010521154671,1324231748561/1010521154671,67048676147/144360164953,-287304163350/1010521154671,-41244721430/1010521154671,20300686079/1010521154671,1233822131/1010521154671,-448497727/1010521154671], x);
G[23] = Polrev([11741389/2279887,9618752/2279887,-30077495/2279887,-23957296/2279887,11825093/2279887,9692232/2279887,-1278470/2279887,-1309968/2279887,33015/2279887,67914/2279887,6/55607,-1188/2279887], x);
G[24] = Polrev([12233718456118/1010521154671,28560767168989/1010521154671,-31702283980950/1010521154671,-69137155763071/1010521154671,7342419555108/1010521154671,29401095058071/1010521154671,111202195368/144360164953,-4198172309298/1010521154671,-206186490180/1010521154671,232966280173/1010521154671,7634561913/1010521154671,-4415513646/1010521154671], x);
G[25] = Polrev([-2192386000799/1010521154671,-8121083447495/1010521154671,2456265470868/1010521154671,19283882208120/1010521154671,1211903670912/1010521154671,-8098293423744/1010521154671,-100973675748/144360164953,1155810304200/1010521154671,85224807063/1010521154671,-64704021520/1010521154671,-2659507092/1010521154671,1242841776/1010521154671], x);
G[26] = Polrev([-2192386000799/1010521154671,-7110562292824/1010521154671,2456265470868/1010521154671,19283882208120/1010521154671,1211903670912/1010521154671,-8098293423744/1010521154671,-100973675748/144360164953,1155810304200/1010521154671,85224807063/1010521154671,-64704021520/1010521154671,-2659507092/1010521154671,1242841776/1010521154671], x);
G[27] = Polrev([339256504707/59442420863,35052628610452/1010521154671,-1291517719769/1010521154671,-97289084318056/1010521154671,-9468751457338/1010521154671,42511472714736/1010521154671,597522657436/144360164953,-6173222493288/1010521154671,-483149309004/1010521154671,20530840090/59442420863,14891179420/1010521154671,-6751576128/1010521154671], x);
G[28] = Polrev([14154660411951/1010521154671,43579325220884/1010521154671,-27954194402439/1010521154671,-118526412673992/1010521154671,1013791434000/1010521154671,3006076873344/59442420863,435619773576/144360164953,-7334464586376/1010521154671,-453882634014/1010521154671,409227733454/1010521154671,15109250056/1010521154671,-7804697736/1010521154671], x);
G[29] = Polrev([13402500496858/1010521154671,18080756105015/1010521154671,-30267720454485/1010521154671,-30437147660170/1010521154671,9748277744910/1010521154671,11650416413772/1010521154671,-90483277139/144360164953,-1591707945936/1010521154671,-28332772680/1010521154671,86047104632/1010521154671,1961776357/1010521154671,-1599266757/1010521154671], x);
G[30] = Polrev([4771529445766/1010521154671,38357818330030/1010521154671,2469993784747/1010521154671,-107365862060844/1010521154671,-11779319844183/1010521154671,47570622224616/1010521154671,694336634658/144360164953,-6948213131532/1010521154671,-554274349380/1010521154671,394128379255/1010521154671,17017754946/1010521154671,-7641680640/1010521154671], x);
G[31] = Polrev([-505976521523/144360164953,-3265080563267/144360164953,2040951811489/144360164953,7857481315574/144360164953,-13160648909/8491774409,-3343265793033/144360164953,-174382001689/144360164953,477563068184/144360164953,28008831263/144360164953,-26557130796/144360164953,-943915835/144360164953,503892549/144360164953], x);
G[32] = Polrev([9241291105008/1010521154671,13710135987809/1010521154671,-20924502713922/1010521154671,-29296090316604/1010521154671,10010928520108/1010521154671,11346791449892/1010521154671,-260997301025/144360164953,-1390863657193/1010521154671,138228362095/1010521154671,58025750278/1010521154671,-3491387884/1010521154671,-661688374/1010521154671], x);
G[33] = Polrev([504884914724896/1010521154671,1400031649238519/1010521154671,-27429238432952/24646857431,-3350680876691444/1010521154671,183373120264016/1010521154671,1415522672008516/1010521154671,8767466417249/144360164953,-201720941786743/1010521154671,-11215921187243/1010521154671,11205183181110/1010521154671,391104610424/1010521154671,-212904774074/1010521154671], x);
G[34] = Polrev([-1,0,0,0,0,0,0,0,0,0,0,0], x);

if(version[1] * 100 + version[2] < 213,
  print("FAIL|needs PARI 2.13 or later for bnfunits and 3-argument bnfisunit");
  quit);
print("PARI|", version[1], ".", version[2], ".", version[3]);

nf = nfinit(M);
print("NF|", nf.sign[1], ",", nf.sign[2], "|", poldegree(M));

\\ The 34 generators above are shipped rather than recomputed, which is only
\\ honest if they generate the same S-unit group a fresh bnfunits would. Express
\\ each one in the freshly computed basis; the change of basis is unimodular
\\ exactly when the two sets generate the same group. Costs 50 milliseconds.
bnfM = bnfinit(M, 1);
SMp = List();
for(ip = 1, #SPR,
  pd = idealprimedec(bnfM, SPR[ip]);
  for(j = 1, #pd, listput(SMp, pd[j])));
SMp = Vec(SMp);
UM = bnfunits(bnfM, SMp);
print("BASE|classno=", bnfM.no, "|primes=", #SMp, "|width=", #UM[1]);
Emat = matrix(#UM[1], n);
for(j = 1, n,
  ex = bnfisunit(bnfM, Mod(G[j], M), UM);
  if(#ex == 0, print("FAIL|shipped generator ", j, " is not an S-unit"); quit);
  for(t = 1, #ex, Emat[t, j] = lift(ex[t])));
dt = if(matsize(Emat)[1] == matsize(Emat)[2], matdet(Emat), 0);
print("SAMEGROUP|", if(abs(dt) == 1, 1, 0), "|det=", dt);
if(abs(dt) != 1,
  print("FAIL|the shipped generators do not generate the S-unit group");
  quit);

\\ Sign block: one row per real place of the base, in the order nfeltsign
\\ returns them, 1 meaning negative. Signature 24 upstairs is the all-zero
\\ right-hand side; signature 24-2k asks for negative at exactly k places.
r1 = nf.sign[1];
SB = matrix(r1, n);
print("SIGNPLACES|", r1);
for(j = 1, n,
  sv = nfeltsign(nf, Mod(G[j], M));
  print1("SB|", j, "|");
  for(t = 1, r1,
    print1(if(sv[t] < 0, "1", "0"));
    SB[t, j] = if(sv[t] < 0, 1, 0));
  print());
print("SBRANK|", matrank(Mod(SB, 2)), "|", r1);

print("DECOY|Q");
for(j = 1, n,
  nu = nfeltnorm(nf, Mod(G[j], M));
  print1("QN|", j, "|");
  if(nu < 0, print1("1"), print1("0"));
  for(t = 1, #SPR, print1(valuation(nu, SPR[t]) % 2));
  print());

SF = nfsubfields(nf, 3);
print("CUBICS|", #SF);
if(#SF == 0, print("FAIL|no cubic"); quit);

for(i = 1, #SF,
  P = SF[i][1];
  em = SF[i][2];
  \\ Every class coordinate below is read in this exact model of the cubic, so a
  \\ different model would silently change them. Assert rather than adapt.
  if(P != Polrev([4337, 2685, -160, 1], x),
    print("FAIL|cubic subfield is not the published model, classes would not transport");
    quit);
  k = poldegree(P);
  print1("SF|", i, "|", k, "|"); emitcoeffs(P, k, x); print();
  \\ The relative norm depends on which embedding of F into the base is used, so
  \\ the target class is only meaningful once one is pinned. Say how many there
  \\ are and which one this run took: everything below is relative to it.
  allem = nfisincl(P, M);
  print("EMB|", i, "|", #allem, "|", em);
  Py = subst(P, x, y);
  bnfF = bnfinit(Py, 1);
  print("SFBNF|", i, "|", bnfF.no, "|", bnfF.sign[1], ",", bnfF.sign[2]);
  SFp = List();
  for(ip = 1, #SPR,
    pd = idealprimedec(bnfF, SPR[ip]);
    for(j = 1, #pd, listput(SFp, pd[j])));
  SFp = Vec(SFp);
  UF = bnfunits(bnfF, SFp);
  wF = #UF[1];
  print("SFW|", i, "|", wF, "|", #SFp);

  fa = nffactor(bnfF, M);
  rel = 0;
  for(j = 1, matsize(fa)[1],
    ff = fa[j, 1];
    if(poldegree(ff, x) != 4, next);
    tt = Mod(0, M);
    for(d = 0, poldegree(ff, x),
      cd = lift(polcoeff(ff, d, x));
      tt = tt + Mod(subst(cd, y, em), M) * Mod(x, M)^d);
    if(tt == 0, rel = ff; break));
  if(rel == 0, print("SFFAIL|", i, "|no relative polynomial"); next);

  print1("REL|", i, "|"); print(lift(rel));
  CB = matrix(wF, n);
  ok = 1;
  for(j = 1, n,
    nu = polresultant(rel, G[j], x);
    ex = bnfisunit(bnfF, nu, UF);
    if(#ex == 0, print("SFFAIL|", i, "|norm of generator ", j, " is not an S-unit"); ok = 0; break);
    print1("NB|", i, "|", j, "|");
    for(t = 1, #ex,
      print1(lift(ex[t]) % 2);
      CB[t, j] = lift(ex[t]) % 2);
    print());
  if(ok == 0, next);

  \\ What the two blocks can and cannot reach together on this base.
  \\ Both blocks are one row per coordinate and one column per generator, so
  \\ the stacked system is (r1 + wF) by n and the right-hand side is a vector
  \\ of r1 wanted signs on top of wF wanted class bits. The NB and SB lines
  \\ above print the transpose, one line per generator, because that is what
  \\ reads on a screen.
  A = matconcat([SB; CB]);
  jr = matrank(Mod(A, 2));
  print("CBRANK|", i, "|", matrank(Mod(CB, 2)), "|", wF);
  print("JOINTRANK|", i, "|", jr, "|", r1 + wF);
  print("KERNEL|", i, "|", n - jr);

  \\ Now aim. Signature 24 is every real place positive, so the top of the
  \\ right-hand side is zero; CTARGET is the class we want rather than the
  \\ one a default build keeps writing.
  \\ The target is named by an element of F rather than by a bit string. Those
  \\ nine coordinates live in whichever generating set bnfunits returned, and a
  \\ generator for a prime is pinned only up to units and sign, so another build
  \\ can shear the basis and make the same string name a different class. The
  \\ element does not move. Its coordinates here read 010000100.
  WREP = Mod(Polrev([10031/8371, 3653/8371, -25/8371], y), Py);
  exw = bnfisunit(bnfF, WREP, UF);
  if(#exw == 0, print("AIMFAIL|", i, "|target element is not an S-unit of F at this S"); next);
  CTARGET = vectorv(wF, t, lift(exw[t]) % 2);
  print1("CTARGET|", i, "|"); for(t = 1, wF, print1(CTARGET[t])); print();
  rhs = concat(vectorv(r1, t, 0), CTARGET);
  sol = matsolvemod(A, 2, rhs, 1);
  if(sol == 0, print("UNREACHABLE|", i); next);
  e = lift(sol[1]);
  u = Mod(1, M);
  for(j = 1, n, if(e[j] % 2 == 1, u = u * Mod(G[j], M)));
  print("AIMED|", i, "|weight=", sum(j = 1, n, e[j] % 2));

  \\ Read the two coordinates back off u itself rather than trusting the solve,
  \\ and then actually compare them. Printing them and moving on is not a check:
  \\ it leaves the reader to notice a mismatch nobody looked for.
  sv = nfeltsign(nf, u);
  usign = 0; print1("USIGN|");
  for(t = 1, r1, b0 = if(sv[t] < 0, 1, 0); print1(b0); if(b0 != 0, usign = 1));
  print();
  exu = bnfisunit(bnfF, polresultant(rel, lift(u), x), UF);
  uclassok = 1; print1("UCLASS|");
  for(t = 1, wF, b0 = lift(exu[t]) % 2; print1(b0); if(b0 != CTARGET[t], uclassok = 0));
  print();
  if(usign != 0, print("VERIFYFAIL|", i, "|solved u is not totally positive"); quit);
  if(uclassok == 0, print("VERIFYFAIL|", i, "|solved u is not in the target class"); quit);
  print("VERIFIED|", i, "|sign and class read back off u match the right-hand side");

  \\ The degree-24 field is the base with a square root of u adjoined. Whether
  \\ that field is new, and whether this particular construction writes it, are
  \\ the two separate tests defined at the top.
  if(issquareinbase(nf, M, u),
    print("AIMFAIL|", i, "|u is a square in the base, so adjoining its root gives no new field");
    next);
  if(!generatesbase(u),
    print("AIMFAIL|", i, "|u lies in a proper subfield, so charpoly(u,x^2) is not the minimal polynomial of its square root; the degree-24 field exists but this construction does not write it, and rnfequation is the route that does");
    next);
  Q = subst(charpoly(u, x), x, x^2);
  if(!polisirreducible(Q),
    print("FAIL|", i, "|charpoly route gave a reducible polynomial for a base generator");
    quit);
  R = polredbest(Q);
  print("FIELD|", i, "|", poldegree(R), "|realroots=", polsturm(R));
  print1("POLY|", i, "|"); emitcoeffs(R, poldegree(R), x); print();

  \\ The solve has a kernel, so this class and this signature name 2^KERNEL
  \\ fields rather than one, and they do not all carry the same group. Walking
  \\ the kernel basis emits eight more of them. Label these and see for
  \\ yourself: the class narrows the label, it does not fix it.
  K = matker(Mod(A, 2));
  nw = min(8, matsize(K)[2]);
  for(b = 1, nw,
    ew = vectorv(n, j, (e[j] + lift(K[j, b])) % 2);
    uw = Mod(1, M);
    for(j = 1, n, if(ew[j] == 1, uw = uw * Mod(G[j], M)));
    if(issquareinbase(nf, M, uw), print("WALKSKIP|", i, "|", b, "|square"); next);
    if(!generatesbase(uw), print("WALKSKIP|", i, "|", b, "|propersubfield"); next);
    Qw = subst(charpoly(uw, x), x, x^2);
    Rw = polredbest(Qw);
    exw = bnfisunit(bnfF, polresultant(rel, lift(uw), x), UF);
    print1("WALK|", i, "|", b, "|");
    for(t = 1, wF, print1(lift(exw[t]) % 2));
    print1("|", polsturm(Rw), "|");
    emitcoeffs(Rw, poldegree(Rw), x); print());
);
print("PREPARED|", n);
}
quit;