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