// Verifies trinomial log file in extended format. // Degree r and input file-name is hardcoded, see below. // RPB, 20040124..20040206 (minor revisions 20070318) r := 132049; // r is hardcoded log := Read("i132049.log-ext"); // as is file name (how to avoid this - // perhaps read filename as a string ? // Following limits are so the program does not take too long or run // out of memory on inputs that are too large for Magma to handle. // They can probably be increased if you have enough memory and/or a // recent version of Magma. // dlim := r; // highest degree to check divisibility by factor ilim := 4000; // highest degree to check irreducibility of factor slim := 100; // highest degree to use IsIrreducible (see below) ulim := 24000; // highest degree to check irreducibility of trinomial GF2 := GF(2,1); P := PolynomialRing(GF2); nextstring := function(s, sind) ss := ""; done := (sind ge #s); s0 := StringToCode("0"); s9 := StringToCode("9"); sa := StringToCode("a"); sz := StringToCode("z"); j := sind; while (not done) do j := j + 1; sv := StringToCode(s[j]); if (sv in [s0..s9]) or (sv in [sa..sz]) then ss := ss * s[j]; else done := true; while (not ((sv in [s0..s9]) or (sv in [sa..sz])) and (j lt #s)) do j := j + 1; sv := StringToCode(s[j]); end while; end if; end while; if (j lt #s) then j := j - 1; end if; return ss, j; end function; MyIsIrreducible := function(Q) // // Tests irreducibility of polynomial Q, // used because built-in IsIrreducible uses too much memory // and may bomb out, at least with Magma V2.6-2 // // This is slower but uses less memory and more reliable // degq := Degree(Q); if (degq le slim) then // Use built in for small degrees return IsIrreducible(Q); end if; if (degq le 1) then return true; end if; Qquo := quo; // work mod Q xq := Qquo!x; if (xq^(2^degq) ne xq) then return false; end if; if IsPrime(degq) then return true; end if; facs := Factorization(degq); for j in [1..#facs] do if (xq^(2^(degq div facs[j][1])) eq xq) then return false; end if; end for; return true; end function; // Convert hex num to polynomial // d is a (strict) upper bound on the degree of the output polynomial HexToPolyFast := function (nf, d) if d le 1 then return P!nf; else k := d div 2; q, r := Quotrem(nf,2^k); h := $$ (q, d - k); l := $$ (r, k); return h * x^k + l; end if; end function; if (not IsPrime(r)) then print "Warning: r = ", r, " is not prime"; end if; dmax := 0; dsum := 0; ktver := 0; ktnot := 0; cutoff := 0; sind := 0; while (sind lt #log) do ss, sind := nextstring(log, sind); s := StringToInteger(ss); ss, sind := nextstring(log, sind); if ((ss eq "u") or (ss eq "primitive") or (ss eq "irreducible")) then if (r le ulim) then T := x^r + x^s + 1; Pquo := quo; xq := Pquo!x; if (MyIsIrreducible(T)) then print r, s, "u"; ktver := ktver + 1; else print r, s, "not irreducible"; quit; end if; else print r, s, "u irreducible, not checked"; // checking these is too slow ktnot := ktnot + 1; end if; else degf := StringToInteger(ss); // Look ahead ss, sindtemp := nextstring(log, sind); if (#ss lt 1) then cutoff := degf + 1; else if ((ss[1] eq "p") or (ss[1] eq "q")) then cutoff := degf; else cutoff := degf + 1; end if; end if; if (degf ge cutoff) then hexnum, sind := nextstring(log, sind); else hexnum := ""; T := x^r + x^s + 1; U := x^(2^degf) - x; G := Gcd(T,U); if (Degree(G) le 1) then print r, s, degf; print "No factor of expected degree (or degf = 1)"; quit; end if; FG := Factorization(G); if (Degree(FG[1][1]) eq degf) then ktver := ktver + 1; dsum := dsum + degf; print r, s, degf; else print r, s, degf; print "wrong degree factor ", FG[1][1]; quit; end if; end if; if (hexnum ne "") then ssv := s; nf := 0; base0 := StringToCode("0"); base9 := StringToCode("9"); basea := StringToCode("a"); if (hexnum[1] ne "p") and (hexnum[1] ne "q") then print r, s, degf, hexnum; print "First character must be p or q:", hexnum; quit; end if; if (hexnum[1] eq "q") then s := r - s; // Consider reverse end if; for j := 2 to #hexnum do dig := StringToCode(hexnum[j]); if (dig le base9) then // convert hex dig := dig - base0; else dig := dig + 10 - basea; end if; if (dig lt 0) or (dig gt 15) then print r, s, degf, hexnum; print "Illegal hex digit:", hexnum[j]; quit; end if; nf := 16*nf + dig; end for; nfs := nf; F := HexToPolyFast (nf, degf + 1); if (Degree(F) ne degf) or ((F mod x) ne 1) then print r, s, degf, hexnum; print "Warning: degrees do not match or constant term zero."; print degf, Degree(F); print F; // In rare cases F may have multiple factors of degree degf, and // we could choose one of these for the log (since March 2007 we // choose the lexicographically least factor) if (not MyIsIrreducible(F)) then print "F is reducible, factors follow"; print Factorisation(F); end if; quit; end if; if (degf gt dlim) then print r, s, degf, "not checked"; // because too slow, degf large ktnot := ktnot + 1; else if (degf le ilim) then // skip next if degf large if (not MyIsIrreducible(F)) then // give warning if F reducible print r, s, degf, hexnum; print "F is reducible, factors follow"; print Factorisation(F); quit; end if; end if; Pquo := quo; // To save time and space, work mod F xq := Pquo!x; if (r le s) then print "r =", r, "should be greater than s =", s; quit; end if; Tq := xq^r + xq^s + 1; // Compute T mod F if (Tq eq 0) then ktver := ktver + 1; dsum := dsum + Degree(F); print r, ssv, Degree(F), hexnum; if (dmax lt Degree(F)) then dmax := Degree(F); end if; else print r, ssv, Degree(F), hexnum, "NOT VERIFIED"; quit; end if; end if; end if; end if; end while; print "Verified", ktver, "and skipped ", ktnot; if (dmax gt 0) then print "Max degree verified in case p/q is", dmax; end if; if (ktver gt 0) then print "Mean degree verified is", RealField(6)!(dsum/(ktver + 0.0)); end if; quit;