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