![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to sp CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / lib|
Return to sp CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / lib |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/lib/sp,v 1.1.1.1 1999/11/10 08:12:30 noro Exp $ */
2: /*
3: sp : functions related to algebraic number fields
4:
5: Revision History:
6:
7: 99/08/24 noro modified for 1999 release version
8: */
9:
10: #include "defs.h"
11:
12: extern ASCENT,GCDTIME,UFTIME,RESTIME,SQTIME,PRINT$
13: extern Ord$
14:
15: def sp(P)
16: {
17: RESTIME=UFTIME=GCDTIME=SQTIME=0;
18: L = flatmf(fctr(P)); X = var(P);
19: AL = []; ADL = [];
20: while ( 1 ) {
21: L = sort_by_deg(L);
22: for ( T = L, H = []; T != []; H = cons(car(T),H), T = cdr(T) )
23: if ( deg(car(T),X) > 1 )
24: break;
25: if ( T == [] ) {
26: if ( dp_gr_print() ) {
27: print(["GCDTIME = ",GCDTIME]);
28: print(["UFTIME = ",UFTIME]);
29: print(["RESTIME = ",RESTIME]);
30: }
31: return [L,ADL];
32: } else {
33: A = newalg(car(T));
34: R = pdiva(car(T),X-A);
35: AL = cons(A,AL);
36: ADL = cons([A,defpoly(A)],ADL);
37: L = aflist(append(H,append([X-A,R],cdr(T))),AL);
38: }
39: }
40: }
41:
42: def aflist(L,AL)
43: {
44: for ( DC = []; L != []; L = cdr(L) ) {
45: T = af_sp(car(L),AL,1);
46: DC = append(DC,T);
47: }
48: return DC;
49: }
50:
51: def sort_by_deg(F)
52: {
53: for ( T = F, S = []; T != []; T = cdr(T) )
54: if ( type(car(T)) != NUM )
55: S = cons(car(T),S);
56: N = length(S); W = newvect(N);
57: for ( I = 0; I < N; I++ )
58: W[I] = S[I];
59: V = var(W[0]);
60: for ( I = 0; I < N; I++ ) {
61: for ( J = I + 1, J0 = I; J < N; J++ )
62: if ( deg(W[J0],V) > deg(W[J],V) )
63: J0 = J;
64: if ( J0 != I ) {
65: T = W[I]; W[I] = W[J0]; W[J0] = T;
66: }
67: }
68: if ( ASCENT )
69: for ( I = N-1, S = []; I >= 0; I-- )
70: S = cons(W[I],S);
71: else
72: for ( I = 0, S = []; I < N; I++ )
73: S = cons(W[I],S);
74: return S;
75: }
76:
77: def flatmf(L) {
78: for ( S = []; L != []; L = cdr(L) )
79: if ( type(F=car(car(L))) != NUM )
80: S = append(S,[F]);
81: return S;
82: }
83:
84: def af(P,AL)
85: {
86: RESTIME=UFTIME=GCDTIME=SQTIME=0;
87: S = reverse(asq(P,AL));
88: for ( L = []; S != []; S = cdr(S) ) {
89: FM = car(S); F = FM[0]; M = FM[1];
90: G = af_sp(F,AL,1);
91: for ( ; G != []; G = cdr(G) )
92: L = cons([car(G),M],L);
93: }
94: if ( dp_gr_print() )
95: print(["GCDTIME = ",GCDTIME,"UFTIME = ",UFTIME,"RESTIME = ",RESTIME,"SQTIME=",SQTIME]);
96: return L;
97: }
98:
99: def af_sp(P,AL,HINT)
100: {
101: if ( !P || type(P) == NUM )
102: return [P];
103: P1 = simpcoef(simpalg(P));
104: return af_spmain(P1,AL,1,HINT,P1,[]);
105: }
106:
107: def af_spmain(P,AL,INIT,HINT,PP,SHIFT)
108: {
109: if ( !P || type(P) == NUM )
110: return [P];
111: P = simpcoef(simpalg(P));
112: if ( DEG(P) == 1 )
113: return [simpalg(P)];
114: if ( AL == [] ) {
115: TTT = time()[0];
116: F = flatmf(ufctrhint_heuristic(P,HINT,PP,SHIFT));
117: UFTIME+=time()[0]-TTT;
118: return F;
119: }
120: A0 = car(AL); P0 = defpoly(A0);
121: V = var(P); V0 = var(P0);
122: P = simpcoef(P);
123: TTT = time()[0];
124: N = simpcoef(sp_norm(A0,V,subst(P,V,V-INIT*A0),AL));
125: RESTIME+=time()[0]-TTT;
126: TTT = time()[0];
127: DCSQ = sortfs(asq(N,AL));
128: SQTIME+=time()[0]-TTT;
129: for ( G = P, A = V+INIT*A0, DCR = []; DCSQ != []; DCSQ = cdr(DCSQ) ) {
130: C = TT(DCSQ); D = TS(DCSQ);
131: if ( !var(C) )
132: continue;
133: if ( D == 1 )
134: DCT = af_spmain(C,cdr(AL),1,HINT*deg(P0,V0),PP,cons([A0,INIT],SHIFT));
135: else
136: DCT = af_spmain(C,cdr(AL),1,1,C,[]);
137: for ( ; DCT != []; DCT = cdr(DCT) ) {
138: if ( !var(car(DCT)) )
139: continue;
140: if ( length(DCSQ) == 1 && length(DCT) == 1 )
141: U = simpcoef(G);
142: else {
143: S = subst(car(DCT),V,A);
144: if ( pra(G,S,AL) )
145: U = cr_gcda(S,G,AL);
146: else
147: U = S;
148: }
149: if ( var(U) == V ) {
150: G = pdiva(G,U);
151: if ( D == 1 )
152: DCR = cons(simpcoef(U),DCR);
153: else {
154: T = af_spmain(U,AL,sp_next(INIT),HINT,PP,SHIFT);
155: DCR = append(DCR,T);
156: }
157: }
158: }
159: }
160: return DCR;
161: }
162:
163: def sp_next(I)
164: {
165: if ( I > 0 )
166: return -I;
167: else
168: return -I+1;
169: }
170:
171: extern USE_RES;
172:
173: def sp_norm(A,V,P,AL)
174: {
175: P = simpcoef(simpalg(P));
176: if (USE_RES)
177: return sp_norm_res(A,V,P,AL);
178: else
179: return sp_norm_ch(A,V,P,AL);
180: }
181:
182: def sp_norm_ch(A,V,P,AL)
183: {
184: Len = length(AL);
185: P0 = defpoly(A); V0 = var(P0);
186: PR = algptorat(P);
187: if ( nmono(P0) == 2 )
188: R = res(V0,PR,P0);
189: else if ( Len == 1 || Len == 3 )
190: R = res_ch1(V0,V,PR,P0);
191: else if ( Len == 2 ) {
192: P1 = defpoly(AL[1]);
193: R = norm_ch1(V0,V,PR,P0,P1);
194: } else
195: R = res(V0,PR,P0);
196: return rattoalgp(R,cdr(AL));
197: }
198:
199: def sp_norm_res(A,V,P,AL)
200: {
201: Len = length(AL);
202: P0 = defpoly(A); V0 = var(P0);
203: PR = algptorat(P);
204: R = res(V0,PR,P0);
205: return rattoalgp(R,cdr(AL));
206: }
207:
208: def simpalg(P) {
209: if ( !P )
210: return 0;
211: else if ( type(P) == NUM )
212: return ntype(P) <= 1 ? P : simpalgn(P);
213: else if ( type(P) == POLY )
214: return simpalgp(P);
215: else if ( type(P) == RAT )
216: return simpalg(nm(P))/simpalg(dn(P));
217: }
218:
219: def simpalgp(P) {
220: for ( V = var(P), I = deg(P,V), T = 0; I >= 0; I-- )
221: if ( C = coef(P,I) )
222: T += simpalg(C)*V^I;
223: return T;
224: }
225:
226: def simpalgn(A) {
227: if ( ntype(A) <= 1 )
228: return A;
229: else if ( type(R=algtorat(A)) == POLY )
230: return simpalgb(A);
231: else
232: return simpalgb(
233: invalgp(simpalgb(rattoalg(dn(R))))
234: *simpalgb(rattoalg(nm(R)))
235: );
236: }
237:
238: def simpalgb(P) {
239: if ( ntype(P) <= 1 )
240: return P;
241: else {
242: A0 = getalg(P);
243: Used = [];
244: while ( A0 != [] ) {
245: S = algtorat(P);
246: for ( A = A0; A != []; A = cdr(A) )
247: S = srem(S,defpoly(car(A)));
248: P = rattoalg(S);
249: Used = append(Used,[car(A0)]);
250: A0 = setminus(getalg(P),Used);
251: }
252: return P;
253: }
254: }
255:
256: def setminus(A,B) {
257: for ( T = reverse(A), R = []; T != []; T = cdr(T) ) {
258: for ( S = B, M = car(T); S != []; S = cdr(S) )
259: if ( car(S) == M )
260: break;
261: if ( S == [] )
262: R = cons(M,R);
263: }
264: return R;
265: }
266:
267: def getalgp(P) {
268: if ( type(P) <= 1 )
269: return getalg(P);
270: else {
271: for ( V = var(P), I = deg(P,V), T = []; I >= 0; I-- )
272: if ( C = coef(P,I) )
273: T = union(T,getalgp(C));
274: return T;
275: }
276: }
277:
278: def union(A,B)
279: {
280: for ( T = B; T != []; T = cdr(T) )
281: A = union1(A,car(T));
282: return A;
283: }
284:
285: def union1(A,E)
286: {
287: if ( A == [] )
288: return [E];
289: else if ( car(A) == E )
290: return A;
291: else
292: return cons(car(A),union1(cdr(A),E));
293: }
294:
295: def invalgp(A)
296: {
297: if ( ntype(A) <= 1 )
298: return 1/A;
299: P0 = defpoly(mainalg(A)); P = algtorat(A);
300: V = var(P0); G1 = P0;
301: G2 = DEG(P)>=DEG(P0)?srem(P,P0):P;
302: for ( H = 1, X = 1, U1 = 0, U2 = 1; deg(G2,V); ) {
303: D = DEG(G1)-DEG(G2); T = LCOEF(G2)^(D+1);
304: L = sqr(G1*T,G2); Q = car(L); R = car(cdr(L));
305: S = U1*T-U2*Q;
306: M = H^D; M1 = M*X;
307: G1 = G2; G2 = sdiv(R,M1);
308: U1 = U2; U2 = sdiv(S,M1);
309: X = LCOEF(G1); H = sdiv(X^D*H,M);
310: }
311: C = invalgp(rattoalg(srem(P*U2,P0)));
312: return C*rattoalg(U2);
313: }
314:
315: def algptorat(P) {
316: if ( type(P) <= 1 )
317: return algtorat(P);
318: else {
319: for ( V = var(P), I = deg(P,V), T = 0; I >= 0; I-- )
320: if ( C = coef(P,I) )
321: T += algptorat(C)*V^I;
322: return T;
323: }
324: }
325:
326: def rattoalgp(P,M) {
327: for ( T = M, S = P; T != []; T = cdr(T) )
328: S = subst(S,algtorat(FIRST(T)),FIRST(T));
329: return S;
330: }
331: def sortfs(L)
332: {
333: #define Factor(a) car(a)
334: #define Mult(a) car(cdr(a))
335: if ( type(TT(L)) == NUM )
336: L = cdr(L);
337: for ( N = 0, T = L; T != []; T = cdr(T), N++ );
338: P = newvect(N); P1 = newvect(N);
339: for ( I = 0, T = L, R = []; T != []; T = cdr(T) )
340: if ( Mult(car(T)) == 1 ) {
341: R = cons(car(T),R); N--;
342: } else {
343: P[I] = car(T); I++;
344: }
345: for ( J = 0, V = var(Factor(P[0])); J < N; J++ ) {
346: for ( K0 = K = J, D = deg(Factor(P[J]),V); K < N; K++ )
347: if ( deg(Factor(P[K]),V) < D ) {
348: K0 = K;
349: D = deg(Factor(P[K]),V);
350: }
351: P1[J] = P[K0];
352: if ( J != K0 )
353: P[K0] = P[J];
354: }
355: for ( I = N - 1; I >= 0; I-- )
356: R = cons(P1[I],R);
357: return R;
358: }
359:
360: def pdiva(P1,P2)
361: {
362: A = union(getalgp(P1),getalgp(P2));
363: P1 = algptorat(P1); P2 = algptorat(P2);
364: return simpalg(rattoalgp(sdiv(P1*LCOEF(P2)^(DEG(P1)-DEG(P2)+1),P2),A));
365: }
366:
367: def pqra(P1,P2)
368: {
369: if ( type(P2) != POLY )
370: return [P1,0];
371: else if ( (type(P1) != POLY) || (deg(P1,var(P1)) < deg(P2,var(P1))) )
372: return [0,P1];
373: else {
374: A = union(getalgp(P1),getalgp(P2));
375: P1 = algptorat(P1); P2 = algptorat(P2);
376: L = sqr(P1*LCOEF(P2)^(DEG(P1)-DEG(P2)+1),P2);
377: return [simpalg(rattoalgp(L[0],A)),simpalg(rattoalgp(L[1],A))];
378: }
379: }
380:
381: def pra(P1,P2,AL)
382: {
383: if ( type(P2) != POLY )
384: return 0;
385: else if ( (type(P1) != POLY) || (deg(P1,var(P1)) < deg(P2,var(P1))) )
386: return P1;
387: else {
388: F1 = algptorat(P1); F2 = algptorat(P2); ML = map(defpoly,AL);
389: B = append(reverse(ML),[F2]);
390: V0 = var(P1);
391: V = cons(V0,map(algtorat,AL));
392: G = srem_by_nf(F1,B,V,2);
393: return simpalg(rattoalgp(G[0]/G[1],AL));
394: }
395: }
396:
397: def sort_alg(VL)
398: {
399: N = length(VL); W = newvect(N,VL);
400: for ( I = 0; I < N; I++ ) {
401: for ( M = I, J = I + 1; J < N; J++ )
402: if ( W[J] > W[M] )
403: M = J;
404: if ( I != M ) {
405: T = W[I]; W[I] = W[M]; W[M] = T;
406: }
407: }
408: for ( I = N-1, L = []; I >= 0; I-- )
409: L = cons(W[I],L);
410: return L;
411: }
412:
413: def asq(P,AL)
414: {
415: P = simpalg(P);
416: if ( type(P) == NUM )
417: return [[1,1]];
418: else if ( getalgp(P) == [] )
419: return sqfr(P);
420: else {
421: V = var(P); N = DEG(P); A = newvect(N+1); B = newvect(N+1);
422: for ( I = 0, F = P; ;I++ ) {
423: if ( type(F) == NUM )
424: break;
425: F1 = diff(F,V);
426: GCD = cr_gcda(F,F1,AL);
427: FLAT = pdiva(F,GCD);
428: if ( type(GCD) == NUM ) {
429: A[I] = F; B[I] = 1;
430: break;
431: }
432: for ( J = 1, F = GCD; ; J++ ) {
433: L = pqra(F,FLAT); Q = L[0]; R = L[1];
434: if ( R )
435: break;
436: else
437: F = Q;
438: }
439: A[I] = FLAT; B[I] = J;
440: }
441: for ( I = 0, J = 0, L = []; A[I]; I++ ) {
442: J += B[I];
443: if ( A[I+1] )
444: C = pdiva(A[I],A[I+1]);
445: else
446: C = A[I];
447: L = cons([C,J],L);
448: }
449: return L;
450: }
451: }
452:
453: def ufctrhint1(P,HINT)
454: {
455: if ( deg(P,var(P)) == 168 ) {
456: SQ = sqfr(P);
457: if ( length(SQ) == 2 && SQ[1][1] == 1 )
458: return [[1,1],[P,1]];
459: else
460: return ufctrhint(P,HINT);
461: } else
462: return ufctrhint(P,HINT);
463: }
464:
465: def simpcoef(P) {
466: return rattoalgp(ptozp(algptorat(P)),getalgp(P));
467: }
468:
469: def ufctrhint_heuristic(P,HINT,PP,SHIFT) {
470: V = var(P); D = deg(P,V);
471: if ( D == HINT )
472: return [[P,1]];
473: for ( S = 0, L = SHIFT, AL = [], K = 1; L != []; L = cdr(L) ) {
474: A = car(L)[0]; S += A*car(L)[1]; AL = cons(A,AL);
475: K *= deg(defpoly(A),algtorat(A));
476: }
477: PPP = simpcoef(simpalg(subst(PP,V,V-S)));
478: for ( T = P-coef(P,D)*V^D, G = D; T; T -= coef(T,DT)*V^DT )
479: G = igcd(G,DT=deg(T,V));
480: if ( G == 1 ) {
481: if ( K*deg(PPP,V) != deg(P,V) )
482: PPP = cr_gcda(PPP,P,AL);
483: return ufctrhint2(P,HINT,PPP,AL);
484: } else {
485: for ( S = 0, T = P; T; T -= coef(T,DT)*V^DT ) {
486: DT = deg(T,V);
487: S += coef(T,DT)*V^(DT/G);
488: }
489: L = fctr(S);
490: for ( DC = [car(L)], L = cdr(L); L != []; L = cdr(L) ) {
491: H = subst(car(car(L)),V,V^G);
492: HH = cr_gcda(PPP,H,AL);
493: T = ufctrhint2(H,HINT,HH,AL);
494: DC = append(DC,T);
495: }
496: return DC;
497: }
498: }
499:
500: def ufctrhint2(P,HINT,PP,AL)
501: {
502: if ( deg(P,var(P)) == HINT )
503: return [[P,1]];
504: if ( AL == [] )
505: return ufctrhint(P,HINT);
506: L = resfctr(algptorat(PP),map(defpoly,AL),map(algtorat,AL),P);
507: for ( T = reverse(L[1]), DL = []; T != []; T = cdr(T) )
508: DL = cons(deg(car(car(T)),a_),DL);
509: return resfmain(P,L[2],L[0],DL);
510: }
511:
512: def res_det(V,P1,P2)
513: {
514: D1 = deg(P1,V); D2 = deg(P2,V);
515: M = newmat(D1+D2,D1+D2);
516: for ( J = 0; J <= D2; J++ )
517: M[0][J] = coef(P2,D2-J,V);
518: for ( I = 1; I < D1; I++ )
519: for ( J = 0; J <= D2; J++ )
520: M[I][I+J] = M[0][J];
521: for ( J = 0; J <= D1; J++ )
522: M[D1][J] = coef(P1,D1-J,V);
523: for ( I = 1; I < D2; I++ )
524: for ( J = 0; J <= D1; J++ )
525: M[D1+I][I+J] = M[D1][J];
526: return det(M);
527: }
528:
529: def norm_ch1(V0,VM,P,P0,PR) {
530: D = deg(P,V0); D0 = deg(P0,V0); DM = deg(P,VM); N = DM*D0;
531: X = newvect(N+1); V = newvect(N+1); U = newvect(N+1);
532: Min = -idiv(N,2);
533: C = coef(P,D,V0);
534: for ( I = J = 0; I <= N; J++ ) {
535: if ( PRINT )
536: print([J,N]);
537: T=J+Min;
538: if ( subst(C,VM,T) ) {
539: U[I] = srem(res(V0,subst(P,VM,T),P0),PR);
540: X[I++] = T;
541: }
542: }
543: for ( I = 1, M = 1, S = V[0] = U[0]; I <= N; I++ ) {
544: for ( J = 0, T = U[I]; J < I; J++ )
545: T = sdiv(T-V[J],X[I]-X[J]);
546: V[I] = T;
547: M *= (VM-X[I-1]);
548: S += T*M;
549: }
550: return S;
551: }
552:
553: def norm_ch2(V0,VM,P,P0,PR) {
554: D = deg(P,V0); D0 = deg(P0,V0); DM = deg(P,VM); N = DM*D0;
555: X = newvect(N+1); V = newvect(N+1); U = newvect(N+1);
556: Min = -idiv(N,2);
557: C = coef(P,D,V0);
558: for ( I = J = 0; I <= N; J++ ) {
559: T=J+Min;
560: if ( subst(C,VM,T) ) {
561: U[I] = srem(res_det(V0,subst(P,VM,T),P0),PR);
562: X[I++] = T;
563: }
564: }
565: for ( I = 1, M = 1, S = V[0] = U[0]; I <= N; I++ ) {
566: for ( J = 0, T = U[I]; J < I; J++ )
567: T = sdiv(T-V[J],X[I]-X[J]);
568: V[I] = T;
569: M *= (VM-X[I-1]);
570: S += T*M;
571: }
572: return S;
573: }
574:
575: def res_ch1(V0,VM,P,P0) {
576: D = deg(P,V0); D0 = deg(P0,V0); N = deg(P,VM)*D0+deg(P0,VM)*D;
577: X = newvect(N+1); V = newvect(N+1); U = newvect(N+1);
578: Min = -idiv(N,2);
579: C = coef(P,D,V0); C0 = coef(P0,D0,V0);
580: for ( I = J = 0; I <= N; J++ ) {
581: if ( PRINT )
582: print([J,N]);
583: T=J+Min;
584: if ( subst(C,VM,T) && subst(C0,VM,T) ) {
585: U[I] = res(V0,subst(P,VM,T),subst(P0,VM,T));
586: X[I++] = T;
587: }
588: }
589: for ( I = 1, M = 1, S = V[0] = U[0]; I <= N; I++ ) {
590: for ( J = 0, T = U[I]; J < I; J++ )
591: T = sdiv(T-V[J],X[I]-X[J]);
592: V[I] = T;
593: M *= (VM-X[I-1]);
594: S += T*M;
595: }
596: return S;
597: }
598:
599: def res_ch(V0,VM,P,P0) {
600: D = deg(P,V0); D0 = deg(P0,V0); N = deg(P,VM)*D0+deg(P0,VM)*D;
601: X = newvect(N+1); V = newvect(N+1); U = newvect(N+1);
602: Min = -idiv(N,2);
603: C = coef(P,D,V0); C0 = coef(P0,D0,V0);
604: for ( I = J = 0; I <= N; J++ ) {
605: T=J+Min;
606: if ( subst(C,VM,T) && subst(C0,VM,T) ) {
607: U[I] = res_det(V0,subst(P,VM,T),subst(P0,VM,T));
608: X[I++] = T;
609: }
610: }
611: for ( I = 1, M = 1, S = V[0] = U[0]; I <= N; I++ ) {
612: for ( J = 0, T = U[I]; J < I; J++ )
613: T = sdiv(T-V[J],X[I]-X[J]);
614: V[I] = T;
615: M *= (VM-X[I-1]);
616: S += T*M;
617: }
618: return S;
619: }
620:
621: def norm_ch2_lag(V,VM,P,P0,PR) {
622: D0 = deg(P0,V); DM = deg(P,VM); N = DM*D0;
623: Min = -idiv(N,2);
624: for ( A = 1, I = 0; I <= N; I++ )
625: A *= (VM-I-Min);
626: for ( I = 0, S = 0; I <= N; I++ ) {
627: R = res_det(V,subst(P,VM,I+Min),P0);
628: R = srem(R,PR);
629: T = sdiv(A,VM-I-Min);
630: S += R*T/subst(T,VM,I+Min);
631: }
632: return S;
633: }
634:
635: def norm_ch_lag(V,VM,P,P0) {
636: D0 = deg(P0,V); DM = deg(P,VM); N = DM*D0;
637: Min = -idiv(N,2);
638: for ( A = 1, I = 0; I <= N; I++ )
639: A *= (VM-I-Min);
640: for ( I = 0, S = 0; I <= N; I++ ) {
641: R = res_det(V,subst(P,VM,I+Min),P0);
642: T = sdiv(A,VM-I-Min);
643: S += R*T/subst(T,VM,I+Min);
644: }
645: return S;
646: }
647:
648: def cr_gcda(P1,P2,EXT)
649: {
650: if ( !(V = var(P1)) || !var(P2) )
651: return 1;
652: AL = union(getalgp(P1),getalgp(P2));
653: if ( AL == [] )
654: return gcd(P1,P2);
655: T = newvect(length(EXT));
656: for ( TAL = AL; TAL != []; TAL = cdr(TAL) ) {
657: A = getalg(car(TAL));
658: for ( TA = A; TA != []; TA = cdr(TA) ) {
659: B = car(TA);
660: for ( TEXT = EXT, I = 0; TEXT != []; TEXT = cdr(TEXT), I++ )
661: if ( car(TEXT) == B )
662: T[I] = B;
663: }
664: }
665: for ( I = length(EXT)-1, S = []; I >= 0; I-- )
666: if ( T[I] )
667: S = cons(T[I],S);
668: EXT = S;
669: NEXT = length(EXT);
670: if ( deg(P1,V) < deg(P2,V) ) {
671: T = P1; P1 = P2; P2 = T;
672: }
673: G1 = ptozp(algptorat(P1)); G2 = ptozp(algptorat(P2));
674: for ( ML = VL = [], T = reverse(EXT); T != []; T = cdr(T) ) {
675: ML = cons(defpoly(car(T)),ML);
676: VL = cons(algptorat(car(T)),VL);
677: }
678: DL = [coef(G1,deg(G1,V),V),coef(G2,deg(G2,V),V)];
679: for ( T = EXT; T != []; T = cdr(T) ) {
680: DL = cons(discr(sp_minipoly(car(T),EXT)),DL);
681: C = LCOEF(defpoly(car(T)));
682: if ( C != 1 && C != -1 )
683: DL = cons(C,DL);
684: }
685: TIME = time()[0];
686: for ( D = deg(P1,V)+1, I = 0; ; I++ ) {
687: MOD = lprime(I);
688: for ( J = 0; J < length(DL); J++ )
689: if ( !(DL[J] % MOD) )
690: break;
691: if ( J != length(DL) )
692: continue;
693: Ord = 2; NOSUGAR = 1;
694: T = ag_mod(G1 % MOD,G2 % MOD,ML,VL,MOD);
695: if ( dp_gr_print() )
696: print(".");
697: if ( !T )
698: continue;
699: T = (T*inv(coef(T,deg(T,V),V),MOD))%MOD;
700: if ( deg(T,V) > D )
701: continue;
702: else if ( deg(T,V) < D ) {
703: IMAGE = T; M = MOD; D = deg(T,V);
704: } else {
705: L = cr(IMAGE,M,T,MOD); IMAGE = L[0]; M = L[1];
706: }
707: F = intptoratp(IMAGE,M,calcb(M));
708: if ( F != [] ) {
709: F = ptozp(F);
710: DIV = rattoalgp(F,EXT);
711: if ( type(DIV) == 1 )
712: return 1;
713: /*
714: if ( srem_simp(G1,F,V,ML) )
715: continue;
716: if ( srem_simp(G2,F,V,ML) )
717: continue;
718: */
719: if ( srem_by_nf(G1,reverse(cons(F,ML)),cons(V,VL),2)[0] )
720: continue;
721: if ( srem_by_nf(G2,reverse(cons(F,ML)),cons(V,VL),2)[0] )
722: continue;
723: TIME = time()[0]-TIME;
724: if ( dp_gr_print() )
725: print([TIME]);
726: GCDTIME += TIME;
727: return DIV;
728: }
729: }
730: }
731:
732: def srem_simp(F1,F2,V,D)
733: {
734: D2 = deg(F2,V); C = coef(F2,D2);
735: while ( (D1 = deg(F1,V)) >= D2 ) {
736: F1 -= coef(F1,D1)/C*V^(D1-D2)*F2;
737: F1 = simp_by_dp(F1,D);
738: }
739: return F1;
740: }
741:
742: def member(E,L)
743: {
744: for ( ; L != []; L = cdr(L) )
745: if ( E == car(L) )
746: return 1;
747: return 0;
748: }
749:
750: def getallalg(A)
751: {
752: T = cdr(vars(defpoly(A)));
753: if ( T == [] )
754: return [A];
755: else {
756: for ( S = [A]; T != []; T = cdr(T) )
757: S = union(S,getallalg(rattoalg(car(T))));
758: return S;
759: }
760: }
761:
762: def discr(P) {
763: V = var(P);
764: return res(V,P,diff(P,V));
765: }
766:
767: def sp_minipoly(A,EXT)
768: {
769: while ( car(EXT) != A )
770: EXT = cdr(EXT);
771: for ( M = x-A; EXT != []; EXT = cdr(EXT) )
772: M = sp_norm(car(EXT),x,M,EXT);
773: F = sqfr(M);
774: return F[1][0];
775: }
776:
777: def cr(F1,M1,F2,M2)
778: {
779: K = inv(M1 % M2,M2);
780: M3 = M1*M2;
781: F3 = (F1 + (F2-(F1%M2))*K*M1) % M3;
782: return [F3,M3];
783: }
784:
785: #define ABS(a) ((a)>=0?(a):(-a))
786:
787: #if 0
788: def calcb(M) {
789: setprec(800);
790: return pari(floor,eval((M/2)^(1/2)));
791: }
792: #endif
793:
794: def calcb(M) {
795: N = 2*M;
796: T = sp_sqrt(N);
797: if ( T^2 <= N && N < (T+1)^2 )
798: return idiv(T,2);
799: else
800: error("afo");
801: }
802:
803: def sp_sqrt(A) {
804: for ( J = 0, T = A; T >= 2^27; J++ ) {
805: T = idiv(T,2^27)+1;
806: }
807: for ( I = 0; T >= 2; I++ ) {
808: S = idiv(T,2);
809: if ( T = S+S )
810: T = S;
811: else
812: T = S+1;
813: }
814: X = (2^27)^idiv(J,2)*2^idiv(I,2);
815: while ( 1 ) {
816: if ( (Y=X^2) < A )
817: X += X;
818: else if ( Y == A )
819: return X;
820: else
821: break;
822: }
823: while ( 1 )
824: if ( (Y = X^2) <= A )
825: return X;
826: else
827: X = idiv(A + Y,2*X);
828: }
829:
830: def intptoratp(P,M,B) {
831: if ( type(P) == 1 ) {
832: L = inttorat(P,M,B);
833: if ( L == 0 )
834: return [];
835: else
836: return L[0]/L[1];
837: } else {
838: V = var(P);
839: S = 0;
840: while ( P ) {
841: D = deg(P,V);
842: C = coef(P,D,V);
843: T = intptoratp(C,M,B);
844: if ( T == [] )
845: return [];
846: S += T*V^D;
847: P -= C*V^D;
848: }
849: return S;
850: }
851: }
852:
853: def ltoalg(L) {
854: F = L[0]; V = reverse(L[1]);
855: N = length(V)-1;
856: for ( I = 0, G = F; I < N; I++ ) {
857: D = car(G);
858: A = newalg(D); V = var(D);
859: for ( G = reverse(cdr(G)), T = []; G != []; G = cdr(G) )
860: T = cons(subst(car(G),V,A),T);
861: G = T;
862: }
863: return G;
864: }
865:
866: /*
867: def ag_mod(F1,F2,D,MOD)
868: {
869: if ( length(D) == 1 )
870: return ag_mod_single(F1,F2,D,MOD);
871: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
872: if ( D1 < D2 ) {
873: T = F1; F1 = F2; F2 = T;
874: T = D1; D1 = D2; D2 = T;
875: }
876: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
877: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
878: for ( T = reverse(D), S = []; T != []; T = cdr(T) ) {
879: U = car(T);
880: S = cons((inv(LCOEF(U),MOD)*U) % MOD,S);
881: }
882: D = S;
883: while ( 1 ) {
884: F = srem_simp_mod(F1,F2,V,D,MOD);
885: if ( !F )
886: return F2;
887: if ( !deg(F,V) )
888: return 1;
889: C = LCOEF(F);
890: INV = inverse_by_gr_mod(C,D,MOD);
891: if ( !INV )
892: return 0;
893: F = simp_by_dp_mod(F*INV,D,MOD);
894: F = (inv(LCOEF(F),MOD)*F) % MOD;
895: F1 = F2; F2 = F;
896: }
897: }
898: */
899:
900: def ag_mod(F1,F2,D,VL,MOD)
901: {
902: if ( length(D) == 1 )
903: return ag_mod_single6(F1,F2,D,MOD);
904: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
905: if ( D1 < D2 ) {
906: T = F1; F1 = F2; F2 = T;
907: T = D1; D1 = D2; D2 = T;
908: }
909: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
910: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
911: for ( T = reverse(D), S = []; T != []; T = cdr(T) ) {
912: U = car(T);
913: S = cons((inv(LCOEF(U),MOD)*U) % MOD,S);
914: }
915: D = S;
916: VL = cons(V,VL); B = append([F1,F2],D); N = length(VL);
917: while ( 1 ) {
918: FLAGS = dp_gr_flags(); dp_gr_flags(["Reverse",1,"NoSugar",1]);
919: G = dp_gr_mod_main(B,VL,0,MOD,Ord);
920: dp_gr_flags(FLAGS);
921: if ( length(G) == 1 )
922: return 1;
923: if ( length(G) == N ) {
924: for ( T = G; T != []; T = cdr(T) )
925: if ( member(V,vars(car(T))) )
926: return car(T);
927: }
928: }
929: }
930:
931: def srem_simp_mod(F1,F2,V,D,MOD)
932: {
933: D2 = deg(F2,V); C = coef(F2,D2);
934: while ( (D1 = deg(F1,V)) >= D2 ) {
935: F1 -= coef(F1,D1)/C*V^(D1-D2)*F2;
936: F1 = simp_by_dp_mod(F1,D,MOD);
937: }
938: return F1;
939: }
940:
941: def ag_mod_single(F1,F2,D,MOD)
942: {
943: TD = TI = TM = 0;
944: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
945: if ( D1 < D2 ) {
946: T = F1; F1 = F2; F2 = T;
947: T = D1; D1 = D2; D2 = T;
948: }
949: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
950: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
951: D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
952: FLAGS = dp_gr_flags(); dp_gr_flags(["Reverse",1,"NoSugar",1]);
953: G = dp_gr_mod_main([F1,F2,D],[V,var(D)],0,MOD,2);
954: dp_gr_flags(FLAGS);
955: if ( length(G) == 1 )
956: return 1;
957: if ( length(G) != 2 )
958: return 0;
959: if ( vars(G[0]) == [var(D)] )
960: return G[1];
961: else
962: return G[0];
963: }
964:
965: def ag_mod_single2(F1,F2,D,MOD)
966: {
967: TD = TI = TM = 0;
968: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
969: if ( D1 < D2 ) {
970: T = F1; F1 = F2; F2 = T;
971: T = D1; D1 = D2; D2 = T;
972: }
973: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
974: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
975: D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
976: while ( 1 ) {
977: T0 = time()[0];
978: F = srem((srem(F1,F2) % MOD),D) % MOD;
979: TD += time()[0] - T0;
980: if ( !F ) {
981: if ( dp_gr_print() )
982: print(["TD",TD,"TM",TM,"TI",TI]);
983: return F2;
984: }
985: if ( !deg(F,V) ) {
986: if ( dp_gr_print() )
987: print(["TD",TD,"TM",TM,"TI",TI]);
988: return 1;
989: }
990: C = LCOEF(F);
991: T0 = time()[0];
992: INV = inva_mod(D,MOD,C);
993: TI += time()[0] - T0;
994: if ( !INV )
995: return 0;
996: T0 = time()[0];
997: F = remc_mod((INV*F) % MOD,D,MOD);
998: TM += time()[0] - T0;
999: F1 = F2; F2 = F;
1000: }
1001: }
1002:
1003: def ag_mod_single3(F1,F2,D,MOD)
1004: {
1005: TD = TI = TM = 0;
1006: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
1007: if ( D1 < D2 ) {
1008: T = F1; F1 = F2; F2 = T;
1009: T = D1; D1 = D2; D2 = T;
1010: }
1011: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
1012: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
1013: D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
1014: while ( 1 ) {
1015: if ( !D2 )
1016: return 1;
1017: while ( D1 >= D2 ) {
1018: F = srem((coef(F2,D2,V)*F1-coef(F1,D1,V)*F2*V^(D1-D2))%MOD,D)%MOD;
1019: F1 = F; D1 = deg(F1,V);
1020: }
1021: if ( !F1 ) {
1022: INV = inva_mod(D,MOD,coef(F2,D2,V));
1023: if ( dp_gr_print() )
1024: print(".");
1025: return srem((INV*F2) % MOD,D)%MOD;
1026: } else {
1027: T = F1; F1 = F2; F2 = T;
1028: T = D1; D1 = D2; D2 = T;
1029: }
1030: }
1031: }
1032:
1033: def ag_mod_single4(F1,F2,D,MOD)
1034: {
1035: if ( !F1 )
1036: return F2;
1037: if ( !F2 )
1038: return F1;
1039: TD = TI = TM = TR = 0;
1040: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
1041: if ( D1 < D2 ) {
1042: T = F1; F1 = F2; F2 = T;
1043: T = D1; D1 = D2; D2 = T;
1044: }
1045: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
1046: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
1047: D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
1048: while ( 1 ) {
1049: T0 = time()[0]; R = srem(F1,F2); TR += time()[0] - T0;
1050: T0 = time()[0]; F = srem(R % MOD,D) % MOD; TD += time()[0] - T0;
1051: if ( !F ) {
1052: if ( dp_gr_print() )
1053: print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
1054: return F2;
1055: }
1056: if ( !deg(F,V) ) {
1057: if ( dp_gr_print() )
1058: print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
1059: return 1;
1060: }
1061: C = LCOEF(F);
1062: T0 = time()[0]; INV = inva_mod(D,MOD,C); TI += time()[0] - T0;
1063: if ( !INV )
1064: return 0;
1065: T0 = time()[0]; F = srem((INV*F) % MOD,D) % MOD; TM += time()[0] - T0;
1066: F1 = F2; F2 = F;
1067: }
1068: }
1069:
1070: def ag_mod_single5(F1,F2,D,MOD)
1071: {
1072: TD = TI = TM = TR = 0;
1073: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
1074: if ( D1 < D2 ) {
1075: T = F1; F1 = F2; F2 = T;
1076: T = D1; D1 = D2; D2 = T;
1077: }
1078: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
1079: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
1080: D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
1081: while ( 1 ) {
1082: T0 = time()[0];
1083: D1 = deg(F1,V); D2 = deg(F2,V); F = F1;
1084: while ( D1 >= D2 ) {
1085: R = (F-coef(F,D1,V)*F2*V^(D1-D2))%MOD;
1086: D1 = deg(R,V); HC = coef(R,D1,V);
1087: F = (R - HC*V^D1) + (srem(HC,D)%MOD)*V^D1;
1088: }
1089: TR += time()[0] - T0;
1090: T0 = time()[0]; F = srem(R % MOD,D) % MOD; TD += time()[0] - T0;
1091: if ( !F ) {
1092: if ( dp_gr_print() )
1093: print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
1094: return F2;
1095: }
1096: if ( !deg(F,V) ) {
1097: if ( dp_gr_print() )
1098: print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
1099: return 1;
1100: }
1101: C = LCOEF(F);
1102: T0 = time()[0]; INV = inva_mod(D,MOD,C); TI += time()[0] - T0;
1103: if ( !INV )
1104: return 0;
1105: T0 = time()[0]; F = srem((INV*F) % MOD,D) % MOD; TM += time()[0] - T0;
1106: F1 = F2; F2 = F;
1107: }
1108: }
1109:
1110: def ag_mod_single6(F1,F2,D,MOD)
1111: {
1112: TD = TI = TM = TR = 0;
1113: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
1114: if ( D1 < D2 ) {
1115: T = F1; F1 = F2; F2 = T;
1116: T = D1; D1 = D2; D2 = T;
1117: }
1118: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
1119: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
1120: D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
1121: while ( 1 ) {
1122: T0 = time()[0];
1123: D1 = deg(F1,V); D2 = deg(F2,V); F = F1;
1124: while ( D1 >= D2 ) {
1125: R = (F-coef(F,D1,V)*F2*V^(D1-D2))%MOD;
1126: D1 = deg(R,V); HC = coef(R,D1,V);
1127: /* F = (R - HC*V^D1) + (srem_mod(HC,D,MOD))*V^D1; */
1128: F = remc_mod(R,D,MOD);
1129: }
1130: TR += time()[0] - T0;
1131: T0 = time()[0]; F = remc_mod(R%MOD,D,MOD); TD += time()[0] - T0;
1132: if ( !F ) {
1133: if ( dp_gr_print() )
1134: print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
1135: return F2;
1136: }
1137: if ( !deg(F,V) ) {
1138: if ( dp_gr_print() )
1139: print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
1140: return 1;
1141: }
1142: C = LCOEF(F);
1143: T0 = time()[0]; INV = inva_mod(D,MOD,C); TI += time()[0] - T0;
1144: if ( !INV )
1145: return 0;
1146: T0 = time()[0]; F = remc_mod((INV*F)%MOD,D,MOD); TM += time()[0] - T0;
1147: F1 = F2; F2 = F;
1148: }
1149: }
1150:
1151: def inverse_by_gr_mod(C,D,MOD)
1152: {
1153: Ord = 2;
1154: dp_gr_flags(["NoSugar",1]);
1155: G = dp_gr_mod_main(cons(x*C-1,D),cons(x,vars(D)),0,MOD,Ord);
1156: dp_gr_flags(["NoSugar",0]);
1157: if ( length(G) == 1 )
1158: return 1;
1159: else if ( length(G) == length(D)+1 ) {
1160: for ( T = G; T != []; T = cdr(T) )
1161: if ( member(x,vars(car(G))) )
1162: break;
1163: T = car(G);
1164: if ( type(coef(T,1,x)) != NUM )
1165: return 0;
1166: else
1167: return coef(T,0,x);
1168: } else
1169: return 0;
1170: }
1171:
1172: def simp_by_dp(F,D)
1173: {
1174: for ( T = D; T != []; T = cdr(T) )
1175: F = srem(F,car(T));
1176: return F;
1177: }
1178:
1179: def simp_by_dp_mod(F,D,MOD)
1180: {
1181: F %= MOD;
1182: for ( T = D; T != []; T = cdr(T) )
1183: F = srem(F,car(T)) % MOD;
1184: return F;
1185: }
1186:
1187: def remc_mod(P,D,M)
1188: {
1189: V = var(P);
1190: if ( !V || V == var(D) )
1191: return srem_mod(P,D,M);
1192: for ( I = deg(P,V), S = 0; I >= 0; I-- )
1193: if ( C = coef(P,I,V) )
1194: S += srem_mod(C,D,M)*V^I;
1195: return S;
1196: }
1197:
1198: def rem_mod(C,D,M)
1199: {
1200: V = var(D);
1201: D2 = deg(D,V);
1202: while ( (D1 = deg(C,V)) >= D2 ) {
1203: C -= (D*V^(D1-D2)*coef(C,D1,V))%M;
1204: C %= M;
1205: }
1206: return C;
1207: }
1208:
1209: def resfctr(F,L,V,N)
1210: {
1211: N = ptozp(N);
1212: V0 = var(N);
1213: DN = diff(N,V0);
1214: for ( I = 0, J = 2, Len = deg(N,V0)+1; I < 5; J++ ) {
1215: M = prime(J);
1216: G = gcd(N,DN,M);
1217: if ( !deg(G,V0) ) {
1218: I++;
1219: T = nfctr_mod(N,M);
1220: if ( T < Len ) {
1221: Len = T; M0 = M;
1222: }
1223: }
1224: }
1225: S = spm(L,V,M0);
1226: T = resfctr_mod(F,S,M0);
1227: return [T,S,M0];
1228: }
1229:
1230: def resfctr_mod(F,L,M)
1231: {
1232: for ( T = L, R = []; T != []; T = cdr(T) ) {
1233: U = car(T); MP = U[0]; W = U[1];
1234: for ( A = W, B = F; A != []; A = cdr(cdr(A)) )
1235: B = sremm(subst(B,A[0],A[1]),MP,M);
1236: C = res(var(MP),B,MP) % M;
1237: R = cons(flatten(cdr(modfctr(C,M))),R);
1238: }
1239: return R;
1240: }
1241:
1242: def flatten(L)
1243: {
1244: for ( T = L, R = []; T != []; T = cdr(T) )
1245: R = cons(car(car(T)),R);
1246: return R;
1247: }
1248:
1249: def spm(L,V,M)
1250: {
1251: if ( length(V) == 1 ) {
1252: U = modfctr(car(L),M);
1253: for ( T = cdr(U), R = []; T != []; T = cdr(T) ) {
1254: S = car(T);
1255: R = cons([subst(S[0],var(S[0]),a_),[var(S[0]),a_]],R);
1256: }
1257: return R;
1258: }
1259: L1 = spm(cdr(L),cdr(V),M);
1260: F0 = car(L); V0 = car(V); VR = cdr(V);
1261: for ( T = L1, R = []; T != []; T = cdr(T) ) {
1262: S = car(T);
1263: F1 = subst(F0,S[1]);
1264: U = fctr_mod(F1,V0,S[0],M);
1265: VS = var(S[0]);
1266: for ( W = U; W != []; W = cdr(W) ) {
1267: A = car(car(W));
1268: if ( deg(A,V0) == 1 ) {
1269: A = monic_mod(A,V0,S[0],M);
1270: R = cons([S[0],append([V0,-coef(A,0,V0)],S[1])],R);
1271: } else {
1272: B = pe_mod(A,S[0],M);
1273: MP = B[0]; VMP = var(MP); NV = B[1];
1274: for ( C = S[1], D = []; C != []; C = cdr(cdr(C)) ) {
1275: G = subst(sremm(subst(C[1],VS,NV[1]),MP,M),VMP,VS);
1276: D = append([C[0],G],D);
1277: }
1278: R = cons([subst(MP,VMP,VS),
1279: append([B[2][0],subst(B[2][1],VMP,VS)],D)],R);
1280: }
1281: }
1282: }
1283: return R;
1284: }
1285:
1286: def pe_mod(F,G,M)
1287: {
1288: V = var(G); W = car(setminus(vars(F),[V]));
1289: NG = deg(G,V); NF = deg(F,W); N = NG*NF;
1290: X = prim;
1291: while ( 1 ) {
1292: D = mrandompoly(N,X,M);
1293: if ( irred_check(D,M) )
1294: break;
1295: }
1296: L = fctr_mod(G,V,D,M);
1297: for ( T = L; T != []; T = cdr(T) ) {
1298: U = car(car(T));
1299: if ( deg(U,V) == 1 )
1300: break;
1301: }
1302: U = monic_mod(U,V,D,M); RV = -coef(U,0,V);
1303: L = fctr_mod(sremm(subst(F,V,RV),D,M),W,D,M);
1304: for ( T = L; T != []; T = cdr(T) ) {
1305: U = car(car(T));
1306: if ( deg(U,W) == 1 )
1307: break;
1308: }
1309: U = monic_mod(U,W,D,M); RW = -coef(U,0,W);
1310: return [D,[V,RV],[W,RW]];
1311: }
1312:
1313: def fctr_mod(F,V,D,M)
1314: {
1315: if ( V != x ) {
1316: F = subst(F,V,x); V0 = V; V = x;
1317: } else
1318: V0 = x;
1319: F = monic_mod(F,V,D,M);
1320: L = sqfr_mod(F,V,D,M);
1321: for ( R = [], T = L; T != []; T = cdr(T) ) {
1322: S = car(T); A = S[0]; E = S[1];
1323: B = ddd_mod(A,V,D,M);
1324: R = append(append_mult(B,E),R);
1325: }
1326: if ( V0 != x ) {
1327: for ( R = reverse(R), T = []; R != []; R = cdr(R) )
1328: T = cons([subst(car(R)[0],x,V0),car(R)[1]],T);
1329: R = T;
1330: }
1331: return R;
1332: }
1333:
1334: def append_mult(L,E)
1335: {
1336: for ( T = L, R = []; T != []; T = cdr(T) )
1337: R = cons([car(T),E],R);
1338: return R;
1339: }
1340:
1341: def sqfr_mod(F,V,D,M)
1342: {
1343: setmod(M);
1344: F = sremm(F,D,M);
1345: F1 = sremm(diff(F,V),D,M);
1346: F1 = sremm(F1*inva_mod(D,M,LCOEF(F1)),D,M);
1347: if ( F1 ) {
1348: F2 = ag_mod_single4(F,F1,[D],M);
1349: FLAT = sremm(sdivm(F,F2,M,V),D,M);
1350: I = 0; L = [];
1351: while ( deg(FLAT,V) ) {
1352: while ( 1 ) {
1353: QR = sqrm(F,FLAT,M,V);
1354: if ( !sremm(QR[1],D,M) ) {
1355: F = sremm(QR[0],D,M); I++;
1356: } else
1357: break;
1358: }
1359: if ( !deg(F,V) )
1360: FLAT1 = 1;
1361: else
1362: FLAT1 = ag_mod_single4(F,FLAT,[D],M);
1363: G = sremm(sdivm(FLAT,FLAT1,M,V),D,M);
1364: FLAT = FLAT1;
1365: L = cons([G,I],L);
1366: }
1367: }
1368: if ( deg(F,V) ) {
1369: T = sqfr_mod(pthroot_p_mod(F,V,D,M),V,D,M);
1370: for ( R = []; T != []; T = cdr(T) ) {
1371: H = car(T); R = cons([H[0],M*H[1]],R);
1372: }
1373: } else
1374: R = [];
1375: return append(L,R);
1376: }
1377:
1378: def pthroot_p_mod(F,V,D,M)
1379: {
1380: for ( T = F, R = 0; T; ) {
1381: D1 = deg(T,V); C = coef(T,D1,V); T -= C*V^D1;
1382: R += pthroot_n_mod(C,D,M)*V^idiv(D1,M);
1383: }
1384: return R;
1385: }
1386:
1387: def pthroot_n_mod(C,D,M)
1388: {
1389: pwr_n_mod(C,D,M,deg(D,var(D))-1);
1390: }
1391:
1392: def pwr_n_mod(C,D,M,N)
1393: {
1394: if ( N == 0 )
1395: return 1;
1396: else if ( N == 1 )
1397: return C;
1398: else {
1399: QR = iqr(N,2);
1400: T = pwr_n_mod(C,D,M,QR[0]);
1401: S = sremm(T^2,D,M);
1402: if ( QR[1] )
1403: return sremm(S*C,D,M);
1404: else
1405: return S;
1406: }
1407: }
1408:
1409: def pwr_p_mod(P,A,V,D,M,N)
1410: {
1411: if ( N == 0 )
1412: return 1;
1413: else if ( N == 1 )
1414: return P;
1415: else {
1416: QR = iqr(N,2);
1417: T = pwr_p_mod(P,A,V,D,M,QR[0]);
1418: S = sremm(sremm(sremm(T^2,D,M),A,M,V),D,M);
1419: if ( QR[1] )
1420: return sremm(sremm(sremm(S*P,D,M),A,M,V),D,M);
1421: else
1422: return S;
1423: }
1424: }
1425:
1426: def qmat_mod(F,V,D,M)
1427: {
1428: R = tab_mod(F,V,D,M);
1429: Q = newmat(N,N);
1430: for ( J = 0; J < N; J++ )
1431: for ( I = 0, T = R[J]; I < N; I++ ) {
1432: Q[I][J] = coef(T,I);
1433: }
1434: for ( I = 0; I < N; I++ )
1435: Q[I][I] = (Q[I][I]+(M-1))%M;
1436: return Q;
1437: }
1438:
1439: def tab_mod(F,V,D,M)
1440: {
1441: MD = M^deg(D,var(D));
1442: N = deg(F,V);
1443: F = sremm(F*inva_mod(D,M,coef(F,N,V)),D,M);
1444: R = newvect(N); R[0] = 1;
1445: R[1] = pwr_mod(V,F,V,D,M,MD);
1446: for ( I = 2; I < N; I++ )
1447: R[I] = sremm(sremm(R[1]*R[I-1],F,M),D,M);
1448: return R;
1449: }
1450:
1451: def ddd_mod(F,V,D,M)
1452: {
1453: if ( deg(F,V) == 1 )
1454: return [F];
1455: TAB = tab_mod(F,V,D,M);
1456: for ( I = 1, W = V, L = []; 2*I <= deg(F,V); I++ ) {
1457: for ( T = 0, K = 0; K <= deg(W,V); K++ )
1458: if ( C = coef(W,K,V) )
1459: T = sremm(T+TAB[K]*C,D,M);
1460: W = T;
1461: GCD = ag_mod_single4(F,monic_mod(W-V,V,D,M),[D],M);
1462: if ( deg(GCD,V) ) {
1463: L = append(berlekamp(GCD,V,I,TAB,D,M),L);
1464: F = sremm(sdivm(F,GCD,M,V),D,M);
1465: W = sremm(sremm(W,F,M,V),D,M);
1466: }
1467: }
1468: if ( deg(F,V) )
1469: return cons(F,L);
1470: else
1471: return L;
1472: }
1473:
1474: def monic_mod(F,V,D,M) {
1475: if ( !F || !deg(F,V) )
1476: return F;
1477: return sremm(F*inva_mod(D,M,coef(F,deg(F,V),V)),D,M);
1478: }
1479:
1480: def berlekamp(F,V,E,TAB,D,M)
1481: {
1482: N = deg(F,V);
1483: Q = newmat(N,N);
1484: for ( J = 0; J < N; J++ ) {
1485: T = sremm(sremm(TAB[J],F,M,V),D,M);
1486: for ( I = 0; I < N; I++ ) {
1487: Q[I][J] = coef(T,I);
1488: }
1489: }
1490: for ( I = 0; I < N; I++ )
1491: Q[I][I] = (Q[I][I]+(M-1))%M;
1492: L = nullspace(Q,D,M); MT = L[0]; IND = L[1];
1493: NF0 = N/E;
1494: PS = null_to_poly(MT,IND,V,M);
1495: R = newvect(NF0); R[0] = monic_mod(F,V,D,M);
1496: for ( I = 1, NF = 1; NF < NF0 && I < NF0; I++ ) {
1497: PSI = PS[I];
1498: MP = minipoly_mod(PSI,F,V,D,M);
1499: ROOT = find_root(MP,V,D,M); NR = length(ROOT);
1500: for ( J = 0; J < NF; J++ ) {
1501: if ( deg(R[J],V) == E )
1502: continue;
1503: for ( K = 0; K < NR; K++ ) {
1504: GCD = ag_mod_single4(R[J],PSI-ROOT[K],[D],M);
1505: if ( deg(GCD,V) > 0 && deg(GCD,V) < deg(R[J],V) ) {
1506: Q = sremm(sdivm(R[J],GCD,M,V),D,M);
1507: R[J] = Q; R[NF++] = GCD;
1508: }
1509: }
1510: }
1511: }
1512: return vtol(R);
1513: }
1514:
1515: def null_to_poly(MT,IND,V,M)
1516: {
1517: N = size(MT)[0];
1518: for ( I = 0, J = 0; I < N; I++ )
1519: if ( IND[I] )
1520: J++;
1521: R = newvect(J);
1522: for ( I = 0, L = 0; I < N; I++ ) {
1523: if ( !IND[I] )
1524: continue;
1525: for ( J = K = 0, T = 0; J < N; J++ )
1526: if ( !IND[J] )
1527: T += MT[K++][I]*V^J;
1528: else if ( J == I )
1529: T += (M-1)*V^I;
1530: R[L++] = T;
1531: }
1532: return R;
1533: }
1534:
1535: def minipoly_mod(P,F,V,D,M)
1536: {
1537: L = [[1,1]]; P0 = P1 = 1;
1538: while ( 1 ) {
1539: P0 *= V;
1540: P1 = sremm(sremm(P*P1,F,M,V),D,M);
1541: L1 = lnf_mod(P0,P1,L,V,D,M); NP0 = L1[0]; NP1 = L1[1];
1542: if ( !NP1 )
1543: return NP0;
1544: else
1545: L = lnf_insert([NP0,NP1],L,V);
1546: }
1547: }
1548:
1549: def lnf_mod(P0,P1,L,V,D,M)
1550: {
1551: NP0 = P0; NP1 = P1;
1552: for ( T = L; T != []; T = cdr(T) ) {
1553: Q = car(T);
1554: D1 = deg(NP1,V);
1555: if ( D1 == deg(Q[1],V) ) {
1556: C = coef(Q[1],D1,V);
1557: INV = inva_mod(D,M,M-C); H = sremm(coef(NP1,D1,V)*INV,D,M);
1558: NP0 = sremm(NP0+Q[0]*H,D,M);
1559: NP1 = sremm(NP1+Q[1]*H,D,M);
1560: }
1561: }
1562: return [NP0,NP1];
1563: }
1564:
1565: def lnf_insert(P,L,V)
1566: {
1567: if ( L == [] )
1568: return [P];
1569: else {
1570: P0 = car(L);
1571: if ( deg(P0[1],V) > deg(P[1],V) )
1572: return cons(P0,lnf_insert(P,cdr(L),V));
1573: else
1574: return cons(P,L);
1575: }
1576: }
1577:
1578: def find_root(P,V,D,M)
1579: {
1580: L = c_z(P,V,1,D,M);
1581: for ( T = L, U = []; T != []; T = cdr(T) ) {
1582: S = monic_mod(car(T),V,D,M); U = cons(-coef(S,0,V),U);
1583: }
1584: return U;
1585: }
1586:
1587: def c_z(F,V,E,D,M)
1588: {
1589: N = deg(F,V);
1590: if ( N == E )
1591: return [F];
1592: Q = M^deg(D,var(D));
1593: K = idiv(N,E);
1594: L = [F];
1595: while ( 1 ) {
1596: W = mrandomgfpoly(2*E,V,D,M);
1597: if ( M == 2 ) {
1598: W = monic_mod(tr_mod(W,F,V,D,M,N-1),V,D,M);
1599: } else {
1600: /* W = monic_mod(pwr_p_mod(W,F,V,D,M,idiv(Q^E-1,2))-1,V,D,M); */
1601: /* T = pwr_p_mod(W,F,V,D,M,idiv(Q^E-1,2)); */
1602: T = pwr_mod(W,F,V,D,M,idiv(Q^E-1,2));
1603: W = monic_mod(T-1,V,D,M);
1604: }
1605: if ( !W )
1606: continue;
1607: G = ag_mod_single4(F,W,[D],M);
1608: if ( deg(G,V) && deg(G,V) < N ) {
1609: L1 = c_z(G,V,E,D,M);
1610: L2 = c_z(sremm(sdivm(F,G,M,V),D,M),V,E,D,M);
1611: return append(L1,L2);
1612: }
1613: }
1614: }
1615:
1616: def tr_mod(P,F,V,D,M,N)
1617: {
1618: for ( I = 1, S = P, W = P; I <= N; I++ ) {
1619: W = sremm(sremm(W^2,F,M,V),D,M);
1620: S = sremm(S+W,D,M);
1621: }
1622: return S;
1623: }
1624:
1625: def mrandomgfpoly(N,V,D,M)
1626: {
1627: W = var(D); ND = deg(D,W);
1628: for ( I = N-2, S = V^(N-1); I >= 0; I-- )
1629: S += randompoly(ND,W,M)*V^I;
1630: return S;
1631: }
1632:
1633: def randompoly(N,V,M)
1634: {
1635: for ( I = 0, S = 0; I < N; I++ )
1636: S += (random()%M)*V^I;
1637: return S;
1638: }
1639:
1640: def mrandompoly(N,V,M)
1641: {
1642: for ( I = N-1, S = V^N; I >=0; I-- )
1643: S += (random()%M)*V^I;
1644: return S;
1645: }
1646:
1647: def srem_by_nf(P,B,V,O) {
1648: dp_ord(O); DP = dp_ptod(P,V);
1649: N = length(B); DB = newvect(N);
1650: for ( I = N-1, IL = []; I >= 0; I-- ) {
1651: DB[I] = dp_ptod(B[I],V);
1652: IL = cons(I,IL);
1653: }
1654: L = dp_true_nf(IL,DP,DB,1);
1655: return [dp_dtop(L[0],V),L[1]];
1656: }
1657: end$