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