![[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 / asir2018 / lib|
Return to sp CVS log | Up to [local] / OpenXM_contrib2 / asir2018 / lib |
1.1 ! noro 1: /*
! 2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
! 3: * All rights reserved.
! 4: *
! 5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
! 6: * non-exclusive and royalty-free license to use, copy, modify and
! 7: * redistribute, solely for non-commercial and non-profit purposes, the
! 8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
! 9: * conditions of this Agreement. For the avoidance of doubt, you acquire
! 10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
! 11: * third party developer retains all rights, including but not limited to
! 12: * copyrights, in and to the SOFTWARE.
! 13: *
! 14: * (1) FLL does not grant you a license in any way for commercial
! 15: * purposes. You may use the SOFTWARE only for non-commercial and
! 16: * non-profit purposes only, such as academic, research and internal
! 17: * business use.
! 18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
! 19: * international copyright treaties. If you make copies of the SOFTWARE,
! 20: * with or without modification, as permitted hereunder, you shall affix
! 21: * to all such copies of the SOFTWARE the above copyright notice.
! 22: * (3) An explicit reference to this SOFTWARE and its copyright owner
! 23: * shall be made on your publication or presentation in any form of the
! 24: * results obtained by use of the SOFTWARE.
! 25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
! 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
! 27: * for such modification or the source code of the modified part of the
! 28: * SOFTWARE.
! 29: *
! 30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
! 31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
! 32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
! 33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
! 34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
! 35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
! 36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
! 37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
! 38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
! 39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
! 40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
! 41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
! 42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
! 43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
! 44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
! 45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
! 46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
! 47: *
! 48: * $OpenXM$
! 49: */
! 50: /*
! 51: sp : functions related to algebraic number fields
! 52:
! 53: Revision History:
! 54:
! 55: 2001/10/12 noro if USE_PARI_FACTOR is nonzero, pari factor is called
! 56: 2000/03/10 noro fixed several bugs around gathering algebraic numbers
! 57: 1999/08/24 noro modified for 1999 release version
! 58: */
! 59:
! 60: #include "defs.h"
! 61:
! 62: extern ASCENT,GCDTIME,UFTIME,RESTIME,SQTIME,PRINT$
! 63: extern SpOrd$
! 64: extern USE_PARI_FACTOR$
! 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: }
! 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: }
! 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;
! 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;
! 214: V = var(P);
! 215: LC = coef(P,deg(P,V),V);
! 216: if ( ntype(LC) != 1 )
! 217: P = simpalg(1/LC*P);
! 218: S = reverse(asq(P));
! 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];
! 234: V = var(P);
! 235: LC = coef(P,deg(P,V),V);
! 236: if ( ntype(LC) != 1 )
! 237: P = simpalg(1/LC*P);
! 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];
! 262: DCSQ = sortfs(asq(N));
! 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) )
! 280: U = cr_gcda(S,G);
! 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) )
! 408: T = union_sort(T,getalgp(C));
! 409: return T;
! 410: }
! 411: }
! 412:
! 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: }
! 422: }
! 423:
! 424: /* C = union of A and B; A and B is sorted. C should also be sorted. */
! 425:
! 426: def union_sort(A,B)
! 427: {
! 428: if ( A == [] )
! 429: return B;
! 430: else if ( B == [] )
! 431: return A;
! 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: }
! 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: {
! 511: A = union_sort(getalgp(P1),getalgp(P2));
! 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 {
! 523: A = union_sort(getalgp(P1),getalgp(P2));
! 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:
! 562: def asq(P)
! 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);
! 575: GCD = cr_gcda(F,F1);
! 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) {
! 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(factpol,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) {
! 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) )
! 646: PPP = cr_gcda(PPP,P);
! 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);
! 656: HH = cr_gcda(PPP,H);
! 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);
! 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) */
! 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:
! 821: def cr_gcda(P1,P2)
! 822: {
! 823: if ( !P1 )
! 824: return P2;
! 825: if ( !P2 )
! 826: return P1;
! 827: if ( !var(P1) || !var(P2) )
! 828: return 1;
! 829: V = var(P1);
! 830: EXT = union_sort(getalgtreep(P1),getalgtreep(P2));
! 831: if ( EXT == [] )
! 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;
! 857: SpOrd = 2; NOSUGAR = 1;
! 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]);
! 1071: G = dp_gr_mod_main(B,VL,0,MOD,SpOrd);
! 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: {
! 1305: SpOrd = 2;
! 1306: dp_gr_flags(["NoSugar",1]);
! 1307: G = dp_gr_mod_main(cons(x*C-1,D),cons(x,vars(D)),0,MOD,SpOrd);
! 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);
! 1366: LC = coef(N,deg(N,V0),V0);
! 1367: LCD = coef(DN,deg(DN,V0),V0);
! 1368: J = 2;
! 1369: while ( 1 ) {
! 1370: Len = deg(N,V0)+1;
! 1371: for ( I = 0; I < 5; J++ ) {
! 1372: M = prime(J);
! 1373: if ( !(LC%M) || !(LCD%M))
! 1374: continue;
! 1375: G = gcd(N,DN,M);
! 1376: if ( !deg(G,V0) ) {
! 1377: I++;
! 1378: T = nfctr_mod(N,M);
! 1379: if ( T < Len ) {
! 1380: Len = T; M0 = M;
! 1381: }
! 1382: }
! 1383: }
! 1384: S = spm(L,V,M0);
! 1385: if ( S ) break;
! 1386: }
! 1387: T = resfctr_mod(F,S,M0);
! 1388: return [T,S,M0];
! 1389: }
! 1390:
! 1391: def resfctr_mod(F,L,M)
! 1392: {
! 1393: for ( T = L, R = []; T != []; T = cdr(T) ) {
! 1394: U = car(T); MP = U[0]; W = U[1];
! 1395: for ( A = W, B = F; A != []; A = cdr(cdr(A)) )
! 1396: B = sremm(subst(B,A[0],A[1]),MP,M);
! 1397: C = res(var(MP),B,MP) % M;
! 1398: R = cons(flatten(cdr(modfctr(C,M))),R);
! 1399: }
! 1400: return reverse(R);
! 1401: }
! 1402:
! 1403: def flatten(L)
! 1404: {
! 1405: for ( T = L, R = []; T != []; T = cdr(T) )
! 1406: R = cons(car(car(T)),R);
! 1407: return R;
! 1408: }
! 1409:
! 1410: def spm(L,V,M)
! 1411: {
! 1412: if ( length(V) == 1 ) {
! 1413: U = modfctr(car(L),M);
! 1414: for ( T = cdr(U), R = []; T != []; T = cdr(T) ) {
! 1415: S = car(T);
! 1416: if ( S[1] > 1 ) return 0;
! 1417: R = cons([subst(S[0],var(S[0]),a_),[var(S[0]),a_]],R);
! 1418: }
! 1419: return R;
! 1420: }
! 1421: L1 = spm(cdr(L),cdr(V),M);
! 1422: if ( !L1 ) return 0;
! 1423: F0 = car(L); V0 = car(V); VR = cdr(V);
! 1424: for ( T = L1, R = []; T != []; T = cdr(T) ) {
! 1425: S = car(T);
! 1426: F1 = subst(F0,S[1]);
! 1427: U = fctr_mod(F1,V0,S[0],M);
! 1428: VS = var(S[0]);
! 1429: for ( W = U; W != []; W = cdr(W) ) {
! 1430: if ( car(W)[1] > 1 ) return 0;
! 1431: A = car(car(W));
! 1432: if ( deg(A,V0) == 1 ) {
! 1433: A = monic_mod(A,V0,S[0],M);
! 1434: R = cons([S[0],append([V0,-coef(A,0,V0)],S[1])],R);
! 1435: } else {
! 1436: B = pe_mod(A,S[0],M);
! 1437: MP = B[0]; VMP = var(MP); NV = B[1];
! 1438: for ( C = S[1], D = []; C != []; C = cdr(cdr(C)) ) {
! 1439: G = subst(sremm(subst(C[1],VS,NV[1]),MP,M),VMP,VS);
! 1440: D = append([C[0],G],D);
! 1441: }
! 1442: R = cons([subst(MP,VMP,VS),
! 1443: append([B[2][0],subst(B[2][1],VMP,VS)],D)],R);
! 1444: }
! 1445: }
! 1446: }
! 1447: return R;
! 1448: }
! 1449:
! 1450: def pe_mod(F,G,M)
! 1451: {
! 1452: V = var(G); W = car(setminus(vars(F),[V]));
! 1453: NG = deg(G,V); NF = deg(F,W); N = NG*NF;
! 1454: X = prim;
! 1455: while ( 1 ) {
! 1456: D = mrandompoly(N,X,M);
! 1457: if ( irred_check(D,M) )
! 1458: break;
! 1459: }
! 1460: L = fctr_mod(G,V,D,M);
! 1461: for ( T = L; T != []; T = cdr(T) ) {
! 1462: U = car(car(T));
! 1463: if ( deg(U,V) == 1 )
! 1464: break;
! 1465: }
! 1466: U = monic_mod(U,V,D,M); RV = -coef(U,0,V);
! 1467: L = fctr_mod(sremm(subst(F,V,RV),D,M),W,D,M);
! 1468: for ( T = L; T != []; T = cdr(T) ) {
! 1469: U = car(car(T));
! 1470: if ( deg(U,W) == 1 )
! 1471: break;
! 1472: }
! 1473: U = monic_mod(U,W,D,M); RW = -coef(U,0,W);
! 1474: return [D,[V,RV],[W,RW]];
! 1475: }
! 1476:
! 1477: def fctr_mod(F,V,D,M)
! 1478: {
! 1479: if ( V != x ) {
! 1480: F = subst(F,V,x); V0 = V; V = x;
! 1481: } else
! 1482: V0 = x;
! 1483: F = monic_mod(F,V,D,M);
! 1484: L = sqfr_mod(F,V,D,M);
! 1485: for ( R = [], T = L; T != []; T = cdr(T) ) {
! 1486: S = car(T); A = S[0]; E = S[1];
! 1487: B = ddd_mod(A,V,D,M);
! 1488: R = append(append_mult(B,E),R);
! 1489: }
! 1490: if ( V0 != x ) {
! 1491: for ( R = reverse(R), T = []; R != []; R = cdr(R) )
! 1492: T = cons([subst(car(R)[0],x,V0),car(R)[1]],T);
! 1493: R = T;
! 1494: }
! 1495: return R;
! 1496: }
! 1497:
! 1498: def append_mult(L,E)
! 1499: {
! 1500: for ( T = L, R = []; T != []; T = cdr(T) )
! 1501: R = cons([car(T),E],R);
! 1502: return R;
! 1503: }
! 1504:
! 1505: def sqfr_mod(F,V,D,M)
! 1506: {
! 1507: setmod(M);
! 1508: F = sremm(F,D,M);
! 1509: F1 = sremm(diff(F,V),D,M);
! 1510: F1 = sremm(F1*inva_mod(D,M,LCOEF(F1)),D,M);
! 1511: if ( F1 ) {
! 1512: F2 = ag_mod_single4(F,F1,[D],M);
! 1513: FLAT = sremm(sdivm(F,F2,M,V),D,M);
! 1514: I = 0; L = [];
! 1515: while ( deg(FLAT,V) ) {
! 1516: while ( 1 ) {
! 1517: QR = sqrm(F,FLAT,M,V);
! 1518: if ( !sremm(QR[1],D,M) ) {
! 1519: F = sremm(QR[0],D,M); I++;
! 1520: } else
! 1521: break;
! 1522: }
! 1523: if ( !deg(F,V) )
! 1524: FLAT1 = 1;
! 1525: else
! 1526: FLAT1 = ag_mod_single4(F,FLAT,[D],M);
! 1527: G = sremm(sdivm(FLAT,FLAT1,M,V),D,M);
! 1528: FLAT = FLAT1;
! 1529: L = cons([G,I],L);
! 1530: }
! 1531: }
! 1532: if ( deg(F,V) ) {
! 1533: T = sqfr_mod(pthroot_p_mod(F,V,D,M),V,D,M);
! 1534: for ( R = []; T != []; T = cdr(T) ) {
! 1535: H = car(T); R = cons([H[0],M*H[1]],R);
! 1536: }
! 1537: } else
! 1538: R = [];
! 1539: return append(L,R);
! 1540: }
! 1541:
! 1542: def pthroot_p_mod(F,V,D,M)
! 1543: {
! 1544: for ( T = F, R = 0; T; ) {
! 1545: D1 = deg(T,V); C = coef(T,D1,V); T -= C*V^D1;
! 1546: R += pthroot_n_mod(C,D,M)*V^idiv(D1,M);
! 1547: }
! 1548: return R;
! 1549: }
! 1550:
! 1551: def pthroot_n_mod(C,D,M)
! 1552: {
! 1553: pwr_n_mod(C,D,M,deg(D,var(D))-1);
! 1554: }
! 1555:
! 1556: def pwr_n_mod(C,D,M,N)
! 1557: {
! 1558: if ( N == 0 )
! 1559: return 1;
! 1560: else if ( N == 1 )
! 1561: return C;
! 1562: else {
! 1563: QR = iqr(N,2);
! 1564: T = pwr_n_mod(C,D,M,QR[0]);
! 1565: S = sremm(T^2,D,M);
! 1566: if ( QR[1] )
! 1567: return sremm(S*C,D,M);
! 1568: else
! 1569: return S;
! 1570: }
! 1571: }
! 1572:
! 1573: def pwr_p_mod(P,A,V,D,M,N)
! 1574: {
! 1575: if ( N == 0 )
! 1576: return 1;
! 1577: else if ( N == 1 )
! 1578: return P;
! 1579: else {
! 1580: QR = iqr(N,2);
! 1581: T = pwr_p_mod(P,A,V,D,M,QR[0]);
! 1582: S = sremm(sremm(sremm(T^2,D,M),A,M,V),D,M);
! 1583: if ( QR[1] )
! 1584: return sremm(sremm(sremm(S*P,D,M),A,M,V),D,M);
! 1585: else
! 1586: return S;
! 1587: }
! 1588: }
! 1589:
! 1590: def qmat_mod(F,V,D,M)
! 1591: {
! 1592: R = tab_mod(F,V,D,M);
! 1593: Q = newmat(N,N);
! 1594: for ( J = 0; J < N; J++ )
! 1595: for ( I = 0, T = R[J]; I < N; I++ ) {
! 1596: Q[I][J] = coef(T,I);
! 1597: }
! 1598: for ( I = 0; I < N; I++ )
! 1599: Q[I][I] = (Q[I][I]+(M-1))%M;
! 1600: return Q;
! 1601: }
! 1602:
! 1603: def tab_mod(F,V,D,M)
! 1604: {
! 1605: MD = M^deg(D,var(D));
! 1606: N = deg(F,V);
! 1607: F = sremm(F*inva_mod(D,M,coef(F,N,V)),D,M);
! 1608: R = newvect(N); R[0] = 1;
! 1609: R[1] = pwr_mod(V,F,V,D,M,MD);
! 1610: for ( I = 2; I < N; I++ )
! 1611: R[I] = sremm(sremm(R[1]*R[I-1],F,M),D,M);
! 1612: return R;
! 1613: }
! 1614:
! 1615: def ddd_mod(F,V,D,M)
! 1616: {
! 1617: if ( deg(F,V) == 1 )
! 1618: return [F];
! 1619: TAB = tab_mod(F,V,D,M);
! 1620: for ( I = 1, W = V, L = []; 2*I <= deg(F,V); I++ ) {
! 1621: for ( T = 0, K = 0; K <= deg(W,V); K++ )
! 1622: if ( C = coef(W,K,V) )
! 1623: T = sremm(T+TAB[K]*C,D,M);
! 1624: W = T;
! 1625: GCD = ag_mod_single4(F,monic_mod(W-V,V,D,M),[D],M);
! 1626: if ( deg(GCD,V) ) {
! 1627: L = append(berlekamp(GCD,V,I,TAB,D,M),L);
! 1628: F = sremm(sdivm(F,GCD,M,V),D,M);
! 1629: W = sremm(sremm(W,F,M,V),D,M);
! 1630: }
! 1631: }
! 1632: if ( deg(F,V) )
! 1633: return cons(F,L);
! 1634: else
! 1635: return L;
! 1636: }
! 1637:
! 1638: def monic_mod(F,V,D,M) {
! 1639: if ( !F || !deg(F,V) )
! 1640: return F;
! 1641: return sremm(F*inva_mod(D,M,coef(F,deg(F,V),V)),D,M);
! 1642: }
! 1643:
! 1644: def berlekamp(F,V,E,TAB,D,M)
! 1645: {
! 1646: N = deg(F,V);
! 1647: Q = newmat(N,N);
! 1648: for ( J = 0; J < N; J++ ) {
! 1649: T = sremm(sremm(TAB[J],F,M,V),D,M);
! 1650: for ( I = 0; I < N; I++ ) {
! 1651: Q[I][J] = coef(T,I);
! 1652: }
! 1653: }
! 1654: for ( I = 0; I < N; I++ )
! 1655: Q[I][I] = (Q[I][I]+(M-1))%M;
! 1656: L = nullspace(Q,D,M); MT = L[0]; IND = L[1];
! 1657: NF0 = N/E;
! 1658: PS = null_to_poly(MT,IND,V,M);
! 1659: R = newvect(NF0); R[0] = monic_mod(F,V,D,M);
! 1660: for ( I = 1, NF = 1; NF < NF0 && I < NF0; I++ ) {
! 1661: PSI = PS[I];
! 1662: MP = minipoly_mod(PSI,F,V,D,M);
! 1663: ROOT = find_root(MP,V,D,M); NR = length(ROOT);
! 1664: for ( J = 0; J < NF; J++ ) {
! 1665: if ( deg(R[J],V) == E )
! 1666: continue;
! 1667: for ( K = 0; K < NR; K++ ) {
! 1668: GCD = ag_mod_single4(R[J],PSI-ROOT[K],[D],M);
! 1669: if ( deg(GCD,V) > 0 && deg(GCD,V) < deg(R[J],V) ) {
! 1670: Q = sremm(sdivm(R[J],GCD,M,V),D,M);
! 1671: R[J] = Q; R[NF++] = GCD;
! 1672: }
! 1673: }
! 1674: }
! 1675: }
! 1676: return vtol(R);
! 1677: }
! 1678:
! 1679: def null_to_poly(MT,IND,V,M)
! 1680: {
! 1681: N = size(MT)[0];
! 1682: for ( I = 0, J = 0; I < N; I++ )
! 1683: if ( IND[I] )
! 1684: J++;
! 1685: R = newvect(J);
! 1686: for ( I = 0, L = 0; I < N; I++ ) {
! 1687: if ( !IND[I] )
! 1688: continue;
! 1689: for ( J = K = 0, T = 0; J < N; J++ )
! 1690: if ( !IND[J] )
! 1691: T += MT[K++][I]*V^J;
! 1692: else if ( J == I )
! 1693: T += (M-1)*V^I;
! 1694: R[L++] = T;
! 1695: }
! 1696: return R;
! 1697: }
! 1698:
! 1699: def minipoly_mod(P,F,V,D,M)
! 1700: {
! 1701: L = [[1,1]]; P0 = P1 = 1;
! 1702: while ( 1 ) {
! 1703: P0 *= V;
! 1704: P1 = sremm(sremm(P*P1,F,M,V),D,M);
! 1705: L1 = lnf_mod(P0,P1,L,V,D,M); NP0 = L1[0]; NP1 = L1[1];
! 1706: if ( !NP1 )
! 1707: return NP0;
! 1708: else
! 1709: L = lnf_insert([NP0,NP1],L,V);
! 1710: }
! 1711: }
! 1712:
! 1713: def lnf_mod(P0,P1,L,V,D,M)
! 1714: {
! 1715: NP0 = P0; NP1 = P1;
! 1716: for ( T = L; T != []; T = cdr(T) ) {
! 1717: Q = car(T);
! 1718: D1 = deg(NP1,V);
! 1719: if ( D1 == deg(Q[1],V) ) {
! 1720: C = coef(Q[1],D1,V);
! 1721: INV = inva_mod(D,M,M-C); H = sremm(coef(NP1,D1,V)*INV,D,M);
! 1722: NP0 = sremm(NP0+Q[0]*H,D,M);
! 1723: NP1 = sremm(NP1+Q[1]*H,D,M);
! 1724: }
! 1725: }
! 1726: return [NP0,NP1];
! 1727: }
! 1728:
! 1729: def lnf_insert(P,L,V)
! 1730: {
! 1731: if ( L == [] )
! 1732: return [P];
! 1733: else {
! 1734: P0 = car(L);
! 1735: if ( deg(P0[1],V) > deg(P[1],V) )
! 1736: return cons(P0,lnf_insert(P,cdr(L),V));
! 1737: else
! 1738: return cons(P,L);
! 1739: }
! 1740: }
! 1741:
! 1742: def find_root(P,V,D,M)
! 1743: {
! 1744: L = c_z(P,V,1,D,M);
! 1745: for ( T = L, U = []; T != []; T = cdr(T) ) {
! 1746: S = monic_mod(car(T),V,D,M); U = cons(-coef(S,0,V),U);
! 1747: }
! 1748: return U;
! 1749: }
! 1750:
! 1751: def c_z(F,V,E,D,M)
! 1752: {
! 1753: N = deg(F,V);
! 1754: if ( N == E )
! 1755: return [F];
! 1756: Q = M^deg(D,var(D));
! 1757: K = idiv(N,E);
! 1758: L = [F];
! 1759: while ( 1 ) {
! 1760: W = mrandomgfpoly(2*E,V,D,M);
! 1761: if ( M == 2 ) {
! 1762: W = monic_mod(tr_mod(W,F,V,D,M,N-1),V,D,M);
! 1763: } else {
! 1764: /* W = monic_mod(pwr_p_mod(W,F,V,D,M,idiv(Q^E-1,2))-1,V,D,M); */
! 1765: /* T = pwr_p_mod(W,F,V,D,M,idiv(Q^E-1,2)); */
! 1766: T = pwr_mod(W,F,V,D,M,idiv(Q^E-1,2));
! 1767: W = monic_mod(T-1,V,D,M);
! 1768: }
! 1769: if ( !W )
! 1770: continue;
! 1771: G = ag_mod_single4(F,W,[D],M);
! 1772: if ( deg(G,V) && deg(G,V) < N ) {
! 1773: L1 = c_z(G,V,E,D,M);
! 1774: L2 = c_z(sremm(sdivm(F,G,M,V),D,M),V,E,D,M);
! 1775: return append(L1,L2);
! 1776: }
! 1777: }
! 1778: }
! 1779:
! 1780: def tr_mod(P,F,V,D,M,N)
! 1781: {
! 1782: for ( I = 1, S = P, W = P; I <= N; I++ ) {
! 1783: W = sremm(sremm(W^2,F,M,V),D,M);
! 1784: S = sremm(S+W,D,M);
! 1785: }
! 1786: return S;
! 1787: }
! 1788:
! 1789: def mrandomgfpoly(N,V,D,M)
! 1790: {
! 1791: W = var(D); ND = deg(D,W);
! 1792: for ( I = N-2, S = V^(N-1); I >= 0; I-- )
! 1793: S += randompoly(ND,W,M)*V^I;
! 1794: return S;
! 1795: }
! 1796:
! 1797: def randompoly(N,V,M)
! 1798: {
! 1799: for ( I = 0, S = 0; I < N; I++ )
! 1800: S += (random()%M)*V^I;
! 1801: return S;
! 1802: }
! 1803:
! 1804: def mrandompoly(N,V,M)
! 1805: {
! 1806: for ( I = N-1, S = V^N; I >=0; I-- )
! 1807: S += (random()%M)*V^I;
! 1808: return S;
! 1809: }
! 1810:
! 1811: def srem_by_nf(P,B,V,O) {
! 1812: dp_ord(O); DP = dp_ptod(P,V);
! 1813: N = length(B); DB = newvect(N);
! 1814: for ( I = N-1, IL = []; I >= 0; I-- ) {
! 1815: DB[I] = dp_ptod(B[I],V);
! 1816: IL = cons(I,IL);
! 1817: }
! 1818: L = dp_true_nf(IL,DP,DB,1);
! 1819: return [dp_dtop(L[0],V),L[1]];
! 1820: }
! 1821: end$