"Fossies" - the Fresh Open Source Software Archive 
"Fossies" - the Fresh Open Source Software Archive A hint: This file contains one or more very long lines, so maybe it is better readable using the pure text view mode that shows the contents as wrapped lines within the browser window.
1 #_f2:= X^4 - +3*X^3-X^2 + 7; 2 #_k:= pAdicField(2,30); 3 #_g:= X^2 + 2*X + 2; 4 #_ek:= TotallyRamifiedExtension(_k,_g); 5 #_ekx:= PolynomialAlgebra(_ek); 6 # _f2:= Coerce(_ekx,_f2); 7 dd_realPrec:=20; 8 9 10 11 12 #_Z2q:= Quotient(pAdicRing(2,30),2^30); 13 #_Z2qx:= PolynomialAlgebra(_Z2q); 14 #_x:= _Z2qx.1; 15 #_f:= _x^6+6*_x^4+2^5*_x^3+12*_x^2-3*2^6*_x+33*2^3; 16 17 #_Z3:= pAdicRing(3,30); 18 #_Z3q:= Quotient(_Z3,3^30); 19 #_Z3qx:= PolynomialAlgebra(_Z3q); 20 #_y:= _Z3qx.1; 21 #_g:= _y^8+4*_y^6+6*_y^4+7*_y^2+9*_y+13; 22 23 24 #//freeze; 25 26 #///////////////////////////////////////////////////////////////////////////// 27 #// Automorphisms are in class_field.m now 28 #///////////////////////////////////////////////////////////////////////////// 29 #/* 30 #declare attributes Map:LocalAutoGenerators; 31 #declare attributes Map:LocalRingMap; 32 33 #function eval_automorphism(x, gens, pos) 34 # if pos eq #gens then 35 # return x; 36 # end if; 37 38 # R := Parent(x); 39 # P := PolynomialAlgebra(R); 40 # return Evaluate( 41 # P ! [eval_automorphism(i, gens, pos + 1) : i in Eltseq(x)], gens[pos] 42 # ); 43 #end function; 44 45 #intrinsic Automorphisms(L::RngPad) -> SeqEnum 46 #{Return the automorphisms of the local ring L. Automorphisms are returned as a sequence of maps} 47 # if L eq PrimeRing(L) then 48 # curr_gens := [L ! 1]; 49 # curr := map<L -> L | x :-> eval_automorphism(x, curr_gens, 1)>; 50 # curr`LocalAutoGenerators := curr_gens; 51 # return [curr]; 52 # end if; 53 54 # R := BaseRing(L); 55 # maps := Automorphisms(R); 56 # f := DefiningPolynomial(L); 57 58 # res := []; 59 # if IsEisenstein(f) then 60 # P := PolynomialAlgebra(L); 61 # for m in maps do 62 # gens := ChangeUniverse(m`LocalAutoGenerators, L); 63 # for r in Roots(P ! [m(x) : x in Eltseq(f)]) do 64 # curr_gens := Insert(gens, 1, r[1]); 65 # curr := map<L -> L | x :-> eval_automorphism(x, curr_gens, 1)>; 66 # curr`LocalAutoGenerators := curr_gens; 67 # Append(~res, curr); 68 # end for; 69 # end for; 70 # elif IsInertial(f) then 71 # for m in maps do 72 # gens := ChangeUniverse(m`LocalAutoGenerators, L); 73 # r := L.1; 74 # for i in [1..Degree(L)] do 75 # curr_gens := Insert(gens, 1, r); 76 # curr := map<L -> L | x :-> eval_automorphism(x, curr_gens, 1)>; 77 # curr`LocalAutoGenerators := curr_gens; 78 # Append(~res, curr); 79 # r := GaloisImage(r, 1); 80 # end for; 81 # end for; 82 # else 83 # error "Error: Unsupported defining polynomial"; 84 # end if; 85 86 # return res; 87 #end intrinsic; 88 89 #function eval_field_automorphism(x, m) 90 # F := Parent(x); 91 # R := RingOfIntegers(F); 92 # pi := UniformizingElement(F); 93 94 # xv := Valuation(x); 95 # xu := ShiftValuation(x, -xv); 96 # return (F ! m(UniformizingElement(R)))^xv * F ! m(R ! xu); 97 #end function; 98 99 #function construct_field_map(L, m) 100 # res := map<L -> L | x :-> eval_field_automorphism(x, m)>; 101 # res`LocalRingMap := m; 102 # return res; 103 #end function; 104 105 #intrinsic Automorphisms(L::FldPad) -> SeqEnum 106 #{Return the automorphisms of the local ring L. Automorphisms are returned as a sequence of maps} 107 # maps := Automorphisms(RingOfIntegers(L)); 108 # return [construct_field_map(L, m) : m in maps]; 109 #end intrinsic; 110 111 #function auto_eq(m1, m2) 112 # return m1`LocalAutoGenerators eq m2`LocalAutoGenerators; 113 #end function; 114 115 #function auto_mult(m1, m2) 116 # gens1 := m1`LocalAutoGenerators; 117 # gens2 := m2`LocalAutoGenerators; 118 # gens := []; 119 # for i in [1..#gens1] do 120 # Append(~gens, eval_automorphism(gens1[i], gens2, 1)); 121 # end for; 122 123 # L := Universe(gens1); 124 # m := map<L -> L | x :-> eval_automorphism(x, gens, 1)>; 125 # m`LocalAutoGenerators := gens; 126 # return m; 127 #end function; 128 129 #intrinsic AutomorphismGroup(L::RngPad) -> GrpPerm 130 #{Return the automorphism group of L as a permutation group} 131 # maps := Automorphisms(L); 132 133 # G, m1 := GenericGroup( 134 # maps : Mult := auto_mult, Eq := auto_eq, Id := maps[1] 135 # ); 136 # m2, G2 := CosetAction(G, sub<G|>); 137 # m := Inverse(m2) * m1; 138 # return 139 # G2, 140 # map< 141 # G2 -> PowerStructure(Map) | g :-> m(g), h :-> h @@ m 142 # >; 143 #end intrinsic; 144 145 #intrinsic AutomorphismGroup(L::FldPad) -> GrpPerm 146 #{Return the automorphism group of L as a permutation group} 147 # G, m := AutomorphismGroup(RingOfIntegers(L)); 148 149 # return 150 # G, 151 # map< 152 # G -> PowerStructure(Map) | 153 # g :-> construct_field_map(L, m(g)), 154 # h :-> h`LocalRingMap @@ m 155 # >; 156 #end intrinsic; 157 # 158 #/////////////////////////////////////////////////////////////////////// 159 #// general functions 160 #////////////////////////////////////////////////////////////////////// 161 162 #intrinsic BaseRing(L::FldPad) -> FldPad 163 #{The base field of L} 164 # return BaseField(L); 165 #end intrinsic; 166 167 168 #//--------------------------------------------------------------------- 169 #// FACTORIZATION 170 #//--------------------------------------------------------------------- 171 172 # default precision 173 _realPrec := 10; 174 175 #/////////////////////////////////////////////////////////////////////// 176 #// some functions for polynomials that probably should be somewhere 177 #// else 178 179 _ChangePrecisionCoeff:= function(PR,n) 180 local Co,Co2; 181 Co:= CoefficientRing(PR); 182 Co2:= ChangePrecision(Co,n); 183 return PolynomialAlgebra(Co2); 184 end; 185 186 187 188 189 IsMonic:= function(f) 190 return LeadingCoefficient(f)= One(CoefficientRing(Parent(f))); 191 end; 192 193 #Computes a polynomial g such that g mod Lm[i] = La[i], Lm and La are lists of Polynomials 194 ChineseRemainder:= function(Lm,La) 195 local PR,Mod,CR,i,gcd,e2,e1; 196 197 PR := Parent(Lm[1]); 198 Mod := Lm[1]; 199 CR := La[1]; 200 for i in [2..Length(Lm)] do 201 gcd := XGCD(Mod,Lm[i]); 202 203 if Base(gcd)<> One(PR) then 204 return Error( "Error, ChineseRemainderTheorem: modules must be coprime."); 205 fi; 206 207 e2 := Mod*(gcd.ext1); 208 e1 := Lm[i]*(gcd.ext2); 209 210 Mod := Mod*Lm[i]; 211 CR := (e1*CR+e2*La[i]) mod Mod; 212 od; 213 return CR; 214 end; 215 216 PGcd:=function(f,g) 217 local PR,R,pi,prec,S,H,L,h,loss,l,i,H,p,k; 218 #Print("PGcd g: ",g,"\nf ",f,"\n"); 219 #// only works for quotient rings. 220 if Degree(f) < 0 or Degree(g) < 0 then return [One(Parent(f)),0]; fi; 221 # ////////////////assert IsMonic(f) or IsMonic(g); 222 #Print("f: ",f); 223 #Print("\ng: ",g,"\n"); 224 PR := Parent(f); 225 R:= CoefficientRing(PR); 226 pi:= UniformizingElement(R); 227 prec := Precision(R); 228 S := SylvesterMatrix(f,g); 229 H := EchelonForm(S); 230 #Print("\nPGcd HNF",H,"\n"); 231 H:=H.base; 232 L := ElementToSequence(H[Rank(H)]); 233 #Print("L ",L); 234 h := Coerce(PR,Reverse(L)); 235 l := ElementToSequence(h); 236 k:= []; 237 for i in [1..Length(l)] do 238 k[i]:= Valuation(l[i]); 239 od; 240 loss:=Minimum(k); 241 H:= ElementToSequence(h); 242 p:= pi^loss; 243 for i in [1..Length(H)] do 244 H[i]:= Div(H[i],p); 245 od; 246 h := Coerce(PR,H); 247 #Print("\nPGcd ",h," loss ",loss,"\n"); 248 return [h, loss]; 249 end; 250 251 252 #/////////////////////////////////////////////////////// 253 254 is_squarefree:= function(f) 255 local val,df,L,gcd,loss; 256 if Degree(f) <= 0 then 257 return [true, f]; 258 fi; 259 260 # vprint RoundFour,4:" -- disc"; 261 # vtime RoundFour,4: 262 if Is(Type(f),elt-alg^pol/res^pad) or Is(Type(f),elt-alg^pol/ord^pad) then 263 val := Valuation(Discriminant(f)); 264 if val < Precision(CoefficientRing(Parent(f))) then 265 #vprint RoundFour,3:" -- polynomial is squarefree by discriminant"; 266 return [true, f, 0]; 267 fi; 268 fi; 269 df := Derivative(f); 270 #vprint RoundFour,4: " -- gcd"; 271 #vtime RoundFour,4: 272 L := PGcd(f,df); 273 gcd:= L[1]; 274 loss:= L[2]; 275 if Degree(gcd) < 1 then 276 #vprint RoundFour,3:" -- polynomial is squarefree by gcd"; 277 return [true, f, loss]; 278 else 279 return [false,Div(f,gcd), loss]; 280 fi; 281 282 end; 283 284 #//TODO ??? 285 #//intrinsic IsSquareFree(f::RngUPolElt[RngPadResExtElt]) 286 #//{return 287 #// is_squarefree(f) 288 289 290 291 #/////////////////////////////////////////////////////// 292 #// functions for handling polynomials with denominators 293 294 #///////////////////////////////////////////////////////////// 295 #// in the round four elements are represented by polynomials 296 #// and the exponent of pi in the denomiantor 297 #// we also store additional information in the structure 298 #theta has to be a record with .elt and .denval components 299 300 # 'cancels out' powers of the uniformizer pi in the polynomial 301 # and its denominator. 302 303 elt_reduce := function(theta) 304 local seq,vallist,i,minval,pi; 305 if theta.elt = Coerce(Parent(theta.elt),0) then 306 theta.denval:= 0; 307 return theta; 308 elif theta.elt<>"null" then 309 seq := ElementToSequence(theta.elt); 310 vallist:= []; 311 for i in [1..Length(seq)] do 312 vallist[i]:= Valuation(seq[i]); 313 od; 314 minval := Minimum(vallist); 315 if minval > 0 and theta.denval>= minval then 316 pi := UniformizingElement(CoefficientRing(Parent(theta.elt))); 317 theta.elt := Div(theta.elt,pi^minval); 318 theta.denval := theta.denval - minval; 319 fi; 320 fi; 321 return theta; 322 end; 323 324 # constructor 325 eltRec := function(phi, denval) 326 local theta; 327 # ELT := recformat 328 # <elt:RngUPolElt, // the numerator of the element 329 # denval:IntegerRing(), // the denominator of elt is pi^den, 330 # chi:RngUPolElt, // its characteristic polynomial 331 # passes_hensel_test, // true if it passes the hensel test 332 # nu:RngUPolElt, // irreducible factor of chi mod pi 333 # f:IntegerRing(), // deg(nu) 334 # res_fact, // factorization of chi mod pi 335 # passes_newton_test, // true if it passes the newton test 336 # v_star, // v_star valuations of elt 337 # e:IntegerRing() // numerator of vstar 338 # >; 339 theta := rec(elt := phi, denval:= denval, chi:="null", passes_hensel_test:="null", nu:="null",f:="null", res_fact:="null",passes_newton_test:="null",v_star:="null",e:="null"); 340 theta :=elt_reduce(theta); 341 return theta; 342 343 end; 344 345 #////////////////////////////////////////////////// 346 #// Find the inverse of some (invertible) element theta modulo Phi. 347 #// the inverse is given as res/pi^denval 348 elt_inverse := function(theta, Phi) 349 350 local pi,P,R,M,L,E,T,res2,gcd,denval,uni; 351 pi := UniformizingElement(CoefficientRing(Parent(Phi))); 352 P := Parent(theta.elt); 353 R := CoefficientRing(P); 354 M := SylvesterMatrix(theta.elt, Phi); 355 L:=EchelonForm(M); 356 E:= L.base; 357 T:= L.ext1; 358 res2 := (Coerce(P,Reverse(ElementToSequence(ExtractBlock(T, NumberOfRows(T), 1,1,Degree(Phi)))))) mod Phi; 359 gcd := Coerce(CoefficientRing(E),E[NumberOfRows(E)][NumberOfColumns(E)]); 360 denval := Valuation(gcd); 361 uni := Div(gcd,pi^denval); 362 res2 := Div(res2,uni); 363 denval := denval - theta.denval; 364 if denval < 0 then 365 res2 := res2*pi^(-denval); 366 denval := 0; 367 fi; 368 return eltRec(res2,denval); 369 end; 370 371 #////////////////////////////////////////////////// 372 # compute the v_star valuation from a given characteristic polynomial f 373 v_star := function(f) 374 local coeffs,L,i, res2; 375 if IsMonic(f)= false then 376 Print("f not monic\n"); 377 return; 378 fi; 379 coeffs := Reverse(ElementToSequence(f)); 380 Butfirst_(coeffs); 381 L:=[]; 382 for i in [1..Length(coeffs)] do 383 L[i]:= (Valuation(coeffs[i]))/i; 384 od; 385 res2 := Minimum(L); 386 return res2; 387 end; 388 389 390 #//////////////////////////////////////////// 391 #// computation of characteristic polynomials 392 393 #//////////////////////////////////////////////////////////////// 394 #// compute the characteristic polynomials using newton relations 395 #// (cohen 1, p 161) 396 #// the traces_from_poly of tyhe reference polynomials Phi are 397 #// precomputed and applied in th computation of char polys wrt 398 #// Phi. 399 400 traces_from_poly := function(Phi) 401 local n,R,L,inc,tmp_R,tmp_PR,tmp_Phi,t,S,k,Sk,i,char; 402 n := Degree(Phi); 403 R := CoefficientRing(Parent(Phi)); 404 if Is(Type(R),ord^ser) or Is(Type(R),res^ser) then 405 char := Coerce(R,Characteristic(CoefficientRing(R))); 406 else char:= -1; 407 fi; 408 #// precision loss occurs in poly_from_traces by division 409 if n > 0 then 410 L:=[]; 411 for i in [1..n] do 412 if Coerce(R,i) = char then 413 L[i]:= 0; 414 else L[i]:= Valuation(Coerce(R,i)); 415 fi; 416 od; 417 inc := Maximum(L); 418 inc := inc^2; 419 else 420 inc := 0; 421 fi; 422 if not Precision(R)=INFTY then 423 tmp_R := ChangePrecision(R,Precision(R)+inc); 424 else tmp_R:= R; 425 fi; 426 tmp_PR := PolynomialAlgebra(tmp_R); 427 tmp_Phi := Coerce(tmp_PR,Phi); 428 429 t := function(i) 430 if i < 0 then 431 return 0; 432 else 433 return Coefficient(tmp_Phi,i); 434 fi; 435 end; 436 437 S := []; 438 439 for k in [1..n] do 440 Sk := -k*t(n-k); 441 for i in [1..Minimum(n,k-1)] do 442 Sk := Sk - t(n-i)*S[k-i]; 443 od; 444 Append_(S,[Sk]); 445 od; 446 return S; 447 end; 448 449 450 poly_from_traces := function(S) 451 local n,Tn,k,Tnk,i; 452 n := Length(S); 453 Tn := []; 454 455 for k in [1..n] do 456 Tnk := -S[k]; 457 for i in[1..k-1] do 458 Tnk := Tnk - Tn[i]*S[k-i]; 459 od; 460 Tnk := Div(Tnk,k); 461 Append_(Tn,[Tnk]); 462 od; 463 Tn := Reverse(Tn); 464 Append_(Tn,[Coerce(Parent(S[1]),1)]); 465 return Coerce(PolynomialAlgebra(Parent(S[1])),Tn); 466 end; 467 468 469 chi := function(Phi, traces_Phi, theta) 470 local ret,k, L,tmp_PR,tmp_theta,tmp_Phi,n,S,pow2,i,Si,j,M,S,li,T,nu; 471 472 tmp_PR := PolynomialAlgebra(Parent(traces_Phi[1])); 473 #L:= ElementToSequence(theta); 474 #Print(theta,"\n"); 475 #Print(L,"\n",Type(L[2]),"\n"); 476 #M:= []; 477 # for i in [1..Length(L)] do 478 # nu:= ElementToSequence(L[i]); 479 480 # if Length(nu)>1 then 481 # li:= false; 482 # for j in [2..Length(nu)] do 483 # if nu[j] <>0 then 484 # li:= true; 485 # fi; 486 # od; 487 # if li then 488 # T:= []; 489 # for j in [1..Length(nu)] do 490 # T[j]:=Coerce(IntegerRing(),nu[j]); 491 # od; 492 # M[i]:= Element(CoefficientRing(tmp_PR),T); 493 # else 494 # M[i]:=Coerce(IntegerRing(),L[i]); 495 #fi; 496 # else 497 # M[i]:=Coerce(IntegerRing(),L[i]); 498 # fi; 499 #od; 500 501 tmp_theta := Coerce(tmp_PR,theta); 502 tmp_Phi := Coerce(tmp_PR,Phi); 503 #Print("tmp_PR: ",tmp_PR,"\nCoefficientRing(tmp_PR): ",CoefficientRing(tmp_PR),"\nType(theta): ",Type(theta),"\nCoefficientRing(Parent(theta)): ",CoefficientRing(Parent(theta)),"\n"); 504 n := Degree(Phi); 505 S := []; 506 pow2 := 1; 507 for i in [1..n] do 508 pow2 := pow2*tmp_theta mod tmp_Phi; 509 510 Si := 0; 511 Si := n*Coefficient(pow2,0); 512 for j in [1..Degree(pow2)] do 513 if Coefficient(pow2,j) <> 0 then 514 Si :=Si+traces_Phi[j]*Coefficient(pow2,j); 515 fi; 516 od; 517 Append_(S,[Si]); 518 od; 519 ret := Coerce(Parent(Phi),poly_from_traces(S)); 520 return ret; 521 end; 522 523 pow_chi:= function(Phi,theta) 524 local A2,chi,L,j,A,chir,a,b, elseq,pi,poly,prec,n,i,j,E,ech,phi,pa,T,nu,ra,k,loss,p,inv,denv,mult; 525 phi:= theta.elt; 526 denv:= theta.denval; 527 pa:= Parent(Phi); 528 prec:= Precision(CoefficientRing(pa)); 529 pi:= UniformizingElement(CoefficientRing(pa)); 530 n:= Degree(Phi); 531 A:= []; 532 for i in [1..n+1] do 533 a:= pi^(denv*(n-i+1)); 534 b:= phi^(i-1); 535 a:= Coerce(Parent(b),a); 536 elseq:= ElementToSequence(a*b mod Phi); 537 for j in [Length(elseq)+1..n+1] do 538 Append_(elseq,[Coerce(CoefficientRing(pa),0)]); 539 od; 540 A:=Append(A,elseq); 541 od; 542 for i in [1..n+1] do 543 L:=[]; 544 for j in [1..n+1] do 545 if i = j then 546 L[j]:= pi^prec; 547 else L[j]:= Coerce(pa,0); 548 fi; 549 od; 550 A:=Append(A,L); 551 od; 552 A:= Matrix(CoefficientRing(pa),2*(n+1),n+1,A); 553 ech:= EchelonForm(A); 554 E:= ech.base; 555 T:= ech.ext1; 556 ra:= Rank(E); 557 for i in [ra..2*(n+1)] do 558 nu:= true; 559 for j in [1..Length(E[i])] do 560 if E[i][j]<> Coerce(Parent(E[i][j]),0) then 561 nu:=false; 562 j:= Length(E[i]); 563 fi; 564 od; 565 if nu then chi:= Coerce(Parent(Phi),List(T[i])); 566 if chi <> Element(Parent(Phi),0) then 567 chir:=chi; 568 fi; 569 fi; 570 od; 571 elseq:= ElementToSequence(chir); 572 k:= []; 573 poly:=[]; 574 for i in [1..Length(elseq)] do 575 poly[i]:= elseq[i]; 576 k[i]:= Valuation(poly[i]); 577 od; 578 loss:=Minimum(k); 579 p:= pi^loss; 580 chir:= Element(Parent(Phi),poly); 581 chir:= Div(chir,p); 582 inv:= IsInvertible(LeadingCoefficient(chir)); 583 if inv.base then 584 inv:= inv.ext1; 585 else return Error("pow_chi not monic"); 586 fi; 587 chir:= chir*Coerce(pa,inv); 588 return [chir,loss]; 589 end; 590 591 592 #{Find char poly of theta * pi^-theta_denval} 593 elt_chi := function(Phi, traces_Phi, theta) 594 local PR,R,tmp_R,tmp_PR,res2,pi; 595 if Is(Type(Phi),elt-alg^pol/ord^ser) or Is(Type(Phi),elt-alg^pol/res^ser) then 596 return pow_chi(Phi,theta)[1]; 597 fi; 598 #vprint RoundFour,5:" -- char poly"; 599 #vprint RoundFour,6:" -- of",theta`elt; 600 #vprint RoundFour,6:" -- with denominator pi ^",theta`denval; 601 PR := Parent(Phi); 602 R := CoefficientRing(Parent(Phi)); 603 if not Precision(R)=INFTY then 604 tmp_R := ChangePrecision(R,Precision(R)+4*theta.denval*Degree(Phi)); 605 else tmp_R:= R; 606 fi; 607 tmp_PR := PolynomialAlgebra(tmp_R); 608 #vtime RoundFour,5: 609 #Print("Phi: ",Phi,"\nCoefficientRing(Parent(Phi)):",CoefficientRing(Parent(Phi)),"\n"); 610 #Print("traces_Phi: ",traces_Phi,"\nType(traces_Phi[1]): ",Type(traces_Phi[1]),"\n"); 611 #Print("theta.elt:",theta.elt,"\nCoefficientRing(Parent(theta.elt)):",CoefficientRing(Parent(theta.elt)),"\n\n"); 612 res2 := chi(Phi, traces_Phi, theta.elt); 613 #//"vstar",v_star(res); 614 #//"denval",theta`denval; 615 if v_star(res2) < theta.denval then 616 return Error("Factorization:Insufficient precision in characteristic polynomial computation"); 617 fi; 618 res2 := Coerce(tmp_PR,res2); 619 pi := UniformizingElement(tmp_R); 620 res2 := Evaluate(res2, pi^theta.denval * tmp_PR.1); 621 res2:= Div(res2,LeadingCoefficient(res2)); 622 res2 := Coerce(PR,res2); 623 if IsMonic(res2)= false then 624 return Error("f not monic at elt_chi\n"); 625 return; 626 fi; 627 #Print(" testing chi(theta)", elt_evaluate(Phi,eltRec(Coerce(PR,res2),0),theta),"\n"); 628 return Coerce(PR,res2); 629 end; 630 631 632 #// compute invariants of elt, war eigentlich procedure, die phi veraenderte 633 elt_invariants:= function(Phi, traces_Phi, phi) 634 local res_chi,fact,K,pc; 635 636 if phi.chi="null" then 637 if Is(Type(Phi),elt-alg^pol/ord^ser) or Is(Type(Phi),elt-alg^pol/res^ser) then 638 pc:= pow_chi(Phi,phi); 639 phi.chi:= pc[1]; 640 phi.chi:= Coerce(PolynomialAlgebra(ChangePrecision(CoefficientRing(Parent(pc[1])),Precision(CoefficientRing(Parent(pc[1])))+pc[2])),phi.chi); 641 else 642 phi.chi := elt_chi(Phi, traces_Phi, phi); 643 fi; 644 fi; 645 phi.v_star := v_star(phi.chi); 646 #Print("v_star:",phi.v_star," chi:",phi.chi,"\n"); 647 #Print("CoefficientRing(Parent(chi)): ", CoefficientRing(Parent(phi.chi)),"\n"); 648 phi.passes_newton_test := 649 phi.v_star = Valuation(TrailingCoefficient(phi.chi)) / Degree(phi.chi); 650 phi.e := Denominator(phi.v_star); 651 if phi.passes_newton_test or phi.v_star = 0 then 652 if phi.res_fact = "null" then 653 if IsField(CoefficientRing(Parent(phi.chi))) then 654 K:=ResidueClassField(IntegerRing(CoefficientRing(Parent(phi.chi)))); 655 else 656 K := ResidueClassField(CoefficientRing(Parent(phi.chi))); 657 fi; 658 659 #Print("K: ",K,"phi.chi: ", phi.chi,"\n"); 660 res_chi := Coerce(PolynomialAlgebra(K),phi.chi); 661 # Print("Parent(res_chi): ",Parent(res_chi),"res_chi: ",res_chi,"\n"); 662 fact := Factorization(res_chi); 663 #Print("fact: ",fact," Parent(fact[1][1]): ",Parent(fact[1][1]),"\n"); 664 phi.res_fact := fact; 665 fi; 666 phi.passes_hensel_test := (Length(phi.res_fact) = 1); 667 if phi.passes_hensel_test then 668 phi.nu := Coerce(Parent(phi.chi),phi.res_fact[1][1]); 669 phi.f := Degree(phi.nu); 670 fi; 671 else 672 phi.passes_hensel_test := true; 673 fi; 674 return phi; 675 end; 676 677 #// ??? !!! TODO ft(k2,L[14][3]); 678 679 #// returns theta`elt(tau`elt/pi^tau`denval)/pi^theta`denval mod Phi 680 elt_evaluate := function(Phi, theta, tau) 681 local PR,R,inc,tmp_R,tmp_PR,pi,tmp_Phi,tmp_theta,tmp_tau,res2,L,se,i,res_denval; 682 PR := Parent(Phi); 683 R := CoefficientRing(PR); 684 inc := Degree(theta.elt)*tau.denval; 685 if not Precision(R)=INFTY then 686 tmp_R := ChangePrecision(R,Precision(R)+inc); 687 else tmp_R:= R; 688 fi; 689 tmp_PR := PolynomialAlgebra(tmp_R); 690 pi := UniformizingElement(tmp_R); 691 tmp_Phi := Coerce(tmp_PR,Phi); 692 tmp_theta :=Coerce(tmp_PR,theta.elt); 693 tmp_tau := Coerce(tmp_PR,tau.elt); 694 695 696 res2 := tmp_theta*pi^inc; 697 L:= []; 698 se:= ElementToSequence(res2); 699 for i in [1..Degree(res2)+1] do 700 L[i]:= Div(se[i],pi^(tau.denval*(i-1))); 701 od; 702 res2 := Coerce(tmp_PR,L); 703 res2 := Evaluate(res2, tmp_tau) mod tmp_Phi; 704 res_denval := theta.denval+inc; 705 return eltRec(Coerce(PR,res2), res_denval); 706 end; 707 708 #// elt_evaluate 709 #// and transfer invariants 710 #// useful when changing reference polynomial 711 elt_evaluate_trans := function(Phi, traces_Phi, theta, tau) 712 local res2; 713 714 res2 := elt_evaluate(Phi, theta, tau); 715 if theta.chi <> "null" then 716 res2.chi := theta.chi; 717 fi; 718 if theta.res_fact <>"null" then 719 res2.res_fact := theta.res_fact; 720 fi; 721 res2:=elt_invariants(Phi, traces_Phi, res2); 722 return res2; 723 end; 724 725 726 727 #////////////////////////////////////////////////////////////////////////// 728 #// 729 #// returns a list of exponents res such that 730 #// v*(pi^res[1]*pool[1]^res[2]*...*pool[m]^res[m+1]) eq val 731 #// where res[i]<ceils[i-1] (for i>=2) and 732 #// Denominator(val) | E 733 calculate_exponents := function(pool, ceils, val, E) 734 local res2,rem_val,len,now,d,tmpval,tmpool,mud,rcr,sol; 735 res2 := []; 736 rem_val := val; 737 len := Length(pool); 738 now := len; 739 740 d := E; 741 742 while now > 0 do 743 744 d := d/(ceils[now]); 745 tmpval := rem_val * d; 746 tmpool := pool[now] *d; 747 mud := Denominator(tmpool); 748 rcr := ResidueClassRing(Coerce(IntegerRing(),mud)); 749 tmpval := tmpval * mud; 750 tmpool := tmpool * mud; 751 sol := Coerce(IntegerRing(),((Coerce(rcr,tmpval))/(Coerce(rcr,tmpool)))); 752 753 rem_val := rem_val - sol * pool[now]; 754 Append_(res2, [sol]); 755 756 now := now - 1; 757 758 od; 759 760 Reverse_(res2); 761 Insert_(res2,rem_val,1); 762 return res2; 763 end; 764 765 #//{Returns an element poly, such that 766 #// poly has v_star valuation vstar} 767 elt_with_v_star := function(Phi, P, P_v_star, S, E, vstar) 768 local exps,poly_val,poly,i,pi; 769 770 exps := calculate_exponents(P_v_star, S, vstar, E); 771 poly_val := Coerce(IntegerRing(),exps[1]); 772 poly := Coerce(Parent(Phi),1); 773 if Length(P) <> 0 then 774 for i in [1..Length(P)] do 775 poly:= poly*P[i] ^ exps[i + 1]; 776 od; 777 fi; 778 779 780 if poly_val > 0 then 781 pi := UniformizingElement(CoefficientRing(Parent(Phi))); 782 poly := poly*pi^poly_val; 783 poly_val := 0; 784 fi; 785 return eltRec(poly, -poly_val); 786 end; 787 788 #////////////////////////////////////////////////////////////////////////// 789 #// the denominator of theta is pi^theta`denval 790 #// if theta`chi is not assigned or theta`denval gt 0 then chi_theta is computed 791 792 hensel_lift_factorization_sub := function(Phi, theta) 793 local PR,R,tmp_R,tmp_PR,K,P,pi,tmp_Phi,traces_tmp,Phi,tmp_theta,facts,i,j,lifted_facts,lf, eval_facts,res2,traces_tmp_Phi,elm,L,elmi; 794 PR := Parent(Phi); 795 R := CoefficientRing(Parent(Phi)); 796 if not Precision(R) =INFTY then 797 tmp_R := ChangePrecision(R,Precision(R)+Degree(Phi)*(theta.denval)); 798 else tmp_R:= R; 799 fi; 800 tmp_PR := PolynomialAlgebra(tmp_R); 801 K := ResidueClassField(tmp_R); 802 P := PolynomialAlgebra(K); 803 pi := UniformizingElement(tmp_R); 804 tmp_Phi := Coerce(tmp_PR,Phi); 805 traces_tmp_Phi := traces_from_poly(tmp_Phi); 806 tmp_theta := theta; 807 if theta.denval <> 0 or theta.chi ="null" then 808 tmp_theta.chi := elt_chi(tmp_Phi, traces_tmp_Phi, tmp_theta); 809 fi; 810 #Print(" testing chi(theta)", elt_evaluate(tmp_Phi,eltRec(tmp_theta.chi,0),tmp_theta),"\n"); 811 if theta.res_fact <>"null" then 812 tmp_theta.res_fact := theta.res_fact; 813 fi; 814 tmp_theta := elt_invariants(tmp_Phi, traces_tmp_Phi,tmp_theta); 815 facts:=[]; 816 if tmp_theta.res_fact<>"null" then 817 for i in [1..Length(tmp_theta.res_fact)] do 818 facts[i]:=Coerce(tmp_PR,Coerce(P,(tmp_theta.res_fact[i][1]^tmp_theta.res_fact[i][2]))); 819 od; 820 fi; 821 if Length(facts)=1 then 822 if facts[1] = tmp_theta.chi then 823 L:=[]; 824 L[1]:= [tmp_theta.chi,0]; 825 return L; 826 else 827 facts:=Append(facts,[Coerce(Parent(facts[1]),1)]); 828 fi; 829 fi; 830 #Print(facts); 831 lifted_facts := HenselLift(tmp_theta.chi, facts); 832 lf:=lifted_facts[1]; 833 for i in [2..Length(lifted_facts)] do 834 lf:= lf*lifted_facts[i]; 835 od; 836 if Is(Type(R),ord^ser) or Is(Type(R),res^ser) then 837 # assuming precision 20 838 elm:=(ElementToSequence(lf-tmp_theta.chi)); 839 if elm <> [] then 840 for i in [1..Length(elm)] do 841 elmi:=ElementToSequence(elm[i]); 842 for j in [1.._realPrec] do 843 if elmi[j]<>Element(Parent(elmi[i]),0) then return Error("assertion violated\n");fi; 844 od; 845 od; 846 fi; 847 else if (lf- tmp_theta.chi)<>0 then return Error("assertion violated\n"); fi; 848 fi; 849 850 eval_facts:= [1..Length(lifted_facts)]; 851 #Print("henseltest lifted_facts ",lifted_facts); 852 Apply_(eval_facts,x->elt_evaluate(tmp_Phi,eltRec(lifted_facts[x],0),tmp_theta)); 853 #Print("henseltest eval_facts ",eval_facts,"\n"); 854 #vprint RoundFour,4: " -- gcd"; 855 #vtime RoundFour,4: 856 res2:=[1..Length(eval_facts)]; 857 Apply_(res2,x->PGcd(Phi,Coerce(PR,(eval_facts[x].elt)))); 858 return res2; 859 end; 860 861 862 #// the denominator of theta is pi^theta_denval 863 #// if chi_theta eq 0 or theta_denval gt 0 then chi_theta is computed 864 hensel_lift_factorization_chi2 := function(Phi, theta) 865 local flag,PR,R,res2,inc,resProd,i,tmp_R,tmp_theta,res2,tmp_PR,re; 866 #Print("hensel_lift_factorization_chi2"); 867 flag := 1; 868 869 PR := Parent(Phi); 870 R := CoefficientRing(Parent(Phi)); 871 # vprint RoundFour,5:" == Hensel lifting to precision", Precision(R),"..."; 872 res2 := hensel_lift_factorization_sub (Phi, theta); 873 #// if the precion was not sufficient, repeat this step with a higher 874 #// precision. 875 inc := 1; 876 resProd:= res2[1][1]; 877 while (IsList(resProd)) do 878 resProd:= resProd[1]; 879 od; 880 for i in[2..Length(res2)] do 881 if IsList(res2[i]) then 882 re:=res2[i][1]; 883 while IsList(re) do 884 re:= re[1]; 885 od; 886 resProd:= resProd*re; 887 fi; 888 od; 889 #Print("\nResProd ",resProd,"\n"); 890 while (flag = 1 and Element(Parent(Phi),0) <>((resProd)-Phi)) do 891 #Print("difference: ",resProd-Phi,"\n"); 892 #//print inc; 893 if not Precision(R) =INFTY then 894 tmp_R := ChangePrecision(R,Precision(R)+4*Degree(Phi)*inc); 895 else tmp_R:= R; 896 fi; 897 #vprint RoundFour,5:" == ... trying Hensel lift again, precision increased by", 898 # 4*Degree(Phi)*inc,"..."; 899 tmp_PR := PolynomialAlgebra(tmp_R); 900 tmp_theta := eltRec(Coerce(tmp_PR,theta.elt), theta.denval); 901 902 if theta.res_fact<>"null" then 903 904 tmp_theta.res_fact := theta.res_fact; 905 fi; 906 #Print("henselLift: prec ",Precision(tmp_R)); 907 res2:= hensel_lift_factorization_sub(Coerce(tmp_PR,Phi),tmp_theta); 908 for i in [1..Length(res2)] do 909 res2[i][1] := [Coerce(PR,res2[i][1])]; 910 od; 911 resProd:= res2[1][1]; 912 while (IsList(resProd)) do 913 resProd:= resProd[1]; 914 od; 915 for i in[2..Length(res2)] do 916 re:=res2[i][1]; 917 while(IsList(re)) do 918 re:= re[1]; 919 od; 920 resProd:= resProd*re; 921 od; 922 inc :=inc+ 1; 923 od; 924 #vprint RoundFour,5:" == ... done"; 925 return res2; 926 end; 927 928 929 newton_lift_factorization_chi2 := function(Phi, theta) 930 local rs,gamma; 931 # //print "newton theta",theta,theta_denval; 932 #//"theta",theta; 933 if theta.v_star = 0 and theta.passes_hensel_test=FALSE then 934 #//print "newton_lift_factorization_chi2 hensel"; 935 return hensel_lift_factorization_chi2(Phi, theta); 936 fi; 937 938 rs := theta.v_star; 939 #vprint RoundFour,5:" == flattening Newton polygon slope from",rs,"to 0"; 940 gamma := eltRec(Coerce(Parent(Phi),(theta.elt^Denominator(rs)) mod Phi), 941 Numerator(rs)+theta.denval*Denominator(rs)); 942 return hensel_lift_factorization_chi2(Phi, gamma); 943 944 end; 945 946 #///////////////////////////////////////////////////////////////////////////////// 947 948 #//{Find a polynomial delta with Degree(delta)<Degree(nu_phi) 949 #// and v_star(delta-gamma)>0} 950 #// returns: 951 #// true, delta, xxx 952 #// if such an element delta has been found 953 #// false, xxx, rho 954 #// if during the search for delta an element that fails the 955 #// hensel test has occured. 956 957 res_field_rep := function (Phi, traces_Phi, nu, gamma) 958 local pi,K,P,nur,F,PF,facts,i,a,b,delta,rho,Co,p,prec,g,L,M,S,T,j,U,li; 959 960 if Degree(gamma.nu) = 1 then 961 return [true, -TrailingCoefficient(gamma.nu), 0]; 962 fi; 963 964 if gamma.denval = 0 then 965 return [true, gamma mod nu, 0]; 966 fi; 967 Co:= CoefficientRing(Parent(Phi)); 968 969 # if Is(Type(Co),res^pad) then 970 # p:= Size(ResidueClassField(Co).base); 971 # prec:= Coerce(IntegerRing(),(Precision(Co)/RamificationDegree(Co))); 972 # g:= DefiningPolynomial(Co); 973 # Co:= pAdicRing(p,prec); 974 # Co:= TotallyRamifiedExtension(Co,g); 975 # L:= ElementToSequence(g); 976 977 # M:= []; 978 # for i in [1..Length(L)] do 979 # S:= ElementToSequence(L[i]); 980 # if Length(S)>1 then 981 # li:= false; 982 # for i in [2..Length(S)] do 983 # if S[i] <>0 then 984 # li:= true; 985 # fi; 986 # od; 987 # if li then 988 # T:= []; 989 # for j in [1..Length(S)] do 990 # T[i]:=Coerce(IntegerRing(),S[i]); 991 # od; 992 # M[i]:= Element(Co,T); 993 # else 994 # M[i]:=Coerce(IntegerRing(),L[i]); 995 # fi; 996 # else M[i]:=Coerce(IntegerRing(),L[i]); 997 # fi; 998 # od; 999 # gamma.elt:= Coerce(PolynomialAlgebra(Co),M); 1000 # fi; 1001 pi := UniformizingElement(Co); 1002 K := ResidueClassField(Co); 1003 P := PolynomialAlgebra(K); 1004 nur := Coerce(P,nu); 1005 F := Extension(K,nur); 1006 PF:= PolynomialAlgebra(F); 1007 facts := Roots(Coerce(PF,(Coerce(P,gamma.nu)))); 1008 for i in [1..Length(facts)] do 1009 a := ElementToSequence(facts[i][1]); 1010 b := Coerce(PF,a); 1011 delta := Coerce(PolynomialAlgebra(Co),b); 1012 rho := eltRec(gamma.elt-Coerce(Parent(gamma.elt),pi)^gamma.denval*delta,gamma.denval); 1013 rho:=elt_invariants(Phi, traces_Phi, rho); 1014 if rho.passes_hensel_test=FALSE then 1015 return [false, 0, rho]; 1016 fi; 1017 if rho.v_star > 0 then 1018 return [true, Coerce(Parent(Phi),delta), 0]; 1019 fi; 1020 od; 1021 return Error ("Factorization: Round Four Panic, no representing element found."); 1022 end; 1023 1024 #// x gives a residue class field extension of degree F_x 1025 #// gamma/pi^gamma_denval gives a residue class field extension of 1026 #// degree F_gamma 1027 #// find a poltynomial that gives a residue class field extension of 1028 #// degree LCM(F_x,F_gamma) 1029 # //print "kappa",kappa,"/",pi^kappa_denval; 1030 1031 ff_composition_gen := function(Phi, x, F_x, gamma) 1032 local kappa,PR,R,denval,tmp_R,tmp_PR,x,pi,tmp_Phi,traces_tmp_Phi,tmp_gamma,K,F,pow,pow2,coeff,li,y,i; 1033 if gamma.f = LCM(F_x,gamma.f) then 1034 kappa := gamma; 1035 return kappa; 1036 fi; 1037 1038 #//print "gamma",gamma; 1039 PR := Parent(Phi); 1040 R := CoefficientRing(Parent(Phi)); 1041 denval := gamma.denval*gamma.f; 1042 if not Precision(R) = INFTY then 1043 tmp_R := ChangePrecision(R,Precision(R)+denval); 1044 else tmp_R:= R; 1045 fi; 1046 #//tmp_R := pAdicQuotientRing(Prime(R),Precision(R)+denval); 1047 tmp_PR := PolynomialAlgebra(tmp_R); 1048 y:=tmp_PR.1; 1049 pi := UniformizingElement(tmp_R); 1050 tmp_Phi := Coerce(tmp_PR,Phi); 1051 traces_tmp_Phi := traces_from_poly(tmp_Phi); 1052 tmp_gamma := Coerce(tmp_PR,gamma.elt); 1053 1054 K := ResidueClassField(tmp_R); 1055 1056 F := LCM(F_x,gamma.f); 1057 pow:=[1..F_x]; 1058 Apply_(pow,x->(pi^(denval)*y^(x-1)) mod (pi^denval*tmp_Phi)); 1059 pow2:=[1..gamma.f]; 1060 Apply_(pow2,x->pi^(denval-(gamma.denval*(x-1)))*tmp_gamma^(x-1)mod (pi^denval*tmp_Phi)); 1061 pow := Concatenation(pow,pow2); 1062 1063 repeat 1064 coeff:=[1..Length(pow)]; 1065 Apply_(coeff,x->Coerce(tmp_R,Random(K))); 1066 #//"coeff",coeff; 1067 li:= [1..Length(pow)]; 1068 Apply_(li,x->coeff[x]*pow[x]); 1069 kappa:= li[1]; 1070 for i in [2..Length(li)] do 1071 kappa := kappa+li[i]; 1072 od; 1073 kappa:= eltRec(kappa mod tmp_Phi,denval); 1074 kappa:= elt_invariants(tmp_Phi, traces_tmp_Phi,kappa); 1075 until kappa.f = F; 1076 1077 kappa.elt := Coerce(PR,kappa.elt); 1078 1079 return kappa; 1080 1081 end; 1082 1083 1084 1085 1086 #//////////////////////////////////////////////////////////////////////// 1087 1088 1089 #//------------------------- 1090 round_four := function(Phi, with_certificates) 1091 local PR,K,pi,Phihat,traces_Phi,traces_Phihat,P,S,P_v_star,E,phi_old_v_star,phi,F,kappa,nu,Pi,theta, 1092 s,Gamma,BETA,phis,poly,psi,poly_inv,gamma,beta,newPhi,r_L,rho_passes_hensel_test,delta,rho, deltapsi,further_factorization; 1093 #//------------------------- 1094 #// the algorithm used is a hybrid between 1095 #// pauli: factoring polynomials over local fields 1096 #// and 1097 #// ford-pauli-roblot: a fast algorithm for factoring polynomials over Qp 1098 #// 1099 #// step g) [extend the ground field] in the pauli algorithm is replaced by 1100 #// by the combination of step 12 in algorithm 5.1 and step 2,3 in 5.2 in 1101 #// ford-pauli-roblot 1102 #// 1103 #// if certificates is false the information needed for the certificates is not 1104 #// computed 1105 1106 1107 1108 # vprint RoundFour,1:"R4: factoring Phi with deg(Phi) =",Degree(Phi); 1109 # vprint RoundFour,3:"R4: Phi =",Phi; 1110 1111 # FACTOR := recformat 1112 # <factor:RngUPolElt, 1113 # multiplicity:IntegerRing(), 1114 # F:IntegerRing(), // inertia deg 1115 # Gamma:RngUPolElt, // unramified 1116 # Gamma_denval:IntegerRing(), // certificate 1117 # E:IntegerRing(), // ram index 1118 # Pi:RngUPolElt, // ramified 1119 # Pi_denval:IntegerRing(), // certificate 1120 # IdealGen1, // two element representation 1121 # IdealGen2 // of the prime ideal 1122 # >; 1123 1124 further_factorization := function(factors, with_certificates) 1125 local L,i,certs,j,fact; 1126 L := []; 1127 for i in [1..Length(factors)] do 1128 fact:=factors[i][1]; 1129 while IsList(fact) do 1130 fact:= fact[1]; 1131 od; 1132 certs:=round_four(fact, with_certificates); 1133 for j in [1..Length(certs)] do 1134 Append_(L,[certs[j]]); 1135 od; 1136 od; 1137 return L; 1138 end; 1139 1140 PR := Parent(Phi); 1141 K := CoefficientRing(PR); 1142 pi := UniformizingElement(K); 1143 1144 # vprint RoundFour,4:"R4: over",K; 1145 1146 if Degree(Phi) <= 1 then 1147 return [rec(factor:=Phi,multiplicity:= 1, F:=1,Gamma:=Coerce(PR,1),Gamma_denval:=0, E:=1,Pi:=Coerce(PR,pi),Pi_denval:=0)]; 1148 fi; 1149 1150 1151 Phihat := Phi; 1152 traces_Phi := traces_from_poly(Phi); 1153 traces_Phihat := traces_Phi; 1154 #//////////////////////////// RAMIFIED DATA 1155 P := []; #// list of polynomials that allows 1156 #// construction of polynomials with 1157 #// v_star valuation with den up to E 1158 S := []; #// list of maximal exponents to be used 1159 #// for the polynomials in P 1160 P_v_star := []; #// list of v_star valuations of the elements 1161 #// of P 1162 E := 1; #// highest known ramification index 1163 phi_old_v_star := 0; #// old v* valuation of phi, used to check 1164 #// that valuation of phi increases 1165 #//////////////////////////// 1166 BETA := eltRec(PR.1,0); #// we record the change of the reference 1167 #// polynomial Phihat in BETA/pi^BETA_denval 1168 #//////////////////////////// 1169 1170 #////////////////////////////////////////////////////////////////////////// 1171 1172 1173 1174 1175 1176 phi := eltRec(PR.1,0); 1177 phi.chi := Phihat; 1178 phi:=elt_invariants(Phihat, traces_Phihat, phi); 1179 1180 if phi.passes_newton_test=FALSE or phi.passes_hensel_test=FALSE then 1181 #vprint RoundFour,1:"R4: x fails Newton or Hensel test -- factoring"; 1182 return further_factorization(newton_lift_factorization_chi2(Phi, phi), with_certificates); 1183 fi; 1184 1185 #///////////////////// UNRAMIFIED DATA 1186 F := phi.f; #// highest known inertia degree 1187 kappa := phi; #// generator of the unramified extension of degree F 1188 nu := kappa.nu; #// generating polynomial of the unramified 1189 #///////////////////// extension of degree F 1190 1191 phi := eltRec(phi.nu,0); 1192 phi:=elt_invariants(Phihat, traces_Phihat,phi); 1193 1194 #vprint RoundFour,1:"R4: initializing E to",E,"and F to",F; 1195 #vprint RoundFour,3:"R4: inertia element:",BETA`elt; 1196 1197 if F = Degree(Phi) then 1198 #vprint RoundFour,1:"R4: F =",F,"= deg(Phi) => Phi is irreducible"; 1199 #//vprint RoundFour,1:"R4: computing 2 element certificate"; 1200 return [rec(factor:=Phi,multiplicity:=1,F:=F,Gamma:=PR.1,Gamma_denval:=0,E:=1, Pi:=Coerce(PR,pi), Pi_denval:= 0)]; 1201 fi; 1202 1203 #//print "nu_phi",nu_phi; 1204 1205 1206 #// main loop 1207 while true do 1208 # vprint RoundFour,2:"R4| phi =",phi`elt; 1209 #vprint RoundFour,1:"R4| deg(phi) =", Degree(phi`elt),"and v*(phi) =", phi`v_star; 1210 # vprint RoundFour,2:"R4| E =",E,"and F =",F; 1211 #/* Step a ***************************************************************/ 1212 #Print(phi.v_star," ",phi_old_v_star,"\n"); 1213 if phi.v_star <= phi_old_v_star then 1214 return Error ("Factorization: Insufficient precision in Round Four main loop"); 1215 fi; 1216 if Degree(phi.elt) > E*F then 1217 return Error ("Factorization: Degree of phi too large in Round Four main loop"); 1218 fi; 1219 1220 if phi.passes_newton_test=FALSE then 1221 #vprint RoundFour,1:"R4a phi fails Newton test -- factoring"; 1222 theta := elt_evaluate_trans(Phi,traces_Phi,phi,BETA); 1223 return further_factorization(newton_lift_factorization_chi2 1224 (Phi, theta),with_certificates); 1225 #vprint RoundFour,4:"R4a phi passes Newton test"; 1226 fi; 1227 1228 #/* Step b [increase E] **************************************************/ 1229 if E mod phi.e <> 0 then 1230 Append_(P, [phi.elt]); 1231 Append_(P_v_star, [v_star(chi(Phihat, traces_Phihat, phi.elt))]); 1232 s := Div(LCM(E, phi.e),E); 1233 Append_(S, [s]); 1234 #vprint RoundFour,1:"R4b increasing E from",E,"to",s*E; 1235 #vprint RoundFour,4:"R4b list of polynomials giving ramification",P; 1236 #vprint RoundFour,4:"R4b with v* valuations",P_v_star; 1237 #vprint RoundFour,4:"R4b and maximal exponents",S; 1238 E := s * E; 1239 1240 if E*F = Degree(Phihat) then #// Phi is irreducible 1241 #vprint RoundFour,1:"R4b E * F =",E,"*",F,"= deg(Phi) =",Degree(Phihat), 1242 # " => Phi is irreducible"; 1243 if with_certificates then 1244 #vprint RoundFour,1:"R4b computing 2 element certificate"; 1245 #//vprint RoundFour,5:"R4b computing ramified certificate"; 1246 Pi := elt_with_v_star(Phihat, P, P_v_star, S, E, 1/E); 1247 Pi := elt_evaluate(Phi, Pi, BETA); 1248 Pi:=elt_invariants(Phi, traces_Phi, Pi); 1249 if Pi.v_star <> 1/E then 1250 return Error ("Factorization: Insufficient precision in computation of ramified certificate"); 1251 fi; 1252 #//vprint RoundFour,5:"R4b computing unramified certificate"; 1253 if F = 1 then 1254 Gamma := eltRec(Coerce(Parent(Phi),1),0); 1255 else 1256 Gamma := BETA; 1257 fi; 1258 Gamma:=elt_invariants(Phi, traces_Phi, Gamma); 1259 #//print "F",F; 1260 if Gamma.f <> F then 1261 return Error ("Factorization: Insufficient precision in computation of unramified certificate"); 1262 fi; 1263 return [rec(factor:=Phi, multiplicity:=1,F:=F, Gamma:=Gamma.elt, Gamma_denval:=Gamma.denval, E:=E, Pi:=Pi.elt,Pi_denval:=Pi.denval)]; 1264 else #// no certificates 1265 return [rec( factor:=Phi, multiplicity:= 1,F:=F,E:=E)]; 1266 fi; 1267 1268 else #// continue in main loop 1269 #// increase valuation of phi, such that deg(phi)=E*F 1270 phis := eltRec((phi.elt)^s,0); phis.v_star := s*phi.v_star; 1271 phi := phis; 1272 #vprint RoundFour,2:"R4b phi =",phi`elt; 1273 #vprint RoundFour,1:"R4b deg(phi) =", Degree(phi`elt), 1274 # "and v*(phi) =", phi`v_star; 1275 fi; 1276 fi; 1277 1278 #/* Step c **************************************************************/ 1279 #vprint RoundFour,5:"R4c find psi with v*(psi) = v*(phi) =",phi`v_star; 1280 poly := elt_with_v_star(Phi, P, P_v_star, S, E, phi.v_star); 1281 psi := poly.elt * pi^-poly.denval; 1282 1283 #/* Step d *************************************************************/ 1284 #vprint RoundFour,5:"R4d gamma := phi/psi"; 1285 poly_inv := elt_inverse(poly, Phihat); 1286 1287 gamma := eltRec((phi.elt * poly_inv.elt) mod Phihat, poly_inv.denval-poly.denval); 1288 gamma.elt := gamma.elt mod pi^(gamma.denval+1); 1289 gamma:=elt_invariants(Phihat, traces_Phihat, gamma); 1290 1291 if gamma.v_star <> 0 then 1292 #//print vstar_phi,old_vstar_phi; 1293 return Error ("Factorization: Insufficient precision in Round Four main loop"); 1294 fi; 1295 1296 #/* Step e *************************************************************/ 1297 if gamma.passes_hensel_test=FALSE then 1298 #vprint RoundFour,1:"R4e gamma fails Hensel test -- factoring"; 1299 theta := elt_evaluate_trans(Phi,traces_Phi,gamma,BETA); 1300 return further_factorization(hensel_lift_factorization_chi2 1301 (Phi, theta), with_certificates); 1302 #else 1303 #Print ("RoundFour,4:R4e gamma passes Hensel test"); 1304 fi; 1305 1306 #/* Step f ************************************************************/ 1307 #//vprint RoundFour,5:"R4f new inertia by gamma ?"; 1308 if F mod gamma.f <> 0 then #// increase F 1309 #// change reference polynomial Phihat, such that 1310 #// F_x = LCM(F,F_gamma) 1311 #// start the main loop again with Step a 1312 kappa := ff_composition_gen(Phihat, PR.1, F, gamma); 1313 if E * LCM(F, gamma.f) = Degree(Phihat) then #// Phi is irreducible 1314 F := LCM(F, gamma.f); 1315 #vprint RoundFour,1:"R4f E * lcm(F,F_gamma) =",E,"*",F,"= deg(Phi) =",Degree(Phihat), 1316 # " => Phi is irreducible"; 1317 if with_certificates then 1318 #vprint RoundFour,1:"R4f computing 2 element certificate"; 1319 Pi := elt_with_v_star(Phihat, P, P_v_star, S, E, 1/E); 1320 Pi := elt_evaluate(Phi, Pi, BETA); 1321 Pi:=elt_invariants(Phi, traces_Phi, Pi); 1322 if Pi.v_star <> 1/E then 1323 return Error ("Factorization: Insufficient precision in computation of ramified certificate"); 1324 fi; 1325 Gamma := elt_evaluate(Phi, kappa, BETA); 1326 Gamma:=elt_invariants(Phi, traces_Phi, Gamma); 1327 #//print "F",F; 1328 if Gamma.f <> F then 1329 return Error ("Factorization: Insufficient precision in computation of unramified certificate"); 1330 fi; 1331 1332 return [rec(factor := Phi, multiplicity := 1, 1333 F:=F, Gamma := Gamma.elt, Gamma_denval := Gamma.denval, 1334 E:=E, Pi := Pi.elt, Pi_denval := Pi.denval)]; 1335 else 1336 return [rec(factor := Phi, multiplicity := 1,F:=F,E:=E)]; 1337 fi; 1338 #/* Step g ************************************************************/ 1339 else #// find a new reference polynomial Phihat 1340 F := LCM(F, gamma.f); 1341 #vprint RoundFour,1:"R4 increasing F from",F,"to", F, "and resetting E to 1"; 1342 #//"kappa",kappa; 1343 beta := kappa; 1344 #//beta`elt := beta`elt+pi^(beta`denval+1)*PR.1; 1345 #vprint RoundFour,5:"R4g searching beta with F_beta = F =",F, 1346 # "and chi_beta irreducible of Degree deg(Phi) =",Degree(Phi); 1347 newPhi := elt_chi(Phihat, traces_Phihat, beta); 1348 #//while not is_squarefree(newPhi) do 1349 if Discriminant(Phihat) <> 0 then 1350 #// we can use the discriminant for testing whether newPhi is 1351 #// squarefree 1352 while Discriminant(newPhi)= 0 do 1353 #// assure that the new reference polynomial is squarefree 1354 #// then it generates the same extension as Phi 1355 #vprint RoundFour,5:"R4g beta := beta + pi*x"; 1356 beta.elt := beta.elt+pi^(beta.denval+1)*PR.1; 1357 newPhi := elt_chi(Phihat, traces_Phihat, beta); 1358 od; 1359 else 1360 #// better also use the gcd for testing whether newPhi is 1361 #// squarefree as we have a precion problem with the discriminant 1362 #// already 1363 while is_squarefree(newPhi)=FALSE do 1364 #vprint RoundFour,5:"R4g beta := beta + pi*x"; 1365 beta.elt := beta.elt+pi^(beta.denval+1)*PR.1; 1366 newPhi := elt_chi(Phihat, traces_Phihat, beta); 1367 od; 1368 fi; 1369 1370 Phihat := newPhi; 1371 traces_Phihat := traces_from_poly(Phihat); 1372 1373 #// reinitialize 1374 E := 1; 1375 P := []; 1376 S := []; 1377 P_v_star := []; 1378 phi_old_v_star := 0; 1379 1380 BETA := elt_evaluate(Phi, beta, BETA); 1381 #vprint RoundFour,4:"R4g new inertia element beta =",BETA`elt; 1382 #vprint RoundFour,4:"R4g with denominator pi ^",BETA`denval; 1383 1384 phi := eltRec(PR.1,0); 1385 phi.chi := Phihat; 1386 phi:=elt_invariants(Phihat,traces_Phihat,phi); 1387 if (phi.passes_newton_test and phi.passes_hensel_test) =FALSE then 1388 #vprint RoundFour,1:"R4g Newton test or Hensel test failed -- factoring"; 1389 theta := elt_evaluate_trans(Phi,traces_Phi,phi,BETA); 1390 return further_factorization(newton_lift_factorization_chi2 1391 (Phi, theta), with_certificates); 1392 fi; 1393 1394 nu := phi.nu; 1395 phi:= eltRec(phi.nu,0); 1396 phi:=elt_invariants(Phihat,traces_Phihat,phi); 1397 fi; 1398 else 1399 #// no increase of F 1400 #// continue in main loop 1401 #/* Step h *************************************************************/ 1402 # vprint RoundFour,5: 1403 #"R4h find integral delta with chi(delta)=chi(gamma) mod (pi) and deg(delta)<F"; 1404 r_L:= res_field_rep(Phihat,traces_Phihat,nu,gamma); 1405 rho_passes_hensel_test:=r_L[1]; 1406 delta:= r_L[2]; 1407 rho:= r_L[3]; 1408 if rho_passes_hensel_test=FALSE then 1409 #vprint RoundFour,1:"R4h rho fails Hensel test -- factoring"; 1410 theta := elt_evaluate_trans(Phi,traces_Phi,rho,BETA); 1411 return further_factorization(hensel_lift_factorization_chi2 1412 (Phi, theta),with_certificates); 1413 fi; 1414 #vprint RoundFour,3:"R4h delta =",delta; 1415 #/* Step i *************************************************************/ 1416 # vprint RoundFour,5:"R4i phi := phi - delta * psi"; 1417 deltapsi := delta * psi; 1418 #//print "deg", Degree(deltapsi); 1419 phi_old_v_star := phi.v_star; 1420 phi := eltRec(phi.elt-deltapsi,0); 1421 phi:=elt_invariants(Phihat,traces_Phihat,phi); 1422 fi; #// increase F 1423 od; 1424 end; 1425 1426 1427 1428 #//////////////////////////////////////////////////////// 1429 #// auxilary functions for factorization. 1430 1431 #// return true, g, dival, uni, exp such that 1432 #// g := uni*pi^multval*f(pi^exp*x) has integral coefficients and is monic 1433 #// g is over the quotient ring ! 1434 #// if because of insuffient precision such elements cannot be computed, return 1435 #// false, 0, 0, required precision, exp 1436 1437 # // if f is monic and integral do nothing 1438 # //if IsMonic(f) and Minimum([Valuation(a) : a in Eltseq(f)]) ge 0 then 1439 # // 1440 # //return 1441 1442 make_monic_and_integral := function(f) 1443 local PK,K,overfield,R,PR,elseq,vallist,minval,divi,g,multval,PR,pi,lead,vallead,exp,expmult,gList,prec,QR,PQR,uni,SList; 1444 PK := Parent(f); 1445 K := CoefficientRing(PK); 1446 overfield := IsField(K); 1447 if overfield then 1448 R := IntegerRing(K); 1449 PR := PolynomialAlgebra(R); 1450 elseq:= ElementToSequence(f); 1451 vallist:= [1..Length(elseq)]; 1452 Apply_(vallist,x->Valuation(elseq[x])); 1453 minval := Minimum(vallist); 1454 divi := UniformizingElement(K)^(minval); 1455 g := Coerce(PR,(f/divi)); 1456 multval := -minval; 1457 else 1458 R := K; 1459 PR := PK; 1460 elseq:= ElementToSequence(f); 1461 vallist:=[1..Length(elseq)]; 1462 Apply_(vallist,x->Valuation(elseq[x])); 1463 minval := Minimum(vallist); 1464 divi := UniformizingElement(K)^(minval); 1465 g := Coerce(PR,(Div(f,divi))); 1466 multval := -minval; 1467 fi; 1468 #//print g; 1469 PR := PolynomialAlgebra(R); 1470 pi := UniformizingElement(R); 1471 lead := LeadingCoefficient(g); 1472 vallead := Valuation(lead); 1473 elseq := ElementToSequence(g); 1474 exp := -1; 1475 if Degree(g) > 0 then 1476 repeat 1477 exp :=exp+1; 1478 expmult := Degree(g)*exp-vallead; 1479 #//print "expmult",expmult; 1480 vallist:= [2..Length(elseq)]; 1481 Apply_(vallist,i->Valuation(elseq[i])-exp*(i-1)+expmult); 1482 minval:= Minimum(vallist); 1483 until (expmult >= 0 and minval >= 0); 1484 if vallead+expmult >= Precision(R) or exp*Degree(g) >= Precision(R) then 1485 return [false, 0, 0, vallead+expmult+exp*Degree(g), exp]; 1486 fi; 1487 else 1488 exp := 0; 1489 expmult := - vallead; 1490 fi; 1491 g := g * pi^expmult; 1492 #//print "exp",exp; 1493 #//print g; 1494 SList:= ElementToSequence(g); 1495 gList:=[1..(Degree(g)+1)]; 1496 Apply_(gList,i->Div(SList[i],pi^(exp*(i-1)))); 1497 #//print g; 1498 multval := multval+expmult; 1499 if not Is (Type(R),res^pad) and not Is (Type(R),res^ser) then 1500 if Precision(R)=INFTY then prec:= _realPrec; 1501 else prec:= Precision(R); 1502 fi; 1503 if Is(Type(R),ord^ser) or Is(Type(R),ord^pow) then 1504 QR := Quotient(PowerSeriesRing(CoefficientRing(R)),prec); 1505 else 1506 QR := Quotient(R,pi^prec); 1507 fi; 1508 else 1509 QR := R; 1510 fi; 1511 PQR:= PolynomialAlgebra(QR); 1512 g := Coerce(PQR,g); 1513 uni := LeadingCoefficient(g); 1514 #//print uni; 1515 g := Div(g,uni); 1516 if IsMonic(g)=FALSE then return Error("g not monic"); fi; 1517 return [true, g, uni, multval, exp]; 1518 end; 1519 1520 #// make the reverse changes made to f by make_monic_and_integral to h 1521 #// actually only the exp part, the rest goes into the scalar returned 1522 #// with the polynomial factorization. 1523 undo_monic_and_integral := function(f,exp,h) 1524 local g,PK,K,pi,g; 1525 g := h;#//*uni; 1526 1527 PK := Parent(f); 1528 K := CoefficientRing(PK); 1529 pi := UniformizingElement(K); 1530 1531 g := Coerce(PK,g); 1532 g := Evaluate(g,pi^exp*PK.1); 1533 1534 return g; 1535 end; 1536 1537 1538 1539 #///////////////////////////////////////////////////////////////////// 1540 1541 #{For a polynomial f over a local ring or field return a precision at which 1542 # the factorization of f as given by Factorization(f) will be Hensel liftable 1543 # to the correct factorization.} 1544 1545 1546 SuggestedPrecision:= function(f) 1547 local max,seq,vallist,tmp_R,tmpi,mo,ok,h,uni,multval,exp,mmai_prec,PK,PR,P,valList,newloss,loss,gcd, 1548 old_P,dh,old_gcd,pgcd,done,good,hh,val; 1549 #// the trivial cases 1550 if Degree(f) = 0 then 1551 max := Maximum([2,Valuation(TrailingCoefficient(f))+1]); 1552 return max; 1553 elif Degree(f) = -1 then 1554 return 1; 1555 fi; 1556 1557 #// monic and integral 1558 seq:= ElementToSequence(f); 1559 vallist:=[1..Length(seq)]; 1560 Apply_(vallist,x->Valuation(seq[x])); 1561 if ((IsUnit(LeadingCoefficient(f)) and 1562 not IsField(CoefficientRing(Parent(f)))) or 1563 IsMonic(f)) and 1564 Minimum(vallist) >= 0 then 1565 #// make sure it is also squarefree 1566 if Valuation(Discriminant(f)) < 1567 Precision(CoefficientRing(Parent(f))) 1568 and Precision(CoefficientRing(Parent(f))) < INFTY then 1569 max := Maximum([2*RamificationDegree(CoefficientRing(Parent(f))), 1570 2*Valuation(Discriminant(f))]); 1571 return max; 1572 fi; 1573 fi; 1574 tmp_R := CoefficientRing(Parent(f)); 1575 tmpi:=UniformizingElement(tmp_R); 1576 if Is(Type(tmp_R),ord^ser) or Is(Type(tmp_R),res^ser) then 1577 f := Coerce(PolynomialAlgebra( Quotient(tmp_R,(2*Precision(tmp_R)))),f); 1578 else 1579 f := Coerce(PolynomialAlgebra( Quotient(tmp_R,tmpi^(2*Precision(tmpi)))) ,f); 1580 fi; 1581 mo:=make_monic_and_integral(f); 1582 ok:= mo[1]; 1583 h:= mo[2]; 1584 uni:=mo[3]; 1585 multval:= mo[4]; 1586 exp:= mo[5]; 1587 mmai_prec := Valuation(LeadingCoefficient(f))+multval+exp*Degree(f); 1588 if not ok then 1589 #//print "increasing prec"; 1590 mmai_prec := multval+exp*Degree(f); 1591 PK := PolynomialAlgebra(ChangePrecision(CoefficientRing(Parent(f)), 1592 mmai_prec)); 1593 mo :=make_monic_and_integral(Coerce(PK,f)); 1594 good:= mo[1]; 1595 h:= mo[2]; 1596 uni:= mo[3]; 1597 multval:= mo[4]; 1598 exp:= mo[5]; 1599 if not good then 1600 return Error ("I cannot suggest a precision."); 1601 fi; 1602 fi; 1603 1604 PR := Parent(h); 1605 P := CoefficientRing(Parent(h)); 1606 1607 #//print h; 1608 1609 # // correct the precision for the derivative 1610 seq := ElementToSequence(h); 1611 valList:=[2..Length(seq)]; 1612 Apply_(valList,x->Valuation(seq[x])+Valuation(Coerce(P,(x-1)))); 1613 max := Maximum(valList); 1614 #//print "precision increase derivative",max-Precision(P); 1615 newloss := Maximum([max-Precision(P),0]); 1616 #//RamificationDegree(P) eq 1 and InertiaDegree(P) eq 1 then 1617 repeat #// loop for the valuation of the discriminant 1618 loss := 0; 1619 gcd := 0; 1620 repeat 1621 old_P := P; 1622 if not Precision(P) =INFTY then 1623 P := ChangePrecision(P,Precision(P)+(newloss)); 1624 fi; 1625 #//print "precision",Precision(P); 1626 PK := PolynomialAlgebra(ChangePrecision(CoefficientRing(Parent(f)), 1627 Precision(P))); 1628 mo:=make_monic_and_integral(Coerce(PK,f)); 1629 good:= mo[1]; 1630 h:= mo[2]; 1631 uni:=mo[3]; 1632 multval:= mo[4]; 1633 exp:= mo[5]; 1634 #//print "new h",h; 1635 if not good then 1636 return Error ("I cannot suggest a precision."); 1637 fi; 1638 h := Coerce(PolynomialAlgebra(P),h); 1639 dh := Derivative(h); 1640 old_gcd := gcd; 1641 #//print "h",h; 1642 #//print "dh",dh; 1643 #vprint RoundFour,4: "R4: gcd"; 1644 pgcd:=PGcd(h,dh); 1645 gcd:= pgcd[1]; newloss := pgcd[2]; 1646 #//print "newloss",newloss; 1647 loss := loss + newloss; 1648 #//print "gcd",gcd, newloss; 1649 #//print "h mod gcd",h mod gcd; 1650 #//print "h div gcd",h div gcd; 1651 until (old_gcd = Coerce(PolynomialAlgebra(old_P),gcd) and h mod gcd = 0); 1652 hh := Div(h,gcd); 1653 #//print hh; 1654 if RamificationDegree(P) = 1 and InertiaDegree(P) = 1 then 1655 val := Valuation(Discriminant(Coerce(PolynomialAlgebra(RationalField()),hh)), 1656 Prime(P)); 1657 done := val.base <> INFTY; 1658 #//print val; 1659 else 1660 val := Valuation(Discriminant(Coerce(PolynomialAlgebra(P),hh))); 1661 done := val < Precision(P); 1662 fi; 1663 #//print "val",2*val,loss,mmai_prec; 1664 until done and newloss < Precision(P); 1665 max := Maximum([2*val+loss,mmai_prec+loss,2,Precision(P)]); 1666 return max; 1667 end; 1668 1669 #{Factorization of a polynomial f over a local ring or a local field into monic irreducible 1670 # factors together with the appropriate scalar and certificates for the irreducibility of the 1671 # factors if requested } 1672 1673 LocalFactorization:= function(f,Certificates,IsSquarefree,Ideals,Extensions) 1674 local with_certificates,mo,ok,h,uni,multval,exp,PQR,QR,sqrfr,is_sq_free,sq_free_part,loss,r4,res1,g, i,hh,rem,zero,done,sequ,fa,valList,a,j,PR,R,pi,min,sca,K,PK,Kpi,cert,poly,copoly,gcd,G,C,D,Phi, traces_Phi,Rho,U,PU,Pi,tmp_U,g2,tmp_PU,fact,T,eltseq,vallist,res2,qr,prec,un,tmpPR,tmpf; 1675 with_certificates := Certificates or Ideals or Extensions; 1676 1677 if Precision(CoefficientRing(Parent(f))) = INFTY and 1678 Extensions then 1679 return Error ("Factorization: 'Extensions' can only be used over finite precision rings/fields"); 1680 fi; 1681 1682 #//print "Factorization",f; 1683 mo:=make_monic_and_integral(f); 1684 ok:= mo[1]; 1685 h:=mo[2]; 1686 uni:=mo[3]; 1687 multval:=mo[4]; 1688 exp := mo[5]; 1689 if not ok then 1690 return Error ("Factorization: Insufficient precision"); 1691 fi; 1692 #//print h; 1693 PQR := Parent(h); 1694 QR := CoefficientRing(PQR); 1695 1696 sqrfr:= is_squarefree(h); 1697 is_sq_free:= sqrfr[1]; 1698 sq_free_part:=sqrfr[2]; 1699 loss:=sqrfr[3]; 1700 if is_sq_free then 1701 1702 r4 := round_four(h,with_certificates); 1703 res1:=[1..Length(r4)]; 1704 Apply_(res1,x->Tuple([undo_monic_and_integral(f, exp, r4[x].factor),1])); 1705 else 1706 if loss>0 and not Is(Type(f),elt-alg^pol/ord^ser) or Is(Type(f),elt-alg^pol/res^ser) then 1707 if Precision(QR) <>INFTY then 1708 tmpPR:= _ChangePrecisionCoeff(PQR,Precision(QR)+loss); 1709 else 1710 tmpPR:= _ChangePrecisionCoeff(PQR,_realPrec+loss); 1711 fi; 1712 tmpf:= Coerce(tmpPR,f); 1713 h:= make_monic_and_integral(tmpf)[2]; 1714 sqrfr:= is_squarefree(h); 1715 sq_free_part:=sqrfr[2]; 1716 fi; 1717 g := sq_free_part; 1718 fa:= 1; 1719 #//"h",h; 1720 #//"g",g; 1721 while IsMonic(g)=false do 1722 # Print("no mo\n"); 1723 #mo:= make_monic_and_integral(g)[2]; 1724 fa:= fa*LeadingCoefficient(g); 1725 mo := g*(IsInvertible(LeadingCoefficient(g)).ext1); 1726 g:= is_squarefree(mo)[2]; 1727 od; 1728 r4 := round_four(g,with_certificates); 1729 if loss>0 and not Is(Type(f),elt-alg^pol/ord^ser) or Is(Type(f),elt-alg^pol/res^ser) then 1730 for i in [1..Length(r4)[1]] do 1731 if Precision(QR)<>INFTY then 1732 r4[1][i].factor:= Coerce(PQR,r4[1][i].factor); 1733 else tmpPR:= _ChangePrecisionCoeff(PQR,_realPrec); 1734 r4[1][i].factor:= Coerce(tmpPR,r4[1][i].factor); 1735 r4[1][i].factor:= Coerce(PQR,r4[1][i].factor); 1736 fi; 1737 od; 1738 fi; 1739 res1:=[1..Length(r4)]; 1740 Apply_(res1,x->Tuple([r4[x].factor,0])); 1741 #//hh := h; 1742 for i in[1..Length(res1)] do 1743 hh := h; 1744 repeat 1745 #//print "hh pre",hh; 1746 hh := Div(hh,res1[i][1]); 1747 rem := hh mod res1[i][1]; 1748 #Print(" = hh mod res ",rem,"\n"); 1749 res1[i][2] := res1[i][2]+1; 1750 #//print "hh post",hh; 1751 zero := (rem = Coerce(Parent(rem),0)); 1752 done := false; 1753 if not zero then 1754 sequ:= ElementToSequence(rem); 1755 valList:=[1..Length(sequ)]; 1756 Apply_(valList,x->Valuation(sequ[x])); 1757 done := Minimum(valList) < Precision(QR)-loss; 1758 fi; 1759 until done; 1760 od; 1761 #//print res; 1762 #//print "exp sum", &+[a[2]*Degree(a[1]) : a in res],"=?=",Degree(f); 1763 # //print res; 1764 a:= res1[1][2]*Degree(res1[1][1]); 1765 for j in [2..Length(res1)] do 1766 a:= a+res1[j][2]*Degree(res1[j][1]); 1767 od; 1768 if a<> Degree(f) then 1769 Print("a=",a); 1770 #//"res",res; 1771 #//"f",f; 1772 return Error ("Factorization: Insufficient precision in factorization of a non-squarefree polynomial."); 1773 fi; 1774 res2:= [1..Length(res1)]; 1775 Apply_(res2,x->Tuple([undo_monic_and_integral(f, exp, res1[x][1]), res1[x][2]])); 1776 res1:= res2; 1777 #//print res; 1778 fi; #// squarefee by gcd 1779 #// end if; // is_squarefree_easy 1780 1781 PR := Parent(f); 1782 R := CoefficientRing(PR); 1783 pi := UniformizingElement(R); 1784 if not IsField(R) then #// distribute the scalar factor 1/pi^multval 1785 for i in [1..Length(res1)] do #// over the factors 1786 eltseq:= ElementToSequence(res1[i][1]); 1787 vallist:=[1..Length(eltseq)]; 1788 Apply_(vallist,x->Valuation(eltseq[x])); 1789 min := Minimum(vallist); 1790 res1[i][1]:=Div( res1[i][1],pi^min); 1791 multval := multval -res1[i][2]*min; 1792 od; 1793 if multval > 0 then 1794 return Error ("Factorization: Something went wrong when trying to factor a non-monic Polynomial."); 1795 fi; 1796 fi; 1797 sca := Coerce(R,uni) * pi^(-multval); 1798 1799 # //========================================== 1800 if with_certificates then 1801 # CERT := recformat<F:IntegerRing(), // inertia deg 1802 # Rho:RngUPolElt, // unramified certificate 1803 # E:IntegerRing(), // ram index 1804 # Pi:RngUPolElt, // ramified certificate 1805 # IdealGen1, // 2 element representation of 1806 # IdealGen2:RngUPolElt,// the corresponding global ideal 1807 # Extension>; // extension of R given by the 1808 # // factor 1809 if not IsField(R) then 1810 #if Is(Type(R),res^pad) then 1811 #K:= pAdicField(Coerce(Z,UniformizingElement(R)),Precision(R)); 1812 #else 1813 K := FieldOfFractions(R); 1814 #fi; 1815 PK := PolynomialAlgebra(K); 1816 else 1817 K := R; 1818 PK := PR; 1819 fi; 1820 1821 1822 Kpi := UniformizingElement(K); 1823 cert:=[1..Length(r4)]; 1824 Apply_(cert,x->rec(F:= r4[x].F, Rho:=Coerce(PK,((r4[x].Gamma))/Kpi^r4[x].Gamma_denval), 1825 E:=r4[x].E, Pi:=Coerce(PK,Coerce(PK,r4[x].Pi)/Kpi^r4[x].Pi_denval))); 1826 1827 #//----------------------------------------- 1828 if Ideals then 1829 for i in [1..Length(res1)] do 1830 cert[i].IdealGen1 := pi; 1831 poly := Coerce(PK,res1[i][1]); 1832 copoly := 1; 1833 for j in [1..Length(res1)] do 1834 if j <> i then 1835 copoly := copoly *Coerce(PK,res1[j][1]); 1836 fi; 1837 od; 1838 gcd:=XGCD(poly,copoly); 1839 G := gcd.base; 1840 C:= gcd.ext1; 1841 D:= gcd.ext2; 1842 cert[i].IdealGen2 := cert[i].Pi*D*copoly+poly*C; 1843 od; 1844 fi; 1845 #//----------------------------------------- 1846 if Extensions then 1847 #vprint RoundFour,1:" -- computing local extensions"; 1848 for i in [1..Length(res1)] do 1849 Phi := r4[i].factor; 1850 #//print "Phi",Phi; 1851 traces_Phi := traces_from_poly(Phi); 1852 if cert[i].F > 1 then 1853 Rho := eltRec(r4[i].Gamma, r4[i].Gamma_denval); 1854 Rho:=elt_invariants(Phi, traces_Phi, Rho); 1855 # //print "Gamma",r4[i]`Gamma, r4[i]`Gamma_denval; 1856 # //print "chi_Rho",chi_Rho; 1857 if not Rho.passes_hensel_test then 1858 return Error ("Factorization: something went wrong with the unramified certificate"); 1859 fi; 1860 if Rho.f <> cert[i].F then return Error("Error in Extensions");fi; 1861 U := UnramifiedExtension(R,Rho.nu); 1862 else 1863 U := R; 1864 fi; 1865 if cert[i].E > 1 then 1866 Pi := eltRec(r4[i].Pi, r4[i].Pi_denval); 1867 Pi:=elt_invariants(Phi, traces_Phi, Pi); 1868 PU := PolynomialAlgebra(U); 1869 #// unfortunately the precison used to factor f must noty be 1870 #// sufficient for factoring chi_Pi 1871 tmp_U := U; 1872 #//tmp_U := ChangePrecision(U, 1873 #// Maximum([2*Valuation(Discriminant(chi_Pi)),Precision(R)])); 1874 tmp_PU := PolynomialAlgebra(tmp_U); 1875 #//vtime RoundFour,4: 1876 fact := LocalFactorization(Coerce(tmp_PU,(Coerce(PU,(Coerce(PR,Pi.chi))))),false,false,false,false); 1877 #//"fact4",fact; 1878 if Degree(fact[1][1][1]) <> Div(Degree(Phi),cert[i].F) then 1879 return Error("Degree"); 1880 fi; 1881 if (2*Valuation(Discriminant(fact[1][1][1]))+cert[i].E)/cert[i].E 1882 > Precision(U) then 1883 return Error ("Factorization: The Precision was not high enough", 1884 "to assure that the ramified extension is unique"); 1885 fi; 1886 T := TotallyRamifiedExtension(U,Coerce(PU,fact[1][1][1])); 1887 else 1888 T := U; 1889 fi; 1890 cert[i].Extension := T; 1891 od; 1892 fi; 1893 return [res1,sca, cert]; 1894 #//=========================================== 1895 else 1896 return [res1,sca]; 1897 fi; 1898 end; 1899 1900 1901 1902 #{True iff f is irreducible} 1903 IsIrreducible:= function(f) 1904 local F; 1905 F := Factorization(f); 1906 return Length(F) = 1; 1907 end; 1908 1909 1910 #{Factorization of a polynomial f over a local ring or a local field into monic irreducible 1911 # factors together with the appropriate scalar and certifcates for the irreducibility of the 1912 # factors if requested } 1913 1914 #Factorization:= function(f,Certificates,IsSquarefree,Ideals,Extensions) 1915 # return LocalFactorization( 1916 # f ,Certificates, IsSquarefree ,Ideals,Extensions); 1917 #end; 1918 1919 InstallMethod( 1920 rec( 1921 name := "Factorization", 1922 kind := "FUNCTION", 1923 sin := [[elt-alg^pol/any^loc,"Phi"]], 1924 sou := [[list]], 1925 opt := [[elt-alg^boo,"Certificates", 1926 "If TRUE two element certificates for the irreducibilty of factors are returned.", rec(Default:=FALSE)], 1927 [elt-alg^boo,"IsSquarefree", 1928 "If TRUE the polynomial is assumed to be squarefree",rec(Default:=FALSE)], 1929 [elt-alg^boo,"Ideals","",rec(Default:=FALSE)], 1930 [elt-alg^boo,"Extensions", 1931 "If TRUE the extensions generated by the irreducible factors are returned", 1932 rec(Default:=FALSE)]], 1933 short := "The factorization of a polynomial over a p-adic field or a "+ 1934 " p-adic ring. The algorithm used is a combinantion of the algorithms by "+ 1935 "Ford-Pauli-Roblot and by Pauli." 1936 ), 1937 function(arg) 1938 local opt; 1939 #Print("arg",arg,"\n"); 1940 opt := Optarg(arg, 1941 ["Certificates", "IsSquarefree", "Ideals", "Extensions"], 1942 [FALSE,FALSE,FALSE,FALSE]); 1943 1944 #Print("opt",opt,"\n"); 1945 #Print("call",arg[1],opt.Certificates, 1946 # opt.IsSquarefree, 1947 # opt.Ideals, 1948 # opt.Extensions,"\n"); 1949 return LocalFactorization(arg[1],opt.("Certificates"), 1950 opt.("IsSquarefree"), 1951 opt.("Ideals"), 1952 opt.("Extensions")); 1953 # return LocalFactorization(arg[1],TRUE,FALSE,FALSE,FALSE); 1954 end 1955 ); 1956 1957 #q3 := pAdicRing(3,30); 1958 #q3x := PolynomialRing(q3); 1959 #f := (q3x.1^3+3)*(q3x.1^3+6); 1960 #LocalFactorization(f,FALSE,FALSE,FALSE,FALSE); 1961 #ReadLib("locFact"); 1962 #Factorization(f); 1963 #q3 := pAdicRing(3,30); 1964 #q3x := PolynomialRing(q3); 1965