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