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