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