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