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