Annotation of OpenXM_contrib2/asir2000/lib/fff, Revision 1.1
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/lib/fff,v 1.1.1.1 1999/11/10 08:12:31 noro Exp $ */
! 2: /*
! 3: fff : Univariate factorizer over a finite field.
! 4:
! 5: Revision History:
! 6:
! 7: 99/05/18 noro the first official version
! 8: 99/06/11 noro
! 9: 99/07/27 noro
! 10: */
! 11:
! 12: #include "defs.h"
! 13:
! 14: extern TPMOD,TQMOD$
! 15:
! 16: /*
! 17: Input : a univariate polynomial F
! 18: Output: a list [[F1,M1],[F2,M2],...], where
! 19: Fi is a monic irreducible factor, Mi is its multiplicity.
! 20: The leading coefficient of F is abondoned.
! 21: */
! 22:
! 23: def fctr_ff(F)
! 24: {
! 25: F = simp_ff(F);
! 26: F = F/LCOEF(F);
! 27: L = sqfr_ff(F);
! 28: for ( R = [], T = L; T != []; T = cdr(T) ) {
! 29: S = car(T); A = S[0]; E = S[1];
! 30: B = ddd_ff(A);
! 31: R = append(append_mult_ff(B,E),R);
! 32: }
! 33: return R;
! 34: }
! 35:
! 36: /*
! 37: Input : a list of polynomial L; an integer E
! 38: Output: a list s.t. [[L0,E],[L1,E],...]
! 39: where Li = L[i]/leading coef of L[i]
! 40: */
! 41:
! 42: def append_mult_ff(L,E)
! 43: {
! 44: for ( T = L, R = []; T != []; T = cdr(T) )
! 45: R = cons([car(T)/LCOEF(car(T)),E],R);
! 46: return R;
! 47: }
! 48:
! 49: /*
! 50: Input : a polynomial F
! 51: Output: a list [[F1,M1],[F2,M2],...]
! 52: where Fi is a square free factor,
! 53: Mi is its multiplicity.
! 54: */
! 55:
! 56: def sqfr_ff(F)
! 57: {
! 58: V = var(F);
! 59: F1 = diff(F,V);
! 60: L = [];
! 61: /* F=H*Fq^p => F'=H'*Fq^p => gcd(F,F')=gcd(H,H')*Fq^p */
! 62: if ( F1 != 0 ) {
! 63: F1 = F1/LCOEF(F1);
! 64: F2 = ugcd(F,F1);
! 65: /* FLAT = H/gcd(H,H') : square free part of H */
! 66: FLAT = sdiv(F,F2);
! 67: I = 0;
! 68: /* square free factorization of H */
! 69: while ( deg(FLAT,V) ) {
! 70: while ( 1 ) {
! 71: QR = sqr(F,FLAT);
! 72: if ( !QR[1] ) {
! 73: F = QR[0]; I++;
! 74: } else
! 75: break;
! 76: }
! 77: if ( !deg(F,V) )
! 78: FLAT1 = simp_ff(1);
! 79: else
! 80: FLAT1 = ugcd(F,FLAT);
! 81: G = sdiv(FLAT,FLAT1);
! 82: FLAT = FLAT1;
! 83: L = cons([G,I],L);
! 84: }
! 85: }
! 86: /* now F = Fq^p */
! 87: if ( deg(F,V) ) {
! 88: Char = characteristic_ff();
! 89: T = sqfr_ff(pthroot_p_ff(F));
! 90: for ( R = []; T != []; T = cdr(T) ) {
! 91: H = car(T); R = cons([H[0],Char*H[1]],R);
! 92: }
! 93: } else
! 94: R = [];
! 95: return append(L,R);
! 96: }
! 97:
! 98: /*
! 99: Input : a polynomial F
! 100: Output: F^(1/char)
! 101: */
! 102:
! 103: def pthroot_p_ff(F)
! 104: {
! 105: V = var(F);
! 106: DVR = characteristic_ff();
! 107: PWR = DVR^(extdeg_ff()-1);
! 108: for ( T = F, R = 0; T; ) {
! 109: D1 = deg(T,V); C = coef(T,D1,V); T -= C*V^D1;
! 110: R += C^PWR*V^idiv(D1,DVR);
! 111: }
! 112: return R;
! 113: }
! 114:
! 115: /*
! 116: Input : a polynomial F of degree N
! 117: Output: a list [V^Ord mod F,Tab]
! 118: where V = var(F), Ord = field order
! 119: Tab[i] = V^(i*Ord) mod F (i=0,...,N-1)
! 120: */
! 121:
! 122: def tab_ff(F)
! 123: {
! 124: V = var(F);
! 125: N = deg(F,V);
! 126: F = F/LCOEF(F);
! 127: XP = pwrmod_ff(F);
! 128: R = pwrtab_ff(F,XP);
! 129: return R;
! 130: }
! 131:
! 132: /*
! 133: Input : a square free polynomial F
! 134: Output: a list [F1,F2,...]
! 135: where Fi is an irreducible factor of F.
! 136: */
! 137:
! 138: def ddd_ff(F)
! 139: {
! 140: V = var(F);
! 141: if ( deg(F,V) == 1 )
! 142: return [F];
! 143: TAB = tab_ff(F);
! 144: for ( I = 1, W = V, L = []; 2*I <= deg(F,V); I++ ) {
! 145: lazy_lm(1);
! 146: for ( T = 0, K = 0; K <= deg(W,V); K++ )
! 147: if ( C = coef(W,K,V) )
! 148: T += TAB[K]*C;
! 149: lazy_lm(0);
! 150: W = simp_ff(T);
! 151: GCD = ugcd(F,W-V);
! 152: if ( deg(GCD,V) ) {
! 153: L = append(berlekamp_ff(GCD,I,TAB),L);
! 154: F = sdiv(F,GCD);
! 155: W = urem(W,F);
! 156: }
! 157: }
! 158: if ( deg(F,V) )
! 159: return cons(F,L);
! 160: else
! 161: return L;
! 162: }
! 163:
! 164: /*
! 165: Input : a polynomial
! 166: Output: 1 if F is irreducible
! 167: 0 otherwise
! 168: */
! 169:
! 170: def irredcheck_ff(F)
! 171: {
! 172: V = var(F);
! 173: if ( deg(F,V) <= 1 )
! 174: return 1;
! 175: F1 = diff(F,V);
! 176: if ( !F1 )
! 177: return 0;
! 178: F1 = F1/LCOEF(F1);
! 179: if ( deg(ugcd(F,F1),V) > 0 )
! 180: return 0;
! 181: TAB = tab_ff(F);
! 182: for ( I = 1, W = V, L = []; 2*I <= deg(F,V); I++ ) {
! 183: for ( T = 0, K = 0; K <= deg(W,V); K++ )
! 184: if ( C = coef(W,K,V) )
! 185: T += TAB[K]*C;
! 186: W = T;
! 187: GCD = ugcd(F,W-V);
! 188: if ( deg(GCD,V) )
! 189: return 0;
! 190: }
! 191: return 1;
! 192: }
! 193:
! 194: /*
! 195: Input : a square free (canonical) modular polynomial F
! 196: Output: a list of polynomials [LF,CF,XP] where
! 197: LF=the product of all the linear factors
! 198: CF=F/LF
! 199: XP=x^field_order mod CF
! 200: */
! 201:
! 202: def meq_linear_part_ff(F)
! 203: {
! 204: F = simp_ff(F);
! 205: F = F/LCOEF(F);
! 206: V = var(F);
! 207: if ( deg(F,V) == 1 )
! 208: return [F,1,0];
! 209: T0 = time()[0];
! 210: XP = pwrmod_ff(F);
! 211: GCD = ugcd(F,XP-V);
! 212: if ( deg(GCD,V) ) {
! 213: GCD = GCD/LCOEF(GCD);
! 214: F = sdiv(F,GCD);
! 215: XP = srem(XP,F);
! 216: R = GCD;
! 217: } else
! 218: R = 1;
! 219: TPMOD += time()[0]-T0;
! 220: return [R,F,XP];
! 221: }
! 222:
! 223: /*
! 224: Input : a square free polynomial F s.t.
! 225: all the irreducible factors of F
! 226: has the same degree D.
! 227: Output: D
! 228: */
! 229:
! 230: def meq_ed_ff(F,XP)
! 231: {
! 232: T0 = time()[0];
! 233: F = simp_ff(F);
! 234: F = F/LCOEF(F);
! 235: V = var(F);
! 236:
! 237: TAB = pwrtab_ff(F,XP);
! 238:
! 239: D = deg(F,V);
! 240: for ( I = 1, W = V, L = []; 2*I <= D; I++ ) {
! 241: lazy_lm(1);
! 242: for ( T = 0, K = 0; K <= deg(W,V); K++ )
! 243: if ( C = coef(W,K,V) )
! 244: T += TAB[K]*C;
! 245: lazy_lm(0);
! 246: W = simp_ff(T);
! 247: if ( W == V ) {
! 248: D = I; break;
! 249: }
! 250: }
! 251: TQMOD += time()[0]-T0;
! 252: return D;
! 253: }
! 254:
! 255: /*
! 256: Input : a square free polynomial F
! 257: an integer E
! 258: an array TAB
! 259: where all the irreducible factors of F has degree E
! 260: and TAB[i] = V^(i*Ord) mod F. (V=var(F), Ord=field order)
! 261: Output: a list containing all the irreducible factors of F
! 262: */
! 263:
! 264: def berlekamp_ff(F,E,TAB)
! 265: {
! 266: V = var(F);
! 267: N = deg(F,V);
! 268: Q = newmat(N,N);
! 269: for ( J = 0; J < N; J++ ) {
! 270: T = urem(TAB[J],F);
! 271: for ( I = 0; I < N; I++ ) {
! 272: Q[I][J] = coef(T,I);
! 273: }
! 274: }
! 275: for ( I = 0; I < N; I++ )
! 276: Q[I][I] -= simp_ff(1);
! 277: L = nullspace_ff(Q); MT = L[0]; IND = L[1];
! 278: NF0 = N/E;
! 279: PS = null_to_poly_ff(MT,IND,V);
! 280: R = newvect(NF0); R[0] = F/LCOEF(F);
! 281: for ( I = 1, NF = 1; NF < NF0 && I < NF0; I++ ) {
! 282: PSI = PS[I];
! 283: MP = minipoly_ff(PSI,F);
! 284: ROOT = find_root_ff(MP); NR = length(ROOT);
! 285: for ( J = 0; J < NF; J++ ) {
! 286: if ( deg(R[J],V) == E )
! 287: continue;
! 288: for ( K = 0; K < NR; K++ ) {
! 289: GCD = ugcd(R[J],PSI-ROOT[K]);
! 290: if ( deg(GCD,V) > 0 && deg(GCD,V) < deg(R[J],V) ) {
! 291: Q = sdiv(R[J],GCD);
! 292: R[J] = Q; R[NF++] = GCD;
! 293: }
! 294: }
! 295: }
! 296: }
! 297: return vtol(R);
! 298: }
! 299:
! 300: /*
! 301: Input : a matrix MT
! 302: an array IND
! 303: an indeterminate V
! 304: MT is a matrix after Gaussian elimination
! 305: IND[I] = 0 means that i-th column of MT represents a basis
! 306: element of the null space.
! 307: Output: an array R which contains all the basis element of
! 308: the null space of MT. Here, a basis element E is represented
! 309: as a polynomial P of V s.t. coef(P,i) = E[i].
! 310: */
! 311:
! 312: def null_to_poly_ff(MT,IND,V)
! 313: {
! 314: N = size(MT)[0];
! 315: for ( I = 0, J = 0; I < N; I++ )
! 316: if ( IND[I] )
! 317: J++;
! 318: R = newvect(J);
! 319: for ( I = 0, L = 0; I < N; I++ ) {
! 320: if ( !IND[I] )
! 321: continue;
! 322: for ( J = K = 0, T = 0; J < N; J++ )
! 323: if ( !IND[J] )
! 324: T += MT[K++][I]*V^J;
! 325: else if ( J == I )
! 326: T -= V^I;
! 327: R[L++] = simp_ff(T);
! 328: }
! 329: return R;
! 330: }
! 331:
! 332: /*
! 333: Input : a polynomial P, a polynomial F
! 334: Output: a minimal polynomial MP(V) of P mod F.
! 335: */
! 336:
! 337: def minipoly_ff(P,F)
! 338: {
! 339: V = var(P);
! 340: P0 = P1 = simp_ff(1);
! 341: L = [[P0,P0]];
! 342: while ( 1 ) {
! 343: /* P0 = V^K, P1 = P^K mod F */
! 344: P0 *= V;
! 345: P1 = urem(P*P1,F);
! 346: /*
! 347: NP0 = P0-c1L1_0-c2L2_0-...,
! 348: NP1 is a normal form w.r.t. [L1_1,L2_1,...]
! 349: NP1 = P1-c1L1_1-c2L2_1-...,
! 350: NP0 represents the normal form computation.
! 351: */
! 352: L1 = lnf_ff(P0,P1,L,V); NP0 = L1[0]; NP1 = L1[1];
! 353: if ( !NP1 )
! 354: return NP0;
! 355: else
! 356: L = lnf_insert([NP0,NP1],L,V);
! 357: }
! 358: }
! 359:
! 360: /*
! 361: Input ; a list of polynomials [P0,P1] = [V^K,P^K mod F]
! 362: a sorted list L=[[L1_0,L1_1],[L2_0,L2_1],...]
! 363: of previously computed pairs of polynomials
! 364: an indeterminate V
! 365: Output: a list of polynomials [NP0,NP1]
! 366: where NP1 = P1-c1L1_1-c2L2_1-...,
! 367: NP0 = P0-c1L1_0-c2L2_0-...,
! 368: NP1 is a normal form w.r.t. [L1_1,L2_1,...]
! 369: NP0 represents the normal form computation.
! 370: [L1_1,L_2_1,...] is sorted so that it is a triangular
! 371: linear basis s.t. deg(Li_1) > deg(Lj_1) for i<j.
! 372: */
! 373:
! 374: def lnf_ff(P0,P1,L,V)
! 375: {
! 376: NP0 = P0; NP1 = P1;
! 377: for ( T = L; T != []; T = cdr(T) ) {
! 378: Q = car(T);
! 379: D1 = deg(NP1,V);
! 380: if ( D1 == deg(Q[1],V) ) {
! 381: H = -coef(NP1,D1,V)/coef(Q[1],D1,V);
! 382: NP0 += Q[0]*H;
! 383: NP1 += Q[1]*H;
! 384: }
! 385: }
! 386: return [NP0,NP1];
! 387: }
! 388:
! 389: /*
! 390: Input : a pair of polynomial P=[P0,P1],
! 391: a list L,
! 392: an indeterminate V
! 393: Output: a list L1 s.t. L1 contains P and L
! 394: and L1 is sorted in the decreasing order
! 395: w.r.t. the degree of the second element
! 396: of elements in L1.
! 397: */
! 398:
! 399: def lnf_insert(P,L,V)
! 400: {
! 401: if ( L == [] )
! 402: return [P];
! 403: else {
! 404: P0 = car(L);
! 405: if ( deg(P0[1],V) > deg(P[1],V) )
! 406: return cons(P0,lnf_insert(P,cdr(L),V));
! 407: else
! 408: return cons(P,L);
! 409: }
! 410: }
! 411:
! 412: /*
! 413: Input : a square free polynomial F s.t.
! 414: all the irreducible factors of F
! 415: has the degree E.
! 416: Output: a list containing all the irreducible factors of F
! 417: */
! 418:
! 419: def c_z_ff(F,E)
! 420: {
! 421: Type = field_type_ff();
! 422: if ( Type == 1 || Type == 3 )
! 423: return c_z_lm(F,E);
! 424: else
! 425: return c_z_gf2n(F,E);
! 426: }
! 427:
! 428: /*
! 429: Input : a square free polynomial P s.t. P splits into linear factors
! 430: Output: a list containing all the root of P
! 431: */
! 432:
! 433: def find_root_ff(P)
! 434: {
! 435: V = var(P);
! 436: L = c_z_ff(P,1);
! 437: for ( T = L, U = []; T != []; T = cdr(T) ) {
! 438: S = car(T)/LCOEF(car(T)); U = cons(-coef(S,0,V),U);
! 439: }
! 440: return U;
! 441: }
! 442:
! 443: /*
! 444: Input : a square free polynomial F s.t.
! 445: all the irreducible factors of F
! 446: has the degree E.
! 447: Output: an irreducible factor of F
! 448: */
! 449:
! 450: def c_z_one_ff(F,E)
! 451: {
! 452: Type = field_type_ff();
! 453: if ( Type == 1 || Type == 3 )
! 454: return c_z_one_lm(F,E);
! 455: else
! 456: return c_z_one_gf2n(F,E);
! 457: }
! 458:
! 459: /*
! 460: Input : a square free polynomial P s.t. P splits into linear factors
! 461: Output: a list containing a root of P
! 462: */
! 463:
! 464: def find_one_root_ff(P)
! 465: {
! 466: V = var(P);
! 467: LF = c_z_one_ff(P,1);
! 468: U = -coef(LF/LCOEF(LF),0,V);
! 469: return [U];
! 470: }
! 471:
! 472: /*
! 473: Input : an integer N; an indeterminate V
! 474: Output: a polynomial F s.t. var(F) = V, deg(F) < N
! 475: and its coefs are random numbers in
! 476: the ground field.
! 477: */
! 478:
! 479: def randpoly_ff(N,V)
! 480: {
! 481: for ( I = 0, S = 0; I < N; I++ )
! 482: S += random_ff()*V^I;
! 483: return S;
! 484: }
! 485:
! 486: /*
! 487: Input : an integer N; an indeterminate V
! 488: Output: a monic polynomial F s.t. var(F) = V, deg(F) = N-1
! 489: and its coefs are random numbers in
! 490: the ground field except for the leading term.
! 491: */
! 492:
! 493: def monic_randpoly_ff(N,V)
! 494: {
! 495: for ( I = 0, N1 = N-1, S = 0; I < N1; I++ )
! 496: S += random_ff()*V^I;
! 497: return V^N1+S;
! 498: }
! 499:
! 500: /* GF(p) specific functions */
! 501:
! 502: /*
! 503: Input : a square free polynomial F s.t.
! 504: all the irreducible factors of F
! 505: has the degree E.
! 506: Output: a list containing all the irreducible factors of F
! 507: */
! 508:
! 509: def c_z_lm(F,E)
! 510: {
! 511: V = var(F);
! 512: N = deg(F,V);
! 513: if ( N == E )
! 514: return [F];
! 515: M = field_order_ff();
! 516: K = idiv(N,E);
! 517: L = [F];
! 518: while ( 1 ) {
! 519: W = monic_randpoly_ff(2*E,V);
! 520: T = generic_pwrmod_ff(W,F,idiv(M^E-1,2));
! 521: W = T-1;
! 522: if ( !W )
! 523: continue;
! 524: G = ugcd(F,W);
! 525: if ( deg(G,V) && deg(G,V) < N ) {
! 526: L1 = c_z_lm(G,E);
! 527: L2 = c_z_lm(sdiv(F,G),E);
! 528: return append(L1,L2);
! 529: }
! 530: }
! 531: }
! 532:
! 533: /*
! 534: Input : a square free polynomial F s.t.
! 535: all the irreducible factors of F
! 536: has the degree E.
! 537: Output: an irreducible factor of F
! 538: */
! 539:
! 540: def c_z_one_lm(F,E)
! 541: {
! 542: V = var(F);
! 543: N = deg(F,V);
! 544: if ( N == E )
! 545: return F;
! 546: M = field_order_ff();
! 547: K = idiv(N,E);
! 548: while ( 1 ) {
! 549: W = monic_randpoly_ff(2*E,V);
! 550: T = generic_pwrmod_ff(W,F,idiv(M^E-1,2));
! 551: W = T-1;
! 552: if ( W ) {
! 553: G = ugcd(F,W);
! 554: D = deg(G,V);
! 555: if ( D && D < N ) {
! 556: if ( 2*D <= N ) {
! 557: F1 = G; F2 = sdiv(F,G);
! 558: } else {
! 559: F2 = G; F1 = sdiv(F,G);
! 560: }
! 561: return c_z_one_lm(F1,E);
! 562: }
! 563: }
! 564: }
! 565: }
! 566:
! 567: /* GF(2^n) specific functions */
! 568:
! 569: /*
! 570: Input : a square free polynomial F s.t.
! 571: all the irreducible factors of F
! 572: has the degree E.
! 573: Output: a list containing all the irreducible factors of F
! 574: */
! 575:
! 576: def c_z_gf2n(F,E)
! 577: {
! 578: V = var(F);
! 579: N = deg(F,V);
! 580: if ( N == E )
! 581: return [F];
! 582: K = idiv(N,E);
! 583: L = [F];
! 584: while ( 1 ) {
! 585: W = randpoly_ff(2*E,V);
! 586: T = tracemod_gf2n(W,F,E);
! 587: W = T-1;
! 588: if ( !W )
! 589: continue;
! 590: G = ugcd(F,W);
! 591: if ( deg(G,V) && deg(G,V) < N ) {
! 592: L1 = c_z_gf2n(G,E);
! 593: L2 = c_z_gf2n(sdiv(F,G),E);
! 594: return append(L1,L2);
! 595: }
! 596: }
! 597: }
! 598:
! 599: /*
! 600: Input : a square free polynomial F s.t.
! 601: all the irreducible factors of F
! 602: has the degree E.
! 603: Output: an irreducible factor of F
! 604: */
! 605:
! 606: def c_z_one_gf2n(F,E)
! 607: {
! 608: V = var(F);
! 609: N = deg(F,V);
! 610: if ( N == E )
! 611: return F;
! 612: K = idiv(N,E);
! 613: while ( 1 ) {
! 614: W = randpoly_ff(2*E,V);
! 615: T = tracemod_gf2n(W,F,E);
! 616: W = T-1;
! 617: if ( W ) {
! 618: G = ugcd(F,W);
! 619: D = deg(G,V);
! 620: if ( D && D < N ) {
! 621: if ( 2*D <= N ) {
! 622: F1 = G; F2 = sdiv(F,G);
! 623: } else {
! 624: F2 = G; F1 = sdiv(F,G);
! 625: }
! 626: return c_z_one_gf2n(F1,E);
! 627: }
! 628: }
! 629: }
! 630: }
! 631:
! 632: /*
! 633: Input : an integer D
! 634: Output: an irreducible polynomial F over GF(2)
! 635: of degree D.
! 636: */
! 637:
! 638: def defpoly_mod2(D)
! 639: {
! 640: return gf2ntop(irredpoly_up2(D,0));
! 641: }
! 642:
! 643: def dummy_time() {
! 644: return [0,0,0,0];
! 645: }
! 646: end$
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>