![[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.2 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/lib/sp,v 1.1.1.1 1999/12/03 07:39:11 noro Exp $ */
1.1 noro 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;
1.2 ! noro 87: S = reverse(asq(P));
1.1 noro 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];
1.2 ! noro 127: DCSQ = sortfs(asq(N));
1.1 noro 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) )
1.2 ! noro 145: U = cr_gcda(S,G);
1.1 noro 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) )
1.2 ! noro 273: T = union_sort(T,getalgp(C));
1.1 noro 274: return T;
275: }
276: }
277:
1.2 ! noro 278: def getalgtreep(P) {
! 279: if ( type(P) <= 1 )
! 280: return getalgtree(P);
! 281: else {
! 282: for ( V = var(P), I = deg(P,V), T = []; I >= 0; I-- )
! 283: if ( C = coef(P,I) )
! 284: T = union_sort(T,getalgtreep(C));
! 285: return T;
! 286: }
1.1 noro 287: }
288:
1.2 ! noro 289: /* C = union of A and B; A and B is sorted. C should also be sorted. */
! 290:
! 291: def union_sort(A,B)
1.1 noro 292: {
293: if ( A == [] )
1.2 ! noro 294: return B;
! 295: else if ( B == [] )
1.1 noro 296: return A;
1.2 ! noro 297: else {
! 298: A0 = car(A);
! 299: B0 = car(B);
! 300: if ( A0 == B0 )
! 301: return cons(A0,union_sort(cdr(A),cdr(B)));
! 302: else if ( A0 > B0 )
! 303: return cons(A0,union_sort(cdr(A),B));
! 304: else
! 305: return cons(B0,union_sort(A,cdr(B)));
! 306: }
1.1 noro 307: }
308:
309: def invalgp(A)
310: {
311: if ( ntype(A) <= 1 )
312: return 1/A;
313: P0 = defpoly(mainalg(A)); P = algtorat(A);
314: V = var(P0); G1 = P0;
315: G2 = DEG(P)>=DEG(P0)?srem(P,P0):P;
316: for ( H = 1, X = 1, U1 = 0, U2 = 1; deg(G2,V); ) {
317: D = DEG(G1)-DEG(G2); T = LCOEF(G2)^(D+1);
318: L = sqr(G1*T,G2); Q = car(L); R = car(cdr(L));
319: S = U1*T-U2*Q;
320: M = H^D; M1 = M*X;
321: G1 = G2; G2 = sdiv(R,M1);
322: U1 = U2; U2 = sdiv(S,M1);
323: X = LCOEF(G1); H = sdiv(X^D*H,M);
324: }
325: C = invalgp(rattoalg(srem(P*U2,P0)));
326: return C*rattoalg(U2);
327: }
328:
329: def algptorat(P) {
330: if ( type(P) <= 1 )
331: return algtorat(P);
332: else {
333: for ( V = var(P), I = deg(P,V), T = 0; I >= 0; I-- )
334: if ( C = coef(P,I) )
335: T += algptorat(C)*V^I;
336: return T;
337: }
338: }
339:
340: def rattoalgp(P,M) {
341: for ( T = M, S = P; T != []; T = cdr(T) )
342: S = subst(S,algtorat(FIRST(T)),FIRST(T));
343: return S;
344: }
345: def sortfs(L)
346: {
347: #define Factor(a) car(a)
348: #define Mult(a) car(cdr(a))
349: if ( type(TT(L)) == NUM )
350: L = cdr(L);
351: for ( N = 0, T = L; T != []; T = cdr(T), N++ );
352: P = newvect(N); P1 = newvect(N);
353: for ( I = 0, T = L, R = []; T != []; T = cdr(T) )
354: if ( Mult(car(T)) == 1 ) {
355: R = cons(car(T),R); N--;
356: } else {
357: P[I] = car(T); I++;
358: }
359: for ( J = 0, V = var(Factor(P[0])); J < N; J++ ) {
360: for ( K0 = K = J, D = deg(Factor(P[J]),V); K < N; K++ )
361: if ( deg(Factor(P[K]),V) < D ) {
362: K0 = K;
363: D = deg(Factor(P[K]),V);
364: }
365: P1[J] = P[K0];
366: if ( J != K0 )
367: P[K0] = P[J];
368: }
369: for ( I = N - 1; I >= 0; I-- )
370: R = cons(P1[I],R);
371: return R;
372: }
373:
374: def pdiva(P1,P2)
375: {
1.2 ! noro 376: A = union_sort(getalgp(P1),getalgp(P2));
1.1 noro 377: P1 = algptorat(P1); P2 = algptorat(P2);
378: return simpalg(rattoalgp(sdiv(P1*LCOEF(P2)^(DEG(P1)-DEG(P2)+1),P2),A));
379: }
380:
381: def pqra(P1,P2)
382: {
383: if ( type(P2) != POLY )
384: return [P1,0];
385: else if ( (type(P1) != POLY) || (deg(P1,var(P1)) < deg(P2,var(P1))) )
386: return [0,P1];
387: else {
1.2 ! noro 388: A = union_sort(getalgp(P1),getalgp(P2));
1.1 noro 389: P1 = algptorat(P1); P2 = algptorat(P2);
390: L = sqr(P1*LCOEF(P2)^(DEG(P1)-DEG(P2)+1),P2);
391: return [simpalg(rattoalgp(L[0],A)),simpalg(rattoalgp(L[1],A))];
392: }
393: }
394:
395: def pra(P1,P2,AL)
396: {
397: if ( type(P2) != POLY )
398: return 0;
399: else if ( (type(P1) != POLY) || (deg(P1,var(P1)) < deg(P2,var(P1))) )
400: return P1;
401: else {
402: F1 = algptorat(P1); F2 = algptorat(P2); ML = map(defpoly,AL);
403: B = append(reverse(ML),[F2]);
404: V0 = var(P1);
405: V = cons(V0,map(algtorat,AL));
406: G = srem_by_nf(F1,B,V,2);
407: return simpalg(rattoalgp(G[0]/G[1],AL));
408: }
409: }
410:
411: def sort_alg(VL)
412: {
413: N = length(VL); W = newvect(N,VL);
414: for ( I = 0; I < N; I++ ) {
415: for ( M = I, J = I + 1; J < N; J++ )
416: if ( W[J] > W[M] )
417: M = J;
418: if ( I != M ) {
419: T = W[I]; W[I] = W[M]; W[M] = T;
420: }
421: }
422: for ( I = N-1, L = []; I >= 0; I-- )
423: L = cons(W[I],L);
424: return L;
425: }
426:
1.2 ! noro 427: def asq(P)
1.1 noro 428: {
429: P = simpalg(P);
430: if ( type(P) == NUM )
431: return [[1,1]];
432: else if ( getalgp(P) == [] )
433: return sqfr(P);
434: else {
435: V = var(P); N = DEG(P); A = newvect(N+1); B = newvect(N+1);
436: for ( I = 0, F = P; ;I++ ) {
437: if ( type(F) == NUM )
438: break;
439: F1 = diff(F,V);
1.2 ! noro 440: GCD = cr_gcda(F,F1);
1.1 noro 441: FLAT = pdiva(F,GCD);
442: if ( type(GCD) == NUM ) {
443: A[I] = F; B[I] = 1;
444: break;
445: }
446: for ( J = 1, F = GCD; ; J++ ) {
447: L = pqra(F,FLAT); Q = L[0]; R = L[1];
448: if ( R )
449: break;
450: else
451: F = Q;
452: }
453: A[I] = FLAT; B[I] = J;
454: }
455: for ( I = 0, J = 0, L = []; A[I]; I++ ) {
456: J += B[I];
457: if ( A[I+1] )
458: C = pdiva(A[I],A[I+1]);
459: else
460: C = A[I];
461: L = cons([C,J],L);
462: }
463: return L;
464: }
465: }
466:
467: def ufctrhint1(P,HINT)
468: {
469: if ( deg(P,var(P)) == 168 ) {
470: SQ = sqfr(P);
471: if ( length(SQ) == 2 && SQ[1][1] == 1 )
472: return [[1,1],[P,1]];
473: else
474: return ufctrhint(P,HINT);
475: } else
476: return ufctrhint(P,HINT);
477: }
478:
479: def simpcoef(P) {
480: return rattoalgp(ptozp(algptorat(P)),getalgp(P));
481: }
482:
483: def ufctrhint_heuristic(P,HINT,PP,SHIFT) {
484: V = var(P); D = deg(P,V);
485: if ( D == HINT )
486: return [[P,1]];
487: for ( S = 0, L = SHIFT, AL = [], K = 1; L != []; L = cdr(L) ) {
488: A = car(L)[0]; S += A*car(L)[1]; AL = cons(A,AL);
489: K *= deg(defpoly(A),algtorat(A));
490: }
491: PPP = simpcoef(simpalg(subst(PP,V,V-S)));
492: for ( T = P-coef(P,D)*V^D, G = D; T; T -= coef(T,DT)*V^DT )
493: G = igcd(G,DT=deg(T,V));
494: if ( G == 1 ) {
495: if ( K*deg(PPP,V) != deg(P,V) )
1.2 ! noro 496: PPP = cr_gcda(PPP,P);
1.1 noro 497: return ufctrhint2(P,HINT,PPP,AL);
498: } else {
499: for ( S = 0, T = P; T; T -= coef(T,DT)*V^DT ) {
500: DT = deg(T,V);
501: S += coef(T,DT)*V^(DT/G);
502: }
503: L = fctr(S);
504: for ( DC = [car(L)], L = cdr(L); L != []; L = cdr(L) ) {
505: H = subst(car(car(L)),V,V^G);
1.2 ! noro 506: HH = cr_gcda(PPP,H);
1.1 noro 507: T = ufctrhint2(H,HINT,HH,AL);
508: DC = append(DC,T);
509: }
510: return DC;
511: }
512: }
513:
514: def ufctrhint2(P,HINT,PP,AL)
515: {
516: if ( deg(P,var(P)) == HINT )
517: return [[P,1]];
518: if ( AL == [] )
519: return ufctrhint(P,HINT);
520: L = resfctr(algptorat(PP),map(defpoly,AL),map(algtorat,AL),P);
521: for ( T = reverse(L[1]), DL = []; T != []; T = cdr(T) )
522: DL = cons(deg(car(car(T)),a_),DL);
523: return resfmain(P,L[2],L[0],DL);
524: }
525:
526: def res_det(V,P1,P2)
527: {
528: D1 = deg(P1,V); D2 = deg(P2,V);
529: M = newmat(D1+D2,D1+D2);
530: for ( J = 0; J <= D2; J++ )
531: M[0][J] = coef(P2,D2-J,V);
532: for ( I = 1; I < D1; I++ )
533: for ( J = 0; J <= D2; J++ )
534: M[I][I+J] = M[0][J];
535: for ( J = 0; J <= D1; J++ )
536: M[D1][J] = coef(P1,D1-J,V);
537: for ( I = 1; I < D2; I++ )
538: for ( J = 0; J <= D1; J++ )
539: M[D1+I][I+J] = M[D1][J];
540: return det(M);
541: }
542:
543: def norm_ch1(V0,VM,P,P0,PR) {
544: D = deg(P,V0); D0 = deg(P0,V0); DM = deg(P,VM); N = DM*D0;
545: X = newvect(N+1); V = newvect(N+1); U = newvect(N+1);
546: Min = -idiv(N,2);
547: C = coef(P,D,V0);
548: for ( I = J = 0; I <= N; J++ ) {
549: if ( PRINT )
550: print([J,N]);
551: T=J+Min;
552: if ( subst(C,VM,T) ) {
553: U[I] = srem(res(V0,subst(P,VM,T),P0),PR);
554: X[I++] = T;
555: }
556: }
557: for ( I = 1, M = 1, S = V[0] = U[0]; I <= N; I++ ) {
558: for ( J = 0, T = U[I]; J < I; J++ )
559: T = sdiv(T-V[J],X[I]-X[J]);
560: V[I] = T;
561: M *= (VM-X[I-1]);
562: S += T*M;
563: }
564: return S;
565: }
566:
567: def norm_ch2(V0,VM,P,P0,PR) {
568: D = deg(P,V0); D0 = deg(P0,V0); DM = deg(P,VM); N = DM*D0;
569: X = newvect(N+1); V = newvect(N+1); U = newvect(N+1);
570: Min = -idiv(N,2);
571: C = coef(P,D,V0);
572: for ( I = J = 0; I <= N; J++ ) {
573: T=J+Min;
574: if ( subst(C,VM,T) ) {
575: U[I] = srem(res_det(V0,subst(P,VM,T),P0),PR);
576: X[I++] = T;
577: }
578: }
579: for ( I = 1, M = 1, S = V[0] = U[0]; I <= N; I++ ) {
580: for ( J = 0, T = U[I]; J < I; J++ )
581: T = sdiv(T-V[J],X[I]-X[J]);
582: V[I] = T;
583: M *= (VM-X[I-1]);
584: S += T*M;
585: }
586: return S;
587: }
588:
589: def res_ch1(V0,VM,P,P0) {
590: D = deg(P,V0); D0 = deg(P0,V0); N = deg(P,VM)*D0+deg(P0,VM)*D;
591: X = newvect(N+1); V = newvect(N+1); U = newvect(N+1);
592: Min = -idiv(N,2);
593: C = coef(P,D,V0); C0 = coef(P0,D0,V0);
594: for ( I = J = 0; I <= N; J++ ) {
595: if ( PRINT )
596: print([J,N]);
597: T=J+Min;
598: if ( subst(C,VM,T) && subst(C0,VM,T) ) {
599: U[I] = res(V0,subst(P,VM,T),subst(P0,VM,T));
600: X[I++] = T;
601: }
602: }
603: for ( I = 1, M = 1, S = V[0] = U[0]; I <= N; I++ ) {
604: for ( J = 0, T = U[I]; J < I; J++ )
605: T = sdiv(T-V[J],X[I]-X[J]);
606: V[I] = T;
607: M *= (VM-X[I-1]);
608: S += T*M;
609: }
610: return S;
611: }
612:
613: def res_ch(V0,VM,P,P0) {
614: D = deg(P,V0); D0 = deg(P0,V0); N = deg(P,VM)*D0+deg(P0,VM)*D;
615: X = newvect(N+1); V = newvect(N+1); U = newvect(N+1);
616: Min = -idiv(N,2);
617: C = coef(P,D,V0); C0 = coef(P0,D0,V0);
618: for ( I = J = 0; I <= N; J++ ) {
619: T=J+Min;
620: if ( subst(C,VM,T) && subst(C0,VM,T) ) {
621: U[I] = res_det(V0,subst(P,VM,T),subst(P0,VM,T));
622: X[I++] = T;
623: }
624: }
625: for ( I = 1, M = 1, S = V[0] = U[0]; I <= N; I++ ) {
626: for ( J = 0, T = U[I]; J < I; J++ )
627: T = sdiv(T-V[J],X[I]-X[J]);
628: V[I] = T;
629: M *= (VM-X[I-1]);
630: S += T*M;
631: }
632: return S;
633: }
634:
635: def norm_ch2_lag(V,VM,P,P0,PR) {
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: R = srem(R,PR);
643: T = sdiv(A,VM-I-Min);
644: S += R*T/subst(T,VM,I+Min);
645: }
646: return S;
647: }
648:
649: def norm_ch_lag(V,VM,P,P0) {
650: D0 = deg(P0,V); DM = deg(P,VM); N = DM*D0;
651: Min = -idiv(N,2);
652: for ( A = 1, I = 0; I <= N; I++ )
653: A *= (VM-I-Min);
654: for ( I = 0, S = 0; I <= N; I++ ) {
655: R = res_det(V,subst(P,VM,I+Min),P0);
656: T = sdiv(A,VM-I-Min);
657: S += R*T/subst(T,VM,I+Min);
658: }
659: return S;
660: }
661:
1.2 ! noro 662: def cr_gcda(P1,P2)
1.1 noro 663: {
664: if ( !(V = var(P1)) || !var(P2) )
665: return 1;
1.2 ! noro 666: EXT = union_sort(getalgtreep(P1),getalgtreep(P2));
! 667: if ( EXT == [] )
1.1 noro 668: return gcd(P1,P2);
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 discr(P) {
751: V = var(P);
752: return res(V,P,diff(P,V));
753: }
754:
755: def sp_minipoly(A,EXT)
756: {
757: while ( car(EXT) != A )
758: EXT = cdr(EXT);
759: for ( M = x-A; EXT != []; EXT = cdr(EXT) )
760: M = sp_norm(car(EXT),x,M,EXT);
761: F = sqfr(M);
762: return F[1][0];
763: }
764:
765: def cr(F1,M1,F2,M2)
766: {
767: K = inv(M1 % M2,M2);
768: M3 = M1*M2;
769: F3 = (F1 + (F2-(F1%M2))*K*M1) % M3;
770: return [F3,M3];
771: }
772:
773: #define ABS(a) ((a)>=0?(a):(-a))
774:
775: #if 0
776: def calcb(M) {
777: setprec(800);
778: return pari(floor,eval((M/2)^(1/2)));
779: }
780: #endif
781:
782: def calcb(M) {
783: N = 2*M;
784: T = sp_sqrt(N);
785: if ( T^2 <= N && N < (T+1)^2 )
786: return idiv(T,2);
787: else
788: error("afo");
789: }
790:
791: def sp_sqrt(A) {
792: for ( J = 0, T = A; T >= 2^27; J++ ) {
793: T = idiv(T,2^27)+1;
794: }
795: for ( I = 0; T >= 2; I++ ) {
796: S = idiv(T,2);
797: if ( T = S+S )
798: T = S;
799: else
800: T = S+1;
801: }
802: X = (2^27)^idiv(J,2)*2^idiv(I,2);
803: while ( 1 ) {
804: if ( (Y=X^2) < A )
805: X += X;
806: else if ( Y == A )
807: return X;
808: else
809: break;
810: }
811: while ( 1 )
812: if ( (Y = X^2) <= A )
813: return X;
814: else
815: X = idiv(A + Y,2*X);
816: }
817:
818: def intptoratp(P,M,B) {
819: if ( type(P) == 1 ) {
820: L = inttorat(P,M,B);
821: if ( L == 0 )
822: return [];
823: else
824: return L[0]/L[1];
825: } else {
826: V = var(P);
827: S = 0;
828: while ( P ) {
829: D = deg(P,V);
830: C = coef(P,D,V);
831: T = intptoratp(C,M,B);
832: if ( T == [] )
833: return [];
834: S += T*V^D;
835: P -= C*V^D;
836: }
837: return S;
838: }
839: }
840:
841: def ltoalg(L) {
842: F = L[0]; V = reverse(L[1]);
843: N = length(V)-1;
844: for ( I = 0, G = F; I < N; I++ ) {
845: D = car(G);
846: A = newalg(D); V = var(D);
847: for ( G = reverse(cdr(G)), T = []; G != []; G = cdr(G) )
848: T = cons(subst(car(G),V,A),T);
849: G = T;
850: }
851: return G;
852: }
853:
854: /*
855: def ag_mod(F1,F2,D,MOD)
856: {
857: if ( length(D) == 1 )
858: return ag_mod_single(F1,F2,D,MOD);
859: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
860: if ( D1 < D2 ) {
861: T = F1; F1 = F2; F2 = T;
862: T = D1; D1 = D2; D2 = T;
863: }
864: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
865: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
866: for ( T = reverse(D), S = []; T != []; T = cdr(T) ) {
867: U = car(T);
868: S = cons((inv(LCOEF(U),MOD)*U) % MOD,S);
869: }
870: D = S;
871: while ( 1 ) {
872: F = srem_simp_mod(F1,F2,V,D,MOD);
873: if ( !F )
874: return F2;
875: if ( !deg(F,V) )
876: return 1;
877: C = LCOEF(F);
878: INV = inverse_by_gr_mod(C,D,MOD);
879: if ( !INV )
880: return 0;
881: F = simp_by_dp_mod(F*INV,D,MOD);
882: F = (inv(LCOEF(F),MOD)*F) % MOD;
883: F1 = F2; F2 = F;
884: }
885: }
886: */
887:
888: def ag_mod(F1,F2,D,VL,MOD)
889: {
890: if ( length(D) == 1 )
891: return ag_mod_single6(F1,F2,D,MOD);
892: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
893: if ( D1 < D2 ) {
894: T = F1; F1 = F2; F2 = T;
895: T = D1; D1 = D2; D2 = T;
896: }
897: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
898: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
899: for ( T = reverse(D), S = []; T != []; T = cdr(T) ) {
900: U = car(T);
901: S = cons((inv(LCOEF(U),MOD)*U) % MOD,S);
902: }
903: D = S;
904: VL = cons(V,VL); B = append([F1,F2],D); N = length(VL);
905: while ( 1 ) {
906: FLAGS = dp_gr_flags(); dp_gr_flags(["Reverse",1,"NoSugar",1]);
907: G = dp_gr_mod_main(B,VL,0,MOD,Ord);
908: dp_gr_flags(FLAGS);
909: if ( length(G) == 1 )
910: return 1;
911: if ( length(G) == N ) {
912: for ( T = G; T != []; T = cdr(T) )
913: if ( member(V,vars(car(T))) )
914: return car(T);
915: }
916: }
917: }
918:
919: def srem_simp_mod(F1,F2,V,D,MOD)
920: {
921: D2 = deg(F2,V); C = coef(F2,D2);
922: while ( (D1 = deg(F1,V)) >= D2 ) {
923: F1 -= coef(F1,D1)/C*V^(D1-D2)*F2;
924: F1 = simp_by_dp_mod(F1,D,MOD);
925: }
926: return F1;
927: }
928:
929: def ag_mod_single(F1,F2,D,MOD)
930: {
931: TD = TI = TM = 0;
932: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
933: if ( D1 < D2 ) {
934: T = F1; F1 = F2; F2 = T;
935: T = D1; D1 = D2; D2 = T;
936: }
937: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
938: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
939: D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
940: FLAGS = dp_gr_flags(); dp_gr_flags(["Reverse",1,"NoSugar",1]);
941: G = dp_gr_mod_main([F1,F2,D],[V,var(D)],0,MOD,2);
942: dp_gr_flags(FLAGS);
943: if ( length(G) == 1 )
944: return 1;
945: if ( length(G) != 2 )
946: return 0;
947: if ( vars(G[0]) == [var(D)] )
948: return G[1];
949: else
950: return G[0];
951: }
952:
953: def ag_mod_single2(F1,F2,D,MOD)
954: {
955: TD = TI = TM = 0;
956: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
957: if ( D1 < D2 ) {
958: T = F1; F1 = F2; F2 = T;
959: T = D1; D1 = D2; D2 = T;
960: }
961: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
962: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
963: D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
964: while ( 1 ) {
965: T0 = time()[0];
966: F = srem((srem(F1,F2) % MOD),D) % MOD;
967: TD += time()[0] - T0;
968: if ( !F ) {
969: if ( dp_gr_print() )
970: print(["TD",TD,"TM",TM,"TI",TI]);
971: return F2;
972: }
973: if ( !deg(F,V) ) {
974: if ( dp_gr_print() )
975: print(["TD",TD,"TM",TM,"TI",TI]);
976: return 1;
977: }
978: C = LCOEF(F);
979: T0 = time()[0];
980: INV = inva_mod(D,MOD,C);
981: TI += time()[0] - T0;
982: if ( !INV )
983: return 0;
984: T0 = time()[0];
985: F = remc_mod((INV*F) % MOD,D,MOD);
986: TM += time()[0] - T0;
987: F1 = F2; F2 = F;
988: }
989: }
990:
991: def ag_mod_single3(F1,F2,D,MOD)
992: {
993: TD = TI = TM = 0;
994: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
995: if ( D1 < D2 ) {
996: T = F1; F1 = F2; F2 = T;
997: T = D1; D1 = D2; D2 = T;
998: }
999: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
1000: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
1001: D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
1002: while ( 1 ) {
1003: if ( !D2 )
1004: return 1;
1005: while ( D1 >= D2 ) {
1006: F = srem((coef(F2,D2,V)*F1-coef(F1,D1,V)*F2*V^(D1-D2))%MOD,D)%MOD;
1007: F1 = F; D1 = deg(F1,V);
1008: }
1009: if ( !F1 ) {
1010: INV = inva_mod(D,MOD,coef(F2,D2,V));
1011: if ( dp_gr_print() )
1012: print(".");
1013: return srem((INV*F2) % MOD,D)%MOD;
1014: } else {
1015: T = F1; F1 = F2; F2 = T;
1016: T = D1; D1 = D2; D2 = T;
1017: }
1018: }
1019: }
1020:
1021: def ag_mod_single4(F1,F2,D,MOD)
1022: {
1023: if ( !F1 )
1024: return F2;
1025: if ( !F2 )
1026: return F1;
1027: TD = TI = TM = TR = 0;
1028: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
1029: if ( D1 < D2 ) {
1030: T = F1; F1 = F2; F2 = T;
1031: T = D1; D1 = D2; D2 = T;
1032: }
1033: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
1034: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
1035: D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
1036: while ( 1 ) {
1037: T0 = time()[0]; R = srem(F1,F2); TR += time()[0] - T0;
1038: T0 = time()[0]; F = srem(R % MOD,D) % MOD; TD += time()[0] - T0;
1039: if ( !F ) {
1040: if ( dp_gr_print() )
1041: print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
1042: return F2;
1043: }
1044: if ( !deg(F,V) ) {
1045: if ( dp_gr_print() )
1046: print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
1047: return 1;
1048: }
1049: C = LCOEF(F);
1050: T0 = time()[0]; INV = inva_mod(D,MOD,C); TI += time()[0] - T0;
1051: if ( !INV )
1052: return 0;
1053: T0 = time()[0]; F = srem((INV*F) % MOD,D) % MOD; TM += time()[0] - T0;
1054: F1 = F2; F2 = F;
1055: }
1056: }
1057:
1058: def ag_mod_single5(F1,F2,D,MOD)
1059: {
1060: TD = TI = TM = TR = 0;
1061: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
1062: if ( D1 < D2 ) {
1063: T = F1; F1 = F2; F2 = T;
1064: T = D1; D1 = D2; D2 = T;
1065: }
1066: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
1067: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
1068: D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
1069: while ( 1 ) {
1070: T0 = time()[0];
1071: D1 = deg(F1,V); D2 = deg(F2,V); F = F1;
1072: while ( D1 >= D2 ) {
1073: R = (F-coef(F,D1,V)*F2*V^(D1-D2))%MOD;
1074: D1 = deg(R,V); HC = coef(R,D1,V);
1075: F = (R - HC*V^D1) + (srem(HC,D)%MOD)*V^D1;
1076: }
1077: TR += time()[0] - T0;
1078: T0 = time()[0]; F = srem(R % MOD,D) % MOD; TD += time()[0] - T0;
1079: if ( !F ) {
1080: if ( dp_gr_print() )
1081: print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
1082: return F2;
1083: }
1084: if ( !deg(F,V) ) {
1085: if ( dp_gr_print() )
1086: print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
1087: return 1;
1088: }
1089: C = LCOEF(F);
1090: T0 = time()[0]; INV = inva_mod(D,MOD,C); TI += time()[0] - T0;
1091: if ( !INV )
1092: return 0;
1093: T0 = time()[0]; F = srem((INV*F) % MOD,D) % MOD; TM += time()[0] - T0;
1094: F1 = F2; F2 = F;
1095: }
1096: }
1097:
1098: def ag_mod_single6(F1,F2,D,MOD)
1099: {
1100: TD = TI = TM = TR = 0;
1101: V = var(F1); D1 = deg(F1,V); D2 = deg(F2,V);
1102: if ( D1 < D2 ) {
1103: T = F1; F1 = F2; F2 = T;
1104: T = D1; D1 = D2; D2 = T;
1105: }
1106: F1 = (inv(LCOEF(F1),MOD)*F1) % MOD;
1107: F2 = (inv(LCOEF(F2),MOD)*F2) % MOD;
1108: D = (inv(LCOEF(car(D)),MOD)*car(D)) % MOD;
1109: while ( 1 ) {
1110: T0 = time()[0];
1111: D1 = deg(F1,V); D2 = deg(F2,V); F = F1;
1112: while ( D1 >= D2 ) {
1113: R = (F-coef(F,D1,V)*F2*V^(D1-D2))%MOD;
1114: D1 = deg(R,V); HC = coef(R,D1,V);
1115: /* F = (R - HC*V^D1) + (srem_mod(HC,D,MOD))*V^D1; */
1116: F = remc_mod(R,D,MOD);
1117: }
1118: TR += time()[0] - T0;
1119: T0 = time()[0]; F = remc_mod(R%MOD,D,MOD); TD += time()[0] - T0;
1120: if ( !F ) {
1121: if ( dp_gr_print() )
1122: print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
1123: return F2;
1124: }
1125: if ( !deg(F,V) ) {
1126: if ( dp_gr_print() )
1127: print(["TD",TD,"TM",TM,"TI",TI,"TR",TR]);
1128: return 1;
1129: }
1130: C = LCOEF(F);
1131: T0 = time()[0]; INV = inva_mod(D,MOD,C); TI += time()[0] - T0;
1132: if ( !INV )
1133: return 0;
1134: T0 = time()[0]; F = remc_mod((INV*F)%MOD,D,MOD); TM += time()[0] - T0;
1135: F1 = F2; F2 = F;
1136: }
1137: }
1138:
1139: def inverse_by_gr_mod(C,D,MOD)
1140: {
1141: Ord = 2;
1142: dp_gr_flags(["NoSugar",1]);
1143: G = dp_gr_mod_main(cons(x*C-1,D),cons(x,vars(D)),0,MOD,Ord);
1144: dp_gr_flags(["NoSugar",0]);
1145: if ( length(G) == 1 )
1146: return 1;
1147: else if ( length(G) == length(D)+1 ) {
1148: for ( T = G; T != []; T = cdr(T) )
1149: if ( member(x,vars(car(G))) )
1150: break;
1151: T = car(G);
1152: if ( type(coef(T,1,x)) != NUM )
1153: return 0;
1154: else
1155: return coef(T,0,x);
1156: } else
1157: return 0;
1158: }
1159:
1160: def simp_by_dp(F,D)
1161: {
1162: for ( T = D; T != []; T = cdr(T) )
1163: F = srem(F,car(T));
1164: return F;
1165: }
1166:
1167: def simp_by_dp_mod(F,D,MOD)
1168: {
1169: F %= MOD;
1170: for ( T = D; T != []; T = cdr(T) )
1171: F = srem(F,car(T)) % MOD;
1172: return F;
1173: }
1174:
1175: def remc_mod(P,D,M)
1176: {
1177: V = var(P);
1178: if ( !V || V == var(D) )
1179: return srem_mod(P,D,M);
1180: for ( I = deg(P,V), S = 0; I >= 0; I-- )
1181: if ( C = coef(P,I,V) )
1182: S += srem_mod(C,D,M)*V^I;
1183: return S;
1184: }
1185:
1186: def rem_mod(C,D,M)
1187: {
1188: V = var(D);
1189: D2 = deg(D,V);
1190: while ( (D1 = deg(C,V)) >= D2 ) {
1191: C -= (D*V^(D1-D2)*coef(C,D1,V))%M;
1192: C %= M;
1193: }
1194: return C;
1195: }
1196:
1197: def resfctr(F,L,V,N)
1198: {
1199: N = ptozp(N);
1200: V0 = var(N);
1201: DN = diff(N,V0);
1202: for ( I = 0, J = 2, Len = deg(N,V0)+1; I < 5; J++ ) {
1203: M = prime(J);
1204: G = gcd(N,DN,M);
1205: if ( !deg(G,V0) ) {
1206: I++;
1207: T = nfctr_mod(N,M);
1208: if ( T < Len ) {
1209: Len = T; M0 = M;
1210: }
1211: }
1212: }
1213: S = spm(L,V,M0);
1214: T = resfctr_mod(F,S,M0);
1215: return [T,S,M0];
1216: }
1217:
1218: def resfctr_mod(F,L,M)
1219: {
1220: for ( T = L, R = []; T != []; T = cdr(T) ) {
1221: U = car(T); MP = U[0]; W = U[1];
1222: for ( A = W, B = F; A != []; A = cdr(cdr(A)) )
1223: B = sremm(subst(B,A[0],A[1]),MP,M);
1224: C = res(var(MP),B,MP) % M;
1225: R = cons(flatten(cdr(modfctr(C,M))),R);
1226: }
1227: return R;
1228: }
1229:
1230: def flatten(L)
1231: {
1232: for ( T = L, R = []; T != []; T = cdr(T) )
1233: R = cons(car(car(T)),R);
1234: return R;
1235: }
1236:
1237: def spm(L,V,M)
1238: {
1239: if ( length(V) == 1 ) {
1240: U = modfctr(car(L),M);
1241: for ( T = cdr(U), R = []; T != []; T = cdr(T) ) {
1242: S = car(T);
1243: R = cons([subst(S[0],var(S[0]),a_),[var(S[0]),a_]],R);
1244: }
1245: return R;
1246: }
1247: L1 = spm(cdr(L),cdr(V),M);
1248: F0 = car(L); V0 = car(V); VR = cdr(V);
1249: for ( T = L1, R = []; T != []; T = cdr(T) ) {
1250: S = car(T);
1251: F1 = subst(F0,S[1]);
1252: U = fctr_mod(F1,V0,S[0],M);
1253: VS = var(S[0]);
1254: for ( W = U; W != []; W = cdr(W) ) {
1255: A = car(car(W));
1256: if ( deg(A,V0) == 1 ) {
1257: A = monic_mod(A,V0,S[0],M);
1258: R = cons([S[0],append([V0,-coef(A,0,V0)],S[1])],R);
1259: } else {
1260: B = pe_mod(A,S[0],M);
1261: MP = B[0]; VMP = var(MP); NV = B[1];
1262: for ( C = S[1], D = []; C != []; C = cdr(cdr(C)) ) {
1263: G = subst(sremm(subst(C[1],VS,NV[1]),MP,M),VMP,VS);
1264: D = append([C[0],G],D);
1265: }
1266: R = cons([subst(MP,VMP,VS),
1267: append([B[2][0],subst(B[2][1],VMP,VS)],D)],R);
1268: }
1269: }
1270: }
1271: return R;
1272: }
1273:
1274: def pe_mod(F,G,M)
1275: {
1276: V = var(G); W = car(setminus(vars(F),[V]));
1277: NG = deg(G,V); NF = deg(F,W); N = NG*NF;
1278: X = prim;
1279: while ( 1 ) {
1280: D = mrandompoly(N,X,M);
1281: if ( irred_check(D,M) )
1282: break;
1283: }
1284: L = fctr_mod(G,V,D,M);
1285: for ( T = L; T != []; T = cdr(T) ) {
1286: U = car(car(T));
1287: if ( deg(U,V) == 1 )
1288: break;
1289: }
1290: U = monic_mod(U,V,D,M); RV = -coef(U,0,V);
1291: L = fctr_mod(sremm(subst(F,V,RV),D,M),W,D,M);
1292: for ( T = L; T != []; T = cdr(T) ) {
1293: U = car(car(T));
1294: if ( deg(U,W) == 1 )
1295: break;
1296: }
1297: U = monic_mod(U,W,D,M); RW = -coef(U,0,W);
1298: return [D,[V,RV],[W,RW]];
1299: }
1300:
1301: def fctr_mod(F,V,D,M)
1302: {
1303: if ( V != x ) {
1304: F = subst(F,V,x); V0 = V; V = x;
1305: } else
1306: V0 = x;
1307: F = monic_mod(F,V,D,M);
1308: L = sqfr_mod(F,V,D,M);
1309: for ( R = [], T = L; T != []; T = cdr(T) ) {
1310: S = car(T); A = S[0]; E = S[1];
1311: B = ddd_mod(A,V,D,M);
1312: R = append(append_mult(B,E),R);
1313: }
1314: if ( V0 != x ) {
1315: for ( R = reverse(R), T = []; R != []; R = cdr(R) )
1316: T = cons([subst(car(R)[0],x,V0),car(R)[1]],T);
1317: R = T;
1318: }
1319: return R;
1320: }
1321:
1322: def append_mult(L,E)
1323: {
1324: for ( T = L, R = []; T != []; T = cdr(T) )
1325: R = cons([car(T),E],R);
1326: return R;
1327: }
1328:
1329: def sqfr_mod(F,V,D,M)
1330: {
1331: setmod(M);
1332: F = sremm(F,D,M);
1333: F1 = sremm(diff(F,V),D,M);
1334: F1 = sremm(F1*inva_mod(D,M,LCOEF(F1)),D,M);
1335: if ( F1 ) {
1336: F2 = ag_mod_single4(F,F1,[D],M);
1337: FLAT = sremm(sdivm(F,F2,M,V),D,M);
1338: I = 0; L = [];
1339: while ( deg(FLAT,V) ) {
1340: while ( 1 ) {
1341: QR = sqrm(F,FLAT,M,V);
1342: if ( !sremm(QR[1],D,M) ) {
1343: F = sremm(QR[0],D,M); I++;
1344: } else
1345: break;
1346: }
1347: if ( !deg(F,V) )
1348: FLAT1 = 1;
1349: else
1350: FLAT1 = ag_mod_single4(F,FLAT,[D],M);
1351: G = sremm(sdivm(FLAT,FLAT1,M,V),D,M);
1352: FLAT = FLAT1;
1353: L = cons([G,I],L);
1354: }
1355: }
1356: if ( deg(F,V) ) {
1357: T = sqfr_mod(pthroot_p_mod(F,V,D,M),V,D,M);
1358: for ( R = []; T != []; T = cdr(T) ) {
1359: H = car(T); R = cons([H[0],M*H[1]],R);
1360: }
1361: } else
1362: R = [];
1363: return append(L,R);
1364: }
1365:
1366: def pthroot_p_mod(F,V,D,M)
1367: {
1368: for ( T = F, R = 0; T; ) {
1369: D1 = deg(T,V); C = coef(T,D1,V); T -= C*V^D1;
1370: R += pthroot_n_mod(C,D,M)*V^idiv(D1,M);
1371: }
1372: return R;
1373: }
1374:
1375: def pthroot_n_mod(C,D,M)
1376: {
1377: pwr_n_mod(C,D,M,deg(D,var(D))-1);
1378: }
1379:
1380: def pwr_n_mod(C,D,M,N)
1381: {
1382: if ( N == 0 )
1383: return 1;
1384: else if ( N == 1 )
1385: return C;
1386: else {
1387: QR = iqr(N,2);
1388: T = pwr_n_mod(C,D,M,QR[0]);
1389: S = sremm(T^2,D,M);
1390: if ( QR[1] )
1391: return sremm(S*C,D,M);
1392: else
1393: return S;
1394: }
1395: }
1396:
1397: def pwr_p_mod(P,A,V,D,M,N)
1398: {
1399: if ( N == 0 )
1400: return 1;
1401: else if ( N == 1 )
1402: return P;
1403: else {
1404: QR = iqr(N,2);
1405: T = pwr_p_mod(P,A,V,D,M,QR[0]);
1406: S = sremm(sremm(sremm(T^2,D,M),A,M,V),D,M);
1407: if ( QR[1] )
1408: return sremm(sremm(sremm(S*P,D,M),A,M,V),D,M);
1409: else
1410: return S;
1411: }
1412: }
1413:
1414: def qmat_mod(F,V,D,M)
1415: {
1416: R = tab_mod(F,V,D,M);
1417: Q = newmat(N,N);
1418: for ( J = 0; J < N; J++ )
1419: for ( I = 0, T = R[J]; I < N; I++ ) {
1420: Q[I][J] = coef(T,I);
1421: }
1422: for ( I = 0; I < N; I++ )
1423: Q[I][I] = (Q[I][I]+(M-1))%M;
1424: return Q;
1425: }
1426:
1427: def tab_mod(F,V,D,M)
1428: {
1429: MD = M^deg(D,var(D));
1430: N = deg(F,V);
1431: F = sremm(F*inva_mod(D,M,coef(F,N,V)),D,M);
1432: R = newvect(N); R[0] = 1;
1433: R[1] = pwr_mod(V,F,V,D,M,MD);
1434: for ( I = 2; I < N; I++ )
1435: R[I] = sremm(sremm(R[1]*R[I-1],F,M),D,M);
1436: return R;
1437: }
1438:
1439: def ddd_mod(F,V,D,M)
1440: {
1441: if ( deg(F,V) == 1 )
1442: return [F];
1443: TAB = tab_mod(F,V,D,M);
1444: for ( I = 1, W = V, L = []; 2*I <= deg(F,V); I++ ) {
1445: for ( T = 0, K = 0; K <= deg(W,V); K++ )
1446: if ( C = coef(W,K,V) )
1447: T = sremm(T+TAB[K]*C,D,M);
1448: W = T;
1449: GCD = ag_mod_single4(F,monic_mod(W-V,V,D,M),[D],M);
1450: if ( deg(GCD,V) ) {
1451: L = append(berlekamp(GCD,V,I,TAB,D,M),L);
1452: F = sremm(sdivm(F,GCD,M,V),D,M);
1453: W = sremm(sremm(W,F,M,V),D,M);
1454: }
1455: }
1456: if ( deg(F,V) )
1457: return cons(F,L);
1458: else
1459: return L;
1460: }
1461:
1462: def monic_mod(F,V,D,M) {
1463: if ( !F || !deg(F,V) )
1464: return F;
1465: return sremm(F*inva_mod(D,M,coef(F,deg(F,V),V)),D,M);
1466: }
1467:
1468: def berlekamp(F,V,E,TAB,D,M)
1469: {
1470: N = deg(F,V);
1471: Q = newmat(N,N);
1472: for ( J = 0; J < N; J++ ) {
1473: T = sremm(sremm(TAB[J],F,M,V),D,M);
1474: for ( I = 0; I < N; I++ ) {
1475: Q[I][J] = coef(T,I);
1476: }
1477: }
1478: for ( I = 0; I < N; I++ )
1479: Q[I][I] = (Q[I][I]+(M-1))%M;
1480: L = nullspace(Q,D,M); MT = L[0]; IND = L[1];
1481: NF0 = N/E;
1482: PS = null_to_poly(MT,IND,V,M);
1483: R = newvect(NF0); R[0] = monic_mod(F,V,D,M);
1484: for ( I = 1, NF = 1; NF < NF0 && I < NF0; I++ ) {
1485: PSI = PS[I];
1486: MP = minipoly_mod(PSI,F,V,D,M);
1487: ROOT = find_root(MP,V,D,M); NR = length(ROOT);
1488: for ( J = 0; J < NF; J++ ) {
1489: if ( deg(R[J],V) == E )
1490: continue;
1491: for ( K = 0; K < NR; K++ ) {
1492: GCD = ag_mod_single4(R[J],PSI-ROOT[K],[D],M);
1493: if ( deg(GCD,V) > 0 && deg(GCD,V) < deg(R[J],V) ) {
1494: Q = sremm(sdivm(R[J],GCD,M,V),D,M);
1495: R[J] = Q; R[NF++] = GCD;
1496: }
1497: }
1498: }
1499: }
1500: return vtol(R);
1501: }
1502:
1503: def null_to_poly(MT,IND,V,M)
1504: {
1505: N = size(MT)[0];
1506: for ( I = 0, J = 0; I < N; I++ )
1507: if ( IND[I] )
1508: J++;
1509: R = newvect(J);
1510: for ( I = 0, L = 0; I < N; I++ ) {
1511: if ( !IND[I] )
1512: continue;
1513: for ( J = K = 0, T = 0; J < N; J++ )
1514: if ( !IND[J] )
1515: T += MT[K++][I]*V^J;
1516: else if ( J == I )
1517: T += (M-1)*V^I;
1518: R[L++] = T;
1519: }
1520: return R;
1521: }
1522:
1523: def minipoly_mod(P,F,V,D,M)
1524: {
1525: L = [[1,1]]; P0 = P1 = 1;
1526: while ( 1 ) {
1527: P0 *= V;
1528: P1 = sremm(sremm(P*P1,F,M,V),D,M);
1529: L1 = lnf_mod(P0,P1,L,V,D,M); NP0 = L1[0]; NP1 = L1[1];
1530: if ( !NP1 )
1531: return NP0;
1532: else
1533: L = lnf_insert([NP0,NP1],L,V);
1534: }
1535: }
1536:
1537: def lnf_mod(P0,P1,L,V,D,M)
1538: {
1539: NP0 = P0; NP1 = P1;
1540: for ( T = L; T != []; T = cdr(T) ) {
1541: Q = car(T);
1542: D1 = deg(NP1,V);
1543: if ( D1 == deg(Q[1],V) ) {
1544: C = coef(Q[1],D1,V);
1545: INV = inva_mod(D,M,M-C); H = sremm(coef(NP1,D1,V)*INV,D,M);
1546: NP0 = sremm(NP0+Q[0]*H,D,M);
1547: NP1 = sremm(NP1+Q[1]*H,D,M);
1548: }
1549: }
1550: return [NP0,NP1];
1551: }
1552:
1553: def lnf_insert(P,L,V)
1554: {
1555: if ( L == [] )
1556: return [P];
1557: else {
1558: P0 = car(L);
1559: if ( deg(P0[1],V) > deg(P[1],V) )
1560: return cons(P0,lnf_insert(P,cdr(L),V));
1561: else
1562: return cons(P,L);
1563: }
1564: }
1565:
1566: def find_root(P,V,D,M)
1567: {
1568: L = c_z(P,V,1,D,M);
1569: for ( T = L, U = []; T != []; T = cdr(T) ) {
1570: S = monic_mod(car(T),V,D,M); U = cons(-coef(S,0,V),U);
1571: }
1572: return U;
1573: }
1574:
1575: def c_z(F,V,E,D,M)
1576: {
1577: N = deg(F,V);
1578: if ( N == E )
1579: return [F];
1580: Q = M^deg(D,var(D));
1581: K = idiv(N,E);
1582: L = [F];
1583: while ( 1 ) {
1584: W = mrandomgfpoly(2*E,V,D,M);
1585: if ( M == 2 ) {
1586: W = monic_mod(tr_mod(W,F,V,D,M,N-1),V,D,M);
1587: } else {
1588: /* W = monic_mod(pwr_p_mod(W,F,V,D,M,idiv(Q^E-1,2))-1,V,D,M); */
1589: /* T = pwr_p_mod(W,F,V,D,M,idiv(Q^E-1,2)); */
1590: T = pwr_mod(W,F,V,D,M,idiv(Q^E-1,2));
1591: W = monic_mod(T-1,V,D,M);
1592: }
1593: if ( !W )
1594: continue;
1595: G = ag_mod_single4(F,W,[D],M);
1596: if ( deg(G,V) && deg(G,V) < N ) {
1597: L1 = c_z(G,V,E,D,M);
1598: L2 = c_z(sremm(sdivm(F,G,M,V),D,M),V,E,D,M);
1599: return append(L1,L2);
1600: }
1601: }
1602: }
1603:
1604: def tr_mod(P,F,V,D,M,N)
1605: {
1606: for ( I = 1, S = P, W = P; I <= N; I++ ) {
1607: W = sremm(sremm(W^2,F,M,V),D,M);
1608: S = sremm(S+W,D,M);
1609: }
1610: return S;
1611: }
1612:
1613: def mrandomgfpoly(N,V,D,M)
1614: {
1615: W = var(D); ND = deg(D,W);
1616: for ( I = N-2, S = V^(N-1); I >= 0; I-- )
1617: S += randompoly(ND,W,M)*V^I;
1618: return S;
1619: }
1620:
1621: def randompoly(N,V,M)
1622: {
1623: for ( I = 0, S = 0; I < N; I++ )
1624: S += (random()%M)*V^I;
1625: return S;
1626: }
1627:
1628: def mrandompoly(N,V,M)
1629: {
1630: for ( I = N-1, S = V^N; I >=0; I-- )
1631: S += (random()%M)*V^I;
1632: return S;
1633: }
1634:
1635: def srem_by_nf(P,B,V,O) {
1636: dp_ord(O); DP = dp_ptod(P,V);
1637: N = length(B); DB = newvect(N);
1638: for ( I = N-1, IL = []; I >= 0; I-- ) {
1639: DB[I] = dp_ptod(B[I],V);
1640: IL = cons(I,IL);
1641: }
1642: L = dp_true_nf(IL,DP,DB,1);
1643: return [dp_dtop(L[0],V),L[1]];
1644: }
1645: end$