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