Annotation of OpenXM/src/asir-contrib/testing/noro/wishartd.rr, Revision 1.1
1.1 ! noro 1: /* A package for computing PDEs satisfied by a matrix 1F1 on diagonal regions */
! 2: if (version()<20151202) {error("The version of Risa/Asir must be after 20151202 to run this package.");} else{}
! 3: module n_wishartd$
! 4: localf compf1$
! 5: localf my_comp_f$
! 6: localf compc1, compc, compd, compt, addpf, subpf, addc, ftorat, ctorat, mulf$
! 7: localf lcmf, mulf, divf, mulc, mulcpf, diffc, muldpf, normalizef1, normalizec1, normalizepf$
! 8: localf wsetup, mtot, wishartpf, degf, removef, subst_f, lopitalpf, simpop, simppf$
! 9: localf eliminatetop, diagpf, diag0pf, adjop1, adjop, reducepf, reduceotherspf$
! 10: localf print_f, printc, printt, printpf, ftop, ctord, comppfrd, multpf, spolypf, nfpf$
! 11: localf nfpf0$
! 12: localf substc$
! 13: localf dpf,dpf2,dpf3,abs$
! 14: localf pftord$
! 15: localf lopital_others_pf$
! 16: localf pftop$
! 17: localf dnc,pftozpf,addzpf,zpftopf$
! 18: localf mul_lopitalpf$
! 19: localf indep_eqs$
! 20: localf mulpf$
! 21: localf pftoeuler$
! 22: localf gammam,prc,prc2$
! 23: localf ps,psvalue,act_op,decomp_op,create_homo,bpsvalue$
! 24: localf pfaffian,gen_basis$
! 25: localf evalpf,eval_pfaffian,genrk45$
! 26: localf solve_leq,ptom$
! 27: localf msubst$
! 28: localf vton,encodef1,encodef,encodec1,encodec$
! 29: localf vton,ftotex,ctotex,pftotex,optotex,eqtotex$
! 30: localf genprc,prob_by_hgm,prob_by_ps$
! 31: localf partition, all_partition, coef_partition, count_dup,simp_power1$
! 32: localf lop1,lopn,loph,mul_lopitalpf2$
! 33: localf muldmpf$
! 34: localf diagcheck,eliminatecheck$
! 35: localf getchi$
! 36: localf message$
! 37: localf pfaffian2, eval_pfaffian2$
! 38:
! 39: /*
! 40: * pfpoly = [[C,<<...>>],...]
! 41: * C = [[N,L],...] L = [[F,M],...]
! 42: * lc(F) = 1
! 43: */
! 44:
! 45: static Print$
! 46: static PF_N,PF_V,PF_DV$
! 47: static Tlopital,Tred,Tlineq,Tactmul$
! 48: static Tp,Tps,Trk$
! 49:
! 50: /* XXX execute ord([y1,y2,...,dy1,dy2,...]) */
! 51:
! 52: /* F1>F2 <=> var(F1)>var(F2) || var(F1)=var(F2) && rest(F1)>rest(F2) */
! 53: /* F1^M1 > F2^M2 <=> F1>F2 || F1=F2 && M1>M2 */
! 54:
! 55: def abs(A) { return A>0 ? A : -A; }
! 56:
! 57: def compf1(E1,E2)
! 58: {
! 59: F1 = E1[0]; F2 = E2[0];
! 60: if ( F1>F2 ) return 1;
! 61: if ( F1<F2 ) return -1;
! 62: M1 = E1[1]; M2 = E2[1];
! 63: if ( M1 > M2 ) return 1;
! 64: if ( M1 < M2 ) return -1;
! 65: return 0;
! 66: }
! 67:
! 68: def compc1(ND1,ND2)
! 69: {
! 70: if ( R = comp_f(ND1[1],ND2[1]) ) return R;
! 71: N1 = ND1[0]; N2 = ND2[0];
! 72: if ( N1 > N2 ) return 1;
! 73: if ( N1 < N2 ) return -1;
! 74: return 0;
! 75: }
! 76:
! 77: /* compare ND lists */
! 78: def compc(L1,L2)
! 79: {
! 80: for ( ; L1 != [] && L2 != []; L1 = cdr(L1), L2 = cdr(L2) ) {
! 81: E1 = car(L1); E2 = car(L2);
! 82: if ( R = compc1(E1,E2) ) return R;
! 83: }
! 84: if ( L1 != [] ) return 1;
! 85: if ( L2 != [] ) return -1;
! 86: return 0;
! 87: }
! 88:
! 89: def compd(D1,D2)
! 90: {
! 91: if ( D1 > D2 ) return 1;
! 92: if ( D1 < D2 ) return -1;
! 93: return 0;
! 94: }
! 95:
! 96: /* compare terms */
! 97: def compt(T1,T2)
! 98: {
! 99: if ( R = compd(T1[1],T2[1]) ) return R;
! 100: if ( R = compc(T1[0],T2[0]) ) return R;
! 101: return 0;
! 102: }
! 103:
! 104: def addpf(F1,F2)
! 105: {
! 106: R = [];
! 107: while ( F1 != [] && F2 != [] ) {
! 108: T1 = car(F1); T2 = car(F2);
! 109: C = compd(T1[1],T2[1]);
! 110: if ( C > 0 ) {
! 111: R = cons(T1,R); F1 = cdr(F1);
! 112: } else if ( C < 0 ) {
! 113: R = cons(T2,R); F2 = cdr(F2);
! 114: } else {
! 115: S = addc(T1[0],T2[0]);
! 116: if ( S != [] )
! 117: R = cons([S,T1[1]],R);
! 118: F1 = cdr(F1); F2 = cdr(F2);
! 119: }
! 120: }
! 121: R = reverse(R);
! 122: if ( F1 != [] ) R = append(R,F1);
! 123: else if ( F2 != [] ) R = append(R,F2);
! 124: return R;
! 125: }
! 126:
! 127: def subpf(F1,F2)
! 128: {
! 129: T = mulcpf([[-1,[]]],F2);
! 130: T = addpf(F1,T);
! 131: return T;
! 132: }
! 133:
! 134: def addc(F1,F2)
! 135: {
! 136: R = [];
! 137: while ( F1 != [] && F2 != [] ) {
! 138: T1 = car(F1); T2 = car(F2);
! 139: C = comp_f(T1[1],T2[1]);
! 140: if ( !T1[0] || !T2[0] )
! 141: error("afo");
! 142: if ( C > 0 ) {
! 143: R = cons(T1,R); F1 = cdr(F1);
! 144: } else if ( C < 0 ) {
! 145: R = cons(T2,R); F2 = cdr(F2);
! 146: } else {
! 147: S = T1[0]+T2[0];
! 148: if ( S )
! 149: R = cons([S,T1[1]],R);
! 150: F1 = cdr(F1); F2 = cdr(F2);
! 151: }
! 152: }
! 153: R = reverse(R);
! 154: if ( F1 != [] ) R = append(R,F1);
! 155: else if ( F2 != [] ) R = append(R,F2);
! 156: return R;
! 157: }
! 158:
! 159: def ftorat(F)
! 160: {
! 161: R = 1;
! 162: for ( ; F != []; F = cdr(F) ) {
! 163: F0 = car(F);
! 164: R *= F0[0]^F0[1];
! 165: }
! 166: return R;
! 167: }
! 168:
! 169: def ctorat(C)
! 170: {
! 171: R = 0;
! 172: for ( ; C != []; C = cdr(C) ) {
! 173: C0 = car(C);
! 174: R += C0[0]/ftorat(C0[1]);
! 175: R = red(R);
! 176: }
! 177: return R;
! 178: }
! 179:
! 180: def mulf(L1,L2)
! 181: {
! 182: R = [];
! 183: while ( L1 != [] && L2 != [] ) {
! 184: A1 = car(L1); A2 = car(L2);
! 185: if ( A1[0] > A2[0] ) {
! 186: R = cons(A1,R); L1 = cdr(L1);
! 187: } else if ( A1[0] < A2[0] ) {
! 188: R = cons(A2,R); L2 = cdr(L2);
! 189: } else {
! 190: R = cons([A1[0],A1[1]+A2[1]],R);
! 191: L1 = cdr(L1); L2 = cdr(L2);
! 192: }
! 193: }
! 194: R = reverse(R);
! 195: if ( L2 == [] ) return append(R,L1);
! 196: else if ( L1 == [] ) return append(R,L2);
! 197: else return R;
! 198: }
! 199:
! 200: def lcmf(L1,L2)
! 201: {
! 202: R = [];
! 203: while ( L1 != [] && L2 != [] ) {
! 204: A1 = car(L1); A2 = car(L2);
! 205: if ( A1[0] > A2[0] ) {
! 206: R = cons(A1,R); L1 = cdr(L1);
! 207: } else if ( A1[0] < A2[0] ) {
! 208: R = cons(A2,R); L2 = cdr(L2);
! 209: } else {
! 210: M = A1[1]>A2[1]?A1[1]:A2[1];
! 211: R = cons([A1[0],M],R);
! 212: L1 = cdr(L1); L2 = cdr(L2);
! 213: }
! 214: }
! 215: R = reverse(R);
! 216: if ( L2 == [] ) return append(R,L1);
! 217: else if ( L1 == [] ) return append(R,L2);
! 218: else return R;
! 219: }
! 220:
! 221: def mulf(L1,L2)
! 222: {
! 223: R = [];
! 224: while ( L1 != [] && L2 != [] ) {
! 225: A1 = car(L1); A2 = car(L2);
! 226: if ( A1[0] > A2[0] ) {
! 227: R = cons(A1,R); L1 = cdr(L1);
! 228: } else if ( A1[0] < A2[0] ) {
! 229: R = cons(A2,R); L2 = cdr(L2);
! 230: } else {
! 231: R = cons([A1[0],A1[1]+A2[1]],R);
! 232: L1 = cdr(L1); L2 = cdr(L2);
! 233: }
! 234: }
! 235: R = reverse(R);
! 236: if ( L2 == [] ) return append(R,L1);
! 237: else if ( L1 == [] ) return append(R,L2);
! 238: else return R;
! 239: }
! 240:
! 241: def divf(L1,L2)
! 242: {
! 243: R = [];
! 244: while ( L1 != [] && L2 != [] ) {
! 245: A1 = car(L1); A2 = car(L2);
! 246: if ( A1[0] > A2[0] ) {
! 247: R = cons(A1,R); L1 = cdr(L1);
! 248: } else if ( A1[0] < A2[0] ) {
! 249: error("divf : cannot happen");
! 250: } else {
! 251: M = A1[1]-A2[1];
! 252: if ( M < 0 )
! 253: error("divf : cannot happen");
! 254: if ( M > 0 )
! 255: R = cons([A1[0],M],R);
! 256: L1 = cdr(L1); L2 = cdr(L2);
! 257: }
! 258: }
! 259: R = reverse(R);
! 260: if ( L2 == [] ) return append(R,L1);
! 261: else if ( L1 == [] ) return append(R,L2);
! 262: else return R;
! 263: }
! 264:
! 265: def mulc(C1,C2)
! 266: {
! 267: R = [];
! 268: C1 = reverse(C1);
! 269: for ( ; C1 != []; C1 = cdr(C1) ) {
! 270: S = [];
! 271: A1 = car(C1);
! 272: for ( T = C2 ; T != []; T = cdr(T) ) {
! 273: A2 = car(T);
! 274: S = cons([A1[0]*A2[0],mulf(A1[1],A2[1])],S);
! 275: }
! 276: S = reverse(S);
! 277: R = addc(R,S);
! 278: }
! 279: return R;
! 280: }
! 281:
! 282: def vton(V)
! 283: {
! 284: for ( I = 1; I <= PF_N; I++ )
! 285: if ( V == PF_V[I] ) return I;
! 286: error("vton : no such variable");
! 287: }
! 288:
! 289: def encodef1(F)
! 290: {
! 291: A = F[0];
! 292: V1 = var(A);
! 293: N1 = vton(V1);
! 294: R = A-V1;
! 295: if ( !R )
! 296: return [N1,F[1]];
! 297: else
! 298: return [N1*65536+vton(var(R)),F[1]];
! 299: }
! 300:
! 301: def encodef(F)
! 302: {
! 303: return map(encodef1,F);
! 304: }
! 305:
! 306: def encodec1(C)
! 307: {
! 308: return cons(C[0],encodef(C[1]));
! 309: }
! 310:
! 311: def encodec(C)
! 312: {
! 313: R = map(encodec1,C);
! 314: }
! 315:
! 316: def mulcpf(C,F)
! 317: {
! 318: R = [];
! 319: for ( ; F != []; F = cdr(F) ) {
! 320: F0 = car(F);
! 321: C1 = mulc(C,F0[0]);
! 322: R = cons([C1,F0[1]],R);
! 323: }
! 324: return reverse(R);
! 325: }
! 326:
! 327: def mulpf(F1,F2)
! 328: {
! 329: R = [];
! 330: N = length(PF_V);
! 331: for ( T = reverse(F1); T != []; T = cdr(T) ) {
! 332: T0 = car(T); C = T0[0]; Op = T0[1];
! 333: E = dp_etov(Op);
! 334: S = F2;
! 335: for ( I = 0; I < N; I++ )
! 336: for ( J = 0; J < E[I]; J++ ) S = muldpf(PF_V[I],S);
! 337: S = mulcpf(C,S);
! 338: R = addpf(R,S);
! 339: }
! 340: return S;
! 341: }
! 342:
! 343: def diffc(X,C)
! 344: {
! 345: R = [];
! 346: for ( T = C; T != []; T = cdr(T) ) {
! 347: T0 = car(T);
! 348: N = T0[0];
! 349: S = [];
! 350: for ( L = T0[1]; L != []; S = cons(car(L),S), L = cdr(L) ) {
! 351: L0 = car(L);
! 352: F = L0[0]; M = L0[1];
! 353: DF = diff(F,X);
! 354: if ( DF ) {
! 355: V = cons([F,M+1],cdr(L));
! 356: for ( U = S; U != []; U = cdr(U) ) V = cons(car(U),V);
! 357: R = addc([[-M*N*DF,V]],R);
! 358: }
! 359: }
! 360: }
! 361: return R;
! 362: }
! 363:
! 364: def muldpf(X,F)
! 365: {
! 366: R = [];
! 367: DX = dp_ptod(strtov("d"+rtostr(X)),PF_DV);
! 368: for ( T = F; T != []; T = cdr(T) ) {
! 369: T0 = car(T);
! 370: C = diffc(X,T0[0]);
! 371: if ( C != [] )
! 372: R = addpf(R,[[T0[0],T0[1]*DX],[C,T0[1]]]);
! 373: else
! 374: R = addpf(R,[[T0[0],T0[1]*DX]]);
! 375: }
! 376: return R;
! 377: }
! 378:
! 379: /* d^m*c = sum_{i=0}^m m!/i!/(m-i)!*c^(i)*d^(m-i) */
! 380: def muldmpf(X,M,F)
! 381: {
! 382: R = [];
! 383: DX = strtov("d"+rtostr(X));
! 384: FM = fac(M);
! 385: for ( T = F; T != []; T = cdr(T) ) {
! 386: T0 = car(T);
! 387: C = T0[0]; Op = T0[1];
! 388: U = [];
! 389: for ( I = 0; I <= M; I++ ) {
! 390: if ( C == [] ) break;
! 391: C0 = [[FM/fac(I)/fac(M-I),[]]];
! 392: U = cons([mulc(C0,C),dp_ptod(DX^(M-I),PF_DV)*Op],U);
! 393: C = diffc(X,C);
! 394: }
! 395: U = reverse(U);
! 396: R = addpf(U,R);
! 397: }
! 398: return R;
! 399: }
! 400:
! 401: def normalizef1(F)
! 402: {
! 403: if ( coef(F,1,var(F)) < 0 ) F = -F;
! 404: return F;
! 405: }
! 406:
! 407: def normalizec1(C)
! 408: {
! 409: N = C[0]; L = C[1];
! 410: S = 1;
! 411: R = [];
! 412: for ( ; L != []; L = cdr(L) ) {
! 413: L0 = car(L);
! 414: F = L0[0]; M = L0[1];
! 415: if ( coef(F,1,var(F)) < 0 ) {
! 416: F = -F;
! 417: if ( M%2 ) S = -S;
! 418: }
! 419: R = cons([F,M],R);
! 420: }
! 421: return [S*N,reverse(qsort(R,n_wishartd.compf1))];
! 422: }
! 423:
! 424: def normalizepf(F)
! 425: {
! 426: R = [];
! 427: for ( ; F != []; F = cdr(F) ) {
! 428: F0 = car(F);
! 429: C = map(normalizec1,F[0]);
! 430: R = cons([C,F[1]],R);
! 431: }
! 432: return reverse(R);
! 433: }
! 434:
! 435: def substc(C,Y2,Y1)
! 436: {
! 437: R = [];
! 438: for ( T = C; T != []; T = cdr(T) ) {
! 439: T0 = car(T); N = T0[0]; D = T0[1];
! 440: Z = subst_f(D,Y2,Y1);
! 441: R = addc([[Z[0]*N,Z[1]]],R);
! 442: }
! 443: return R;
! 444: }
! 445:
! 446: /* dy0 : dummy variable for dy */
! 447:
! 448: def wsetup(N)
! 449: {
! 450: V = [];
! 451: DV = [];
! 452: for ( I = N; I>= 0; I-- ) {
! 453: YI = strtov("y"+rtostr(I));
! 454: DYI = strtov("dy"+rtostr(I));
! 455: V = cons(YI,V);
! 456: DV = cons(DYI,DV);
! 457: }
! 458: PF_N = N;
! 459: PF_V = V;
! 460: PF_DV = DV;
! 461: ord(append(V,DV));
! 462: }
! 463:
! 464: def mtot(M,Dn)
! 465: {
! 466: D = dp_ptod(M,PF_DV);
! 467: F = fctr(Dn); C = car(F)[0];
! 468: Dn = reverse(qsort(cdr(F),n_wishartd.compf1));
! 469: return [[[dp_hc(D)/C,Dn]],dp_ht(D)];
! 470: }
! 471:
! 472: def wishartpf(N,I)
! 473: {
! 474: YI = PF_V[I]; DYI = PF_DV[I];
! 475: R = [mtot(DYI^2,1)];
! 476: R = addpf([mtot(c*DYI,YI)],R);
! 477: R = addpf([mtot(-DYI,1)],R);
! 478: for ( J = 1; J <= N; J++ ) {
! 479: if ( J == I ) continue;
! 480: YJ = PF_V[J]; DYJ = PF_DV[J];
! 481: R = addpf([mtot(1/2*DYI,YI-YJ)],R);
! 482: R = addpf([mtot(-1/2*DYJ,YI-YJ)],R);
! 483: R = addpf([mtot(-1/2*DYI,YI)],R);
! 484: R = addpf([mtot(1/2*DYJ,YI)],R);
! 485: }
! 486: R = addpf([mtot(-a,YI)],R);
! 487: return R;
! 488: }
! 489:
! 490: def degf(F,D)
! 491: {
! 492: for ( ; F != []; F = cdr(F) ) {
! 493: F0 = car(F);
! 494: if ( F0[0] == D ) return F0[1];
! 495: }
! 496: return 0;
! 497: }
! 498:
! 499: def removef(F,D)
! 500: {
! 501: R = [];
! 502: for ( ; F != []; F = cdr(F) ) {
! 503: F0 = car(F);
! 504: if ( F0[0] != D ) R = cons(F0,R);
! 505: }
! 506: return reverse(R);
! 507: }
! 508:
! 509: def subst_f(F,Y2,Y1)
! 510: {
! 511: R = [];
! 512: Sgn = 1;
! 513: for ( ; F != []; F = cdr(F) ) {
! 514: F0 = car(F);
! 515: T = subst(F0[0],Y2,Y1);
! 516: if ( coef(T,1,var(T)) < 0 ) {
! 517: T = -T;
! 518: if ( F0[1]%2 ) Sgn = -Sgn;
! 519: }
! 520: R = cons([T,F0[1]],R);
! 521: }
! 522: if ( R == [] ) return [Sgn,R];
! 523: R = qsort(R,n_wishartd.compf1);
! 524: S = [car(R)]; R = cdr(R);
! 525: while ( R != [] ) {
! 526: R0 = car(R);
! 527: S0 = car(S);
! 528: if ( R0[0] == S0[0] )
! 529: S = cons([S0[0],S0[1]+R0[1]],cdr(S));
! 530: else
! 531: S = cons(R0,S);
! 532: R = cdr(R);
! 533: }
! 534: return [Sgn,S];
! 535: }
! 536:
! 537: /* Y2 -> Y1 */
! 538: def lopitalpf(P,Y1,Y2)
! 539: {
! 540: // if ( !member(Y2,vars(P)) ) return P;
! 541: D = normalizef1(Y2-Y1);
! 542: if ( D == Y2-Y1 ) Sgn = 1;
! 543: else Sgn = -1;
! 544: DY2 = strtov("d"+rtostr(Y2));
! 545: R = [];
! 546: for ( T = reverse(P); T != []; T = cdr(T) ) {
! 547: T0 = car(T);
! 548: C = T0[0]; Op = T0[1];
! 549: for ( U = reverse(C); U != []; U = cdr(U) ) {
! 550: U0 = car(U);
! 551: Nm = U0[0]; Dn = U0[1];
! 552: K = degf(Dn,D);
! 553: if ( !K ) {
! 554: Z = subst_f(Dn,Y2,Y1);
! 555: Sgn1 = Z[0]; Dn1 = Z[1];
! 556: R = addpf([[[[Sgn1*Nm,Dn1]],Op]],R);
! 557: } else {
! 558: Dn1 = removef(Dn,D);
! 559: C1 = [[1,Dn1]];
! 560: Z = [];
! 561: for ( J = 0; J <= K; J++ ) {
! 562: B = fac(K)/fac(J)/fac(K-J);
! 563: C1s = substc(C1,Y2,Y1);
! 564: if ( C1s != [] ) {
! 565: W = [[C1s,dp_ptod(DY2^(K-J),PF_DV)*Op]];
! 566: W = mulcpf([[B,[]]],W);
! 567: Z = addpf(W,Z);
! 568: }
! 569: C1 = diffc(Y2,C1);
! 570: if ( C1 == [] ) break;
! 571: }
! 572: Z = mulcpf([[Sgn^K*Nm/fac(K),[]]],Z);
! 573: R = addpf(Z,R);
! 574: }
! 575: }
! 576: }
! 577: return R;
! 578: }
! 579:
! 580: def lopital_others_pf(P,L) {
! 581: Q = P;
! 582: for ( T = L; T != []; T = cdr(T) ) {
! 583: T0 = car(T);
! 584: I0 = T0[0]; I1 = T0[1]; I = I1-I0+1;
! 585: for ( M = I0; M <= I1; M++ ) {
! 586: Q = lopitalpf(Q,PF_V[I0],PF_V[M]);
! 587: }
! 588: Q = simppf(Q,I0,I);
! 589: Q = adjop(Q,I0,I);
! 590: }
! 591: return Q;
! 592: }
! 593:
! 594: def simpop(Op,I0,NV)
! 595: {
! 596: E = dp_etov(Op);
! 597: R = [];
! 598: for ( J = 0, I = I0; J < NV; I++, J++ ) R = cons(E[I],R);
! 599: R = reverse(qsort(R));
! 600: for ( J = 0, I = I0; J < NV; I++, J++ ) E[I] = R[J];
! 601: return dp_vtoe(E);
! 602: }
! 603:
! 604: def simppf(P,I0,NV)
! 605: {
! 606: R = [];
! 607: for ( T = P; T != []; T = cdr(T) ) {
! 608: T0 = car(T);
! 609: R = addpf([[T0[0],simpop(T0[1],I0,NV)]],R);
! 610: }
! 611: return R;
! 612: }
! 613:
! 614: /* simplify (dy1+...+dyn)^k */
! 615:
! 616: def simp_power1(K,I0,NV)
! 617: {
! 618: P = all_partition(K,NV);
! 619: M = map(coef_partition,P,K,NV,I0);
! 620: for ( R = 0, T = M; T != []; T = cdr(T) )
! 621: R += car(T);
! 622: return R;
! 623: }
! 624:
! 625: def indep_eqs(Eq)
! 626: {
! 627: M = length(Eq);
! 628: D = 0;
! 629: for ( I = 0; I < M; I++ )
! 630: for ( H = Eq[I][0]; H != []; H = cdr(H) ) D += car(H)[1];
! 631: for ( N = 0, T = D; T; T = dp_rest(T), N++ );
! 632: Terms = vector(N);
! 633: for ( I = 0, T = D; T; T = dp_rest(T), I++ ) Terms[I] = dp_ht(T);
! 634: A = matrix(M,N);
! 635: for ( I = 0; I < M; I++ ) {
! 636: for ( H = Eq[I][0]; H != []; H = cdr(H) ) {
! 637: T = car(H)[1];
! 638: for ( K = 0; K < N; K++ )
! 639: if ( Terms[K] == T ) A[I][K] = ctorat(car(H)[0]);
! 640: }
! 641: }
! 642: for ( I = 0; I < M; I++ ) {
! 643: Dn = 1;
! 644: for ( J = 0; J < N; J++ ) if ( A[I][J] ) Dn = ilcm(dn(A[I][J]),Dn);
! 645: for ( J = 0; J < N; J++ ) A[I][J] *= Dn;
! 646: }
! 647: R = indep_rows_mod(A,lprime(0));
! 648: if ( length(R) == N ) {
! 649: L = [];
! 650: for ( I = N-1; I >= 0; I-- )
! 651: L = cons(Eq[R[I]],L);
! 652: return L;
! 653: } else
! 654: return 0;
! 655: }
! 656:
! 657: def eliminatetop(Eq)
! 658: {
! 659: Eq = indep_eqs(Eq);
! 660: if ( !Eq ) return 0;
! 661:
! 662: M = length(Eq);
! 663: R = vector(M);
! 664: D = 0;
! 665: for ( I = 0; I < M; I++ )
! 666: for ( H = Eq[I][0]; H != []; H = cdr(H) ) D += car(H)[1];
! 667: for ( N = 0, T = D; T; T = dp_rest(T), N++ );
! 668: if ( N != M )
! 669: return 0;
! 670: N2 = 2*N;
! 671: Terms = vector(N);
! 672: for ( I = 0, T = D; T; T = dp_rest(T), I++ ) Terms[I] = dp_ht(T);
! 673: A = matrix(N,N2);
! 674: for ( I = 0; I < N; I++ ) {
! 675: for ( H = Eq[I][0]; H != []; H = cdr(H) ) {
! 676: T = car(H)[1];
! 677: for ( K = 0; K < N; K++ )
! 678: if ( Terms[K] == T ) A[I][K] = ctorat(car(H)[0]);
! 679: }
! 680: A[I][N+I] = 1;
! 681: }
! 682: for ( I = 0; I < N; I++ ) {
! 683: Dn = 1;
! 684: for ( J = 0; J < N2; J++ ) if ( A[I][J] ) Dn = ilcm(dn(A[I][J]),Dn);
! 685: for ( J = 0; J < N2; J++ ) A[I][J] *= Dn;
! 686: }
! 687: Ret = generic_gauss_elim(A);
! 688: Num = Ret[0]; Den = Ret[1]; Ind0 = Ret[2]; Ind1 = Ret[3];
! 689: if ( length(Ind0) != N || Ind0[N-1] != N-1 )
! 690: return 0;
! 691: R = [];
! 692: if ( Print ) print(["N=",N]);
! 693: for ( I = 0; I < N; I++ ) {
! 694: if ( Print ) print(["I=",I],2);
! 695: T = [];
! 696: for ( L = 0; L < N; L++ )
! 697: if ( Num[I][L] ) T = addpf(T,mulcpf([[Num[I][L]/Den,[]]],Eq[L][1]));
! 698: R = cons([Terms[I],T],R);
! 699: }
! 700: if ( Print ) print("");
! 701: R = qsort(R,n_wishartd.compt);
! 702: return reverse(R);
! 703: }
! 704:
! 705: def eliminatecheck(Eq,Top)
! 706: {
! 707: Eq = indep_eqs(Eq);
! 708: if ( !Eq ) return 0;
! 709:
! 710: M = length(Eq);
! 711: R = vector(M);
! 712: D = 0;
! 713: for ( I = 0; I < M; I++ )
! 714: for ( H = Eq[I][0]; H != []; H = cdr(H) ) D += car(H)[1];
! 715: for ( N = 0, T = D; T; T = dp_rest(T), N++ );
! 716: if ( N != M )
! 717: return 0;
! 718: N2 = 2*N;
! 719: Terms = vector(N);
! 720: for ( I = 0, T = D; T; T = dp_rest(T), I++ ) Terms[I] = dp_ht(T);
! 721: A = matrix(N,N2);
! 722: for ( I = 0; I < N; I++ ) {
! 723: for ( H = Eq[I][0]; H != []; H = cdr(H) ) {
! 724: T = car(H)[1];
! 725: for ( K = 0; K < N; K++ )
! 726: if ( Terms[K] == T ) A[I][K] = ctorat(car(H)[0]);
! 727: }
! 728: A[I][N+I] = 1;
! 729: }
! 730: for ( I = 0; I < N; I++ ) {
! 731: Dn = 1;
! 732: for ( J = 0; J < N2; J++ ) if ( A[I][J] ) Dn = ilcm(dn(A[I][J]),Dn);
! 733: for ( J = 0; J < N2; J++ ) A[I][J] *= Dn;
! 734: }
! 735: if ( Top ) {
! 736: for ( I = 0; I < N; I++ ) {
! 737: for ( T = J = 0; J < N; J++ ) if ( J != I ) T += abs(A[I][J]);
! 738: if ( abs(A[I][I]) < T )
! 739: if ( Print ) print(["ng",I]);
! 740: }
! 741: } else {
! 742: for ( I = 1; I < N; I++ ) {
! 743: for ( T = J = 0; J < N-2; J++ ) if ( J != I-1 ) T += abs(A[I][J]);
! 744: if ( abs(A[I][I-1]) < T )
! 745: if ( Print ) print(["ng",I]);
! 746: }
! 747: }
! 748:
! 749: Ret = generic_gauss_elim_mod(A,lprime(0));
! 750: Num = Ret[0]; Ind0 = Ret[1]; Ind1 = Ret[2];
! 751: if ( length(Ind0) != N || Ind0[N-1] != N-1 )
! 752: return 0;
! 753: else
! 754: return 1;
! 755: }
! 756:
! 757: def mul_lopitalpf(Z,N,K,I0,I1,E)
! 758: {
! 759: I = I1-I0+1;
! 760: for ( J = I0; J <= N; J++ )
! 761: for ( L = 0; L < E[J]; L++ ) Z = muldpf(PF_V[J],Z);
! 762: for ( J = I0+1; J <= I1; J++ ) {
! 763: Z = lopitalpf(Z,PF_V[I0],PF_V[J]);
! 764: }
! 765: S = simppf(Z,I0,I);
! 766: H = []; L = [];
! 767: for ( ; S != []; S = cdr(S) ) {
! 768: if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
! 769: else L = cons(car(S),L);
! 770: }
! 771: return [reverse(H),reverse(L)];
! 772: }
! 773:
! 774: /* Blocks = [B1,B2,...] */
! 775: /* Bi=[s,e] ys=...=ye f */
! 776:
! 777: def diag0pf(N,Blocks)
! 778: {
! 779: Tlopital = Tred = Tlineq = 0;
! 780: wsetup(N);
! 781: P = vector(N+1);
! 782: for ( J = 1; J <= N; J++ ) P[J] = wishartpf(N,J);
! 783: Simp = [];
! 784: Done = [];
! 785: Len = length(Blocks);
! 786: for ( A = 0; A < Len; A++ ) {
! 787: B = Blocks[A];
! 788: Blocks0 = setminus(Blocks,[B]);
! 789: I0 = B[0];
! 790: I1 = B[1];
! 791: I = I1-I0+1;
! 792: for ( K = I0; K <= I1; K++ )
! 793: P[K] = lopital_others_pf(P[K],Blocks0);
! 794: Rall = [];
! 795: for ( K = I+1; K >= 2; K-- ) {
! 796: if ( Print ) print(["K=",K],2);
! 797: DK = simp_power1(K,I0,I);
! 798: Eq = [];
! 799: TlopitalK = 0;
! 800: for ( T = DK; T; T = dp_rest(T) ) {
! 801: E = dp_etov(T);
! 802: for ( U = I0; U <= I1; U++ )
! 803: if ( E[U] >= 2 )
! 804: break;
! 805: if ( U > I1 ) continue;
! 806: E[U] -= 2;
! 807: Ret = mul_lopitalpf2(P,U,N,K,I0,I1,E);
! 808: Eq = cons(Ret,Eq);
! 809: }
! 810: Eq = reverse(Eq);
! 811:
! 812: for ( Eq0 = [], T = DK; T; T = dp_rest(T) )
! 813: Eq0 = addpf([[[[dp_hc(T),[]]],dp_ht(T)]],Eq0);
! 814:
! 815: DY = mtot(-dy0^K,1);
! 816: if ( K == I+1 ) {
! 817: EqTop = addpf([DY],Eq0);
! 818: } else {
! 819: H = []; L = [];
! 820: for ( S = Eq0 ; S != []; S = cdr(S) ) {
! 821: if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
! 822: else L = cons(car(S),L);
! 823: }
! 824: L = reverse(L); H = reverse(H);
! 825: Eq = cons([H,addpf([DY],L)],Eq);
! 826: }
! 827: T0 = time();
! 828: R = eliminatetop(Eq);
! 829: if ( R )
! 830: Rall = cons(R,Rall);
! 831: else
! 832: return [];
! 833: T1 = time(); Tlineq += T1[0]-T0[0];
! 834: if ( Print ) print(["lineq",T1[0]-T0[0],"sec"]);
! 835: }
! 836: Rall = reverse(Rall);
! 837:
! 838: /* EqTop : dyi0 -> dy -> dyi0 */
! 839: for ( T = Rall; T != []; T = cdr(T) ) {
! 840: if ( Print ) print(["red",length(T)+1],2);
! 841: T0 = time();
! 842: EqTop = reducepf(EqTop,car(T));
! 843: T1 = time(); Tred += T1[0]-T0[0];
! 844: if ( Print ) print([T1[0]-T0[0],"sec"]);
! 845: }
! 846: EqTop = adjop(EqTop,I0,I);
! 847: Done = cons(EqTop,Done);
! 848:
! 849: }
! 850: if ( Print ) print(["lopital",Tlopital,"lineq",Tlineq,"red",Tred]);
! 851: Done = ltov(Done);
! 852: Len = length(Done);
! 853: for ( I = 0; I < Len; I++ ) {
! 854: for ( G = [], J = 0; J < Len; J++ )
! 855: if ( J != I ) G = cons(Done[J],G);
! 856: Done[I] = nfpf(Done[I],G);
! 857: }
! 858: return vtol(Done);
! 859: }
! 860:
! 861: def diagpf(N,Blocks)
! 862: "usage : n_wishartd.diagpf(M,[[S1,E1],[S2,E2],...]),\n where M is the number of variables and [Si,Ei] is a diagonal block and S(i+1)=Ei + 1. For example, n_wishartd.diagpf(10,[[1,9],[10,10]) returns a system of PDEs satisfied by 1F1(y1,...,y1,y10)."
! 863: {
! 864: Tlopital = Tred = Tlineq = 0;
! 865: wsetup(N);
! 866: P = vector(N+1);
! 867: for ( J = 1; J <= N; J++ ) P[J] = wishartpf(N,J);
! 868: Simp = [];
! 869: Done = [];
! 870: Len = length(Blocks);
! 871: for ( A = 0; A < Len; A++ ) {
! 872: B = Blocks[A];
! 873: Blocks0 = setminus(Blocks,[B]);
! 874: I0 = B[0];
! 875: I1 = B[1];
! 876: I = I1-I0+1;
! 877: for ( K = I0; K <= I1; K++ )
! 878: P[K] = lopital_others_pf(P[K],Blocks0);
! 879: Rall = [];
! 880: for ( K = I+1; K >= 2; K-- ) {
! 881: if ( Print ) print(["K=",K],2);
! 882: DK = simp_power1(K,I0,I);
! 883: Eq = [];
! 884: TlopitalK = 0;
! 885: for ( T = DK; T; T = dp_rest(T) ) {
! 886: E = dp_etov(T);
! 887: for ( U = I1; U >= I0; U-- )
! 888: if ( E[U] >= 2 )
! 889: break;
! 890: if ( U < I0 ) continue;
! 891: E[U] -= 2;
! 892: Ret = mul_lopitalpf2(P,U,N,K,I0,I1,E);
! 893: Eq = cons(Ret,Eq);
! 894: }
! 895: Eq = reverse(Eq);
! 896:
! 897: for ( Eq0 = [], T = DK; T; T = dp_rest(T) )
! 898: Eq0 = addpf([[[[dp_hc(T),[]]],dp_ht(T)]],Eq0);
! 899:
! 900: DY = mtot(-dy0^K,1);
! 901: if ( K == I+1 ) {
! 902: EqTop = addpf([DY],Eq0);
! 903: } else {
! 904: H = []; L = [];
! 905: for ( S = Eq0 ; S != []; S = cdr(S) ) {
! 906: if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
! 907: else L = cons(car(S),L);
! 908: }
! 909: L = reverse(L); H = reverse(H);
! 910: Eq = cons([H,addpf([DY],L)],Eq);
! 911: }
! 912: T0 = time();
! 913: R = eliminatetop(Eq);
! 914: if ( R )
! 915: Rall = cons(R,Rall);
! 916: else {
! 917: if ( Print ) print(["lineq retry.."],2);
! 918: Eq1 = [];
! 919: Terms = [];
! 920: for ( T = dp_rest(DK); T; T = dp_rest(T) )
! 921: Terms = cons(dp_ht(T),Terms);
! 922: while ( Terms != [] ) {
! 923: for ( Q = 0; Terms != [] && Q < 10; Q++, Terms = cdr(Terms) ) {
! 924: E = dp_etov(car(Terms));
! 925: if ( E[I0] >= 2 ) {
! 926: E[I0] -= 2;
! 927: Ret = mul_lopitalpf(P[I0],N,K,I0,I1,E);
! 928: Eq1 = cons(Ret,Eq1);
! 929: }
! 930: }
! 931: Eq = append(Eq,Eq1);
! 932: R = eliminatetop(Eq);
! 933: if ( R ) break;
! 934: }
! 935: if ( !R )
! 936: error("diagpf : failed to find a solvable linear system");
! 937: Rall = cons(R,Rall);
! 938: }
! 939: T1 = time(); Tlineq += T1[0]-T0[0];
! 940: if ( Print ) print(["lineq",T1[0]-T0[0],"sec"]);
! 941: }
! 942: Rall = reverse(Rall);
! 943:
! 944: /* EqTop : dyi0 -> dy -> dyi0 */
! 945: for ( T = Rall; T != []; T = cdr(T) ) {
! 946: if ( Print ) print(["red",length(T)+1],2);
! 947: T0 = time();
! 948: EqTop = reducepf(EqTop,car(T));
! 949: T1 = time(); Tred += T1[0]-T0[0];
! 950: if ( Print ) print([T1[0]-T0[0],"sec"]);
! 951: }
! 952: EqTop = adjop(EqTop,I0,I);
! 953: Done = cons(EqTop,Done);
! 954:
! 955: }
! 956: if ( Print ) print(["lopital",Tlopital,"lineq",Tlineq,"red",Tred]);
! 957: Done = ltov(Done);
! 958: Len = length(Done);
! 959: for ( I = 0; I < Len; I++ ) {
! 960: for ( G = [], J = 0; J < Len; J++ )
! 961: if ( J != I ) G = cons(Done[J],G);
! 962: Done[I] = nfpf(Done[I],G);
! 963: }
! 964: return vtol(Done);
! 965: }
! 966:
! 967: def diagcheck(N)
! 968: {
! 969: Tmuld = Tlopital = 0;
! 970: MHist = [];
! 971: wsetup(N);
! 972: P = vector(N+1);
! 973: for ( J = 1; J <= N; J++ ) P[J] = wishartpf(N,J);
! 974: Simp = [];
! 975: Done = [];
! 976: B = [1,N];
! 977: I0 = B[0];
! 978: I1 = B[1];
! 979: I = I1-I0+1;
! 980: for ( K = 2; K <= I+1; K++ ) {
! 981: if ( Print ) print(["K=",K],2);
! 982: DK = simp_power1(K,I0,I);
! 983: Eq = [];
! 984: TlopitalK = 0;
! 985: for ( T = DK; T; T = dp_rest(T) ) {
! 986: E = dp_etov(T);
! 987: for ( U = I0; U <= I1; U++ )
! 988: if ( E[U] >= 2 )
! 989: break;
! 990: if ( U > I1 ) continue;
! 991: E[U] -= 2;
! 992: Ret = mul_lopitalpf2(P,U,N,K,I0,I1,E|high=1);
! 993: Eq = cons(Ret,Eq);
! 994: }
! 995: Eq = reverse(Eq);
! 996:
! 997: for ( Eq0 = [], T = DK; T; T = dp_rest(T) )
! 998: Eq0 = addpf([[[[dp_hc(T),[]]],dp_ht(T)]],Eq0);
! 999:
! 1000: DY = mtot(-dy0^K,1);
! 1001: if ( K == I+1 ) {
! 1002: EqTop = addpf([DY],Eq0);
! 1003: } else {
! 1004: H = []; L = [];
! 1005: for ( S = Eq0 ; S != []; S = cdr(S) ) {
! 1006: if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
! 1007: else L = cons(car(S),L);
! 1008: }
! 1009: L = reverse(L); H = reverse(H);
! 1010: Eq = cons([H,addpf([DY],L)],Eq);
! 1011: }
! 1012: R = eliminatecheck(Eq,K==I+1);
! 1013: if ( !R )
! 1014: return 0;
! 1015: }
! 1016: if ( Print ) print(["muld",Tmuld,"lopital",Tlopital]);
! 1017: return 1;
! 1018: }
! 1019:
! 1020: /* dyi -> dyi/K, dy->dyi */
! 1021: def adjop1(T,I,K)
! 1022: {
! 1023: C = T[0];
! 1024: Op = dp_etov(T[1]);
! 1025: if ( Op[I] ) C = mulc([[1/K,[]]],C);
! 1026: Op[I] += Op[0];
! 1027: Op[0] = 0;
! 1028: return [C,dp_vtoe(Op)];
! 1029: }
! 1030:
! 1031: def adjop(P,I,K)
! 1032: {
! 1033: P = map(adjop1,P,I,K);
! 1034: P = qsort(P,n_wishartd.compt);
! 1035: return reverse(P);
! 1036: }
! 1037:
! 1038: def dnc(C)
! 1039: {
! 1040: D = 1;
! 1041: for ( T = C; T != []; T = cdr(T) ) {
! 1042: T0 = car(T); N = T0[0];
! 1043: M = sdiv(ptozp(N),N);
! 1044: D = ilcm(D,M);
! 1045: }
! 1046: return D;
! 1047: }
! 1048:
! 1049: def pftozpf(F)
! 1050: {
! 1051: D = 1;
! 1052: for ( T = F; T != []; T = cdr(T) ) {
! 1053: T0 = car(T);
! 1054: D = ilcm(D,dnc(T0[0]));
! 1055: }
! 1056: return [D,mulcpf([[D,[]]],F)];
! 1057: }
! 1058:
! 1059: def zpftopf(F)
! 1060: {
! 1061: return mulcpf([[1/F[0],[]]],F[1]);
! 1062: }
! 1063:
! 1064: def addzpf(F1,F2)
! 1065: {
! 1066: D1 = F1[0]; D2 = F2[0];
! 1067: L = ilcm(D1,D2); C1 = L/D1; C2 = L/D2;
! 1068: N = addpf(mulcpf([[L/D1,[]]],F1[1]),mulcpf([[L/D2,[]]],F2[1]));
! 1069: return [L,N];
! 1070: }
! 1071:
! 1072: /* R : reducers */
! 1073: def reducepf(Eq,R)
! 1074: {
! 1075: Ret = pftozpf(Eq); Eq = Ret[1]; Dn = Ret[0];
! 1076: S = [1,[]];
! 1077: for ( T = R, U = []; T != []; T = cdr(T) )
! 1078: U = cons([car(T)[0],pftozpf(car(T)[1])],U);
! 1079: R = reverse(U);
! 1080: for ( T = reverse(Eq); T != []; T = cdr(T) ) {
! 1081: T0 = car(T); C = T0[0]; Op = T0[1];
! 1082: for ( U = (R); U != []; U = cdr(U) ) {
! 1083: U0 = car(U);
! 1084: if ( dp_redble(Op,U0[0]) ) {
! 1085: E = dp_etov(dp_subd(Op,U0[0]));
! 1086: Red = U0[1];
! 1087: N = length(E);
! 1088: for ( J = 1; J < N; J++ )
! 1089: for ( L = 0; L < E[J]; L++ ) Red = [Red[0],muldpf(PF_V[J],Red[1])];
! 1090: Red = [Red[0],mulcpf(C,Red[1])];
! 1091: Red = [Red[0],mulcpf([[-1,[]]],Red[1])];
! 1092: S = addzpf(S,Red);
! 1093: break;
! 1094: }
! 1095: }
! 1096: if ( U == [] ) S = addzpf([1,[T0]],S);
! 1097: }
! 1098: S = [S[0]*Dn,S[1]];
! 1099: return zpftopf(S);
! 1100: }
! 1101:
! 1102: def reduceotherspf(Eq,P,I1,N)
! 1103: {
! 1104: S = [];
! 1105: Z = Eq;
! 1106: while ( Z != [] ) {
! 1107: T0 = car(Z); C = T0[0]; Op = T0[1];
! 1108: for ( I = I1; I <= N; I++ ) {
! 1109: U0 = P[I];
! 1110: M = car(U0)[0][0][0];
! 1111: H = car(U0)[1];
! 1112: if ( dp_redble(Op,H) ) {
! 1113: E = dp_etov(dp_subd(Op,H));
! 1114: Red = U0;
! 1115: Len = length(E);
! 1116: for ( J = 0; J < Len; J++ )
! 1117: for ( L = 0; L < E[J]; L++ ) Red = muldpf(PF_V[J],Red);
! 1118: C1 = mulc([[1/M,[]]],C);
! 1119: Red = mulcpf(C1,Red);
! 1120: Z = subpf(Z,Red);
! 1121: break;
! 1122: }
! 1123: }
! 1124: if ( I > N ) {
! 1125: S = cons(T0,S);
! 1126: Z = cdr(Z);
! 1127: }
! 1128: }
! 1129: return reverse(S);
! 1130: }
! 1131:
! 1132: def print_f(F)
! 1133: {
! 1134: for ( ; F != []; F = cdr(F) ) {
! 1135: F0 = car(F);
! 1136: print("(",0); print(F0[0],0); print(")",0);
! 1137: if ( F0[1] > 1 ) {
! 1138: print("^",0); print(F0[1],0);
! 1139: }
! 1140: if ( cdr(F) != [] ) print("*",0);
! 1141: }
! 1142: }
! 1143:
! 1144: def printc(C)
! 1145: {
! 1146: print("(",0);
! 1147: for ( ; C != []; C = cdr(C) ) {
! 1148: C0 = car(C);
! 1149: print("(",0); print(C0[0],0); print(")",0);
! 1150: if ( C0[1] != [] ) {
! 1151: print("/(",0);
! 1152: print_f(C0[1]);
! 1153: print(")",0);
! 1154: }
! 1155: if ( cdr(C) != [] ) print("+",0);
! 1156: }
! 1157: print(")",0);
! 1158: }
! 1159:
! 1160: def printt(T)
! 1161: {
! 1162: printc(T[0]); print("*",0); print(dp_dtop(T[1],PF_DV),0);
! 1163: }
! 1164:
! 1165: def printpf(F)
! 1166: {
! 1167: for ( ; F != []; F = cdr(F) ) {
! 1168: printt(car(F));
! 1169: if ( cdr(F) != [] ) print("+",0);
! 1170: }
! 1171: print("");
! 1172: }
! 1173:
! 1174: def ftop(F)
! 1175: {
! 1176: P = 1;
! 1177: for ( ; F != []; F = cdr(F) ) {
! 1178: F0 = car(F);
! 1179: P *= F0[0]^F0[1];
! 1180: }
! 1181: return P;
! 1182: }
! 1183:
! 1184: def pftord(F)
! 1185: {
! 1186: R = [];
! 1187: for ( T = F; T != []; T = cdr(T) ) {
! 1188: T0 = car(T); C = T0[0]; Op = T0[1];
! 1189: R = cons([ctord(C),Op],R);
! 1190: }
! 1191: return reverse(R);
! 1192: }
! 1193:
! 1194: def pftop(F)
! 1195: {
! 1196: R = pftord(F);
! 1197: Op = F[0][1];
! 1198: N = length(dp_etov(Op));
! 1199: DY = [];
! 1200: for ( I = N; I >= 0; I-- ) DY = cons(strtov("dy"+rtostr(I)),DY);
! 1201: Lcm = [];
! 1202: for ( T = R; T != []; T = cdr(T) )
! 1203: Lcm = lcmf(Lcm,T[0][0][1]);
! 1204: S = 0;
! 1205: for ( T = R; T != []; T = cdr(T) ) {
! 1206: N = T[0][0][0]*ftop(divf(Lcm,T[0][0][1]));
! 1207: Op = dp_dtop(T[0][1],DY);
! 1208: S += N*Op;
! 1209: }
! 1210: return S;
! 1211: }
! 1212:
! 1213: def ctord(C)
! 1214: {
! 1215: N = 0; D = [];
! 1216: for ( T = reverse(C); T != []; T = cdr(T) ) {
! 1217: T0 = car(T); N0 = T0[0]; D0 = T0[1];
! 1218: L = lcmf(D,D0);
! 1219: /* N/D+N0/D0 = (N*(L/D)+N0*(L/D0))/L */
! 1220: N = N*ftop(divf(L,D))+N0*ftop(divf(L,D0));
! 1221: D = L;
! 1222: }
! 1223: R = [];
! 1224: for ( T = D; T != []; T = cdr(T) ) {
! 1225: T0 = car(T); F = T0[0]; M = T0[1];
! 1226: while ( M > 0 && (Q = tdiv(N,F)) ) {
! 1227: N = Q;
! 1228: M--;
! 1229: }
! 1230: if ( M ) R = cons([F,M],R);
! 1231: }
! 1232: D = reverse(R);
! 1233: return [N,D];
! 1234: }
! 1235:
! 1236:
! 1237: def comppfrd(P,R)
! 1238: {
! 1239: P = qsort(P,n_wishartd.compt); P=reverse(P);
! 1240: R = qsort(R,n_wishartd.compt); R=reverse(R);
! 1241: if ( length(P) != length(R) ) return 0;
! 1242: for ( ; P != []; P = cdr(P), R = cdr(R) ) {
! 1243: P0 = car(P); R0 = car(R);
! 1244: C0 = ctord(P0[0]); C1 = R0[0];
! 1245: D0 = ftop(C0[1]); D1 = ftop(C1[1]);
! 1246: if ( (D0 == D1) && C0[0] == -C1[0] ) continue;
! 1247: if ( (D0 == -D1) && C0[0] == C1[0] ) continue;
! 1248: return 0;
! 1249: }
! 1250: return 1;
! 1251: }
! 1252:
! 1253: def multpf(T,F)
! 1254: {
! 1255: E = dp_etov(T[1]);
! 1256: N = length(E);
! 1257: Z = F;
! 1258: for ( J = 1; J < N; J++ )
! 1259: for ( L = 0; L < E[J]; L++ ) Z = muldpf(PF_V[J],Z);
! 1260: Z = mulcpf(T[0],Z);
! 1261: return Z;
! 1262: }
! 1263:
! 1264: def spolypf(F,G)
! 1265: {
! 1266: F0 = car(F); G0 = car(G);
! 1267: DF = F0[1]; DG = G0[1];
! 1268: L = dp_lcm(DF,DG);
! 1269: F1 = multpf([F0[0],dp_subd(L,DF)],F);
! 1270: G1 = multpf([G0[0],dp_subd(L,DG)],G);
! 1271: S = subpf(F1,G1);
! 1272: return S;
! 1273: }
! 1274:
! 1275: def nfpf(F,G)
! 1276: {
! 1277: Z = F;
! 1278: R = [];
! 1279: while ( Z != [] ) {
! 1280: Z0 = car(Z); C = Z0[0]; D = Z0[1];
! 1281: for ( T = G; T != []; T = cdr(T) ) {
! 1282: T0 = car(T);
! 1283: if ( dp_redble(D,T0[0][1]) ) break;
! 1284: }
! 1285: if ( T != [] ) {
! 1286: CG = ctorat(T0[0][0]);
! 1287: C = mulc(C,[[-1/CG,[]]]);
! 1288: S = multpf([C,dp_subd(D,T0[0][1])],T0);
! 1289: Z = addpf(Z,S);
! 1290: } else {
! 1291: R = cons(Z0,R);
! 1292: Z = cdr(Z);
! 1293: }
! 1294: }
! 1295: S = [];
! 1296: for ( T = R; T != []; T = cdr(T) ) {
! 1297: U = ctord(T[0][0]);
! 1298: if ( U[0] )
! 1299: S = cons(T[0],S);
! 1300: }
! 1301: return S;
! 1302: }
! 1303:
! 1304: def nfpf0(F,G)
! 1305: {
! 1306: Z = F;
! 1307: R = [];
! 1308: while ( Z != [] ) {
! 1309: Z0 = car(Z); C = Z0[0]; D = Z0[1];
! 1310: for ( T = G; T != []; T = cdr(T) ) {
! 1311: T0 = car(T);
! 1312: if ( dp_redble(D,T0[0][1]) ) break;
! 1313: }
! 1314: if ( T != [] ) {
! 1315: CG = ctorat(T0[0][0]);
! 1316: C = mulc(C,[[-1/CG,[]]]);
! 1317: S = multpf([C,dp_subd(D,T0[0][1])],T0);
! 1318: Z = addpf(Z,S);
! 1319: } else {
! 1320: R = cons(Z0,R);
! 1321: Z = cdr(Z);
! 1322: }
! 1323: }
! 1324: S = [];
! 1325: for ( T = R; T != []; T = cdr(T) ) {
! 1326: U = ctord(T[0][0]);
! 1327: if ( U[0] )
! 1328: S = cons(T[0],S);
! 1329: }
! 1330: return S;
! 1331: }
! 1332:
! 1333: def pftoeuler(F)
! 1334: {
! 1335: VDV = [y1,dy1];
! 1336: P = pftop(F);
! 1337: Z = dp_ptod(P,VDV);
! 1338: E = dp_etov(dp_ht(Z));
! 1339: S = E[1]-E[0];
! 1340: if ( S < 0 )
! 1341: error("pftoeuler : invalid input");
! 1342: P *= y1^S;
! 1343: Z = dp_ptod(P,VDV);
! 1344: E = dp_etov(dp_ht(Z));
! 1345: D = E[0];
! 1346: C = vector(D+1);
! 1347: C[0] = 1;
! 1348: for ( I = 1; I <= D; I++ )
! 1349: C[I] = C[I-1]*(t-I+1);
! 1350: R = 0;
! 1351: for ( T = Z; T; T = dp_rest(T) ) {
! 1352: E = dp_etov(dp_ht(T));
! 1353: S = E[0]-E[1];
! 1354: if ( S < 0 )
! 1355: error("pftoeuler : invalid input");
! 1356: R += dp_hc(T)*y1^S*C[E[1]];
! 1357: }
! 1358: return R;
! 1359: }
! 1360:
! 1361: def gammam(M,A)
! 1362: {
! 1363: R = 1;
! 1364: for ( I = 1; I <= M; I++ )
! 1365: R *= mpfr_gamma(A-(I-1)/2);
! 1366: R *= eval(@pi^(1/4*M*(M-1)));
! 1367: return R;
! 1368: }
! 1369:
! 1370: def genprc(M,N,PL,SL,X)
! 1371: {
! 1372: A = (M+1)/2; C = (N+M+1)/2;
! 1373: DetS = 1;
! 1374: TrIS = 0;
! 1375: for ( T = PL, S = SL; T != []; T = cdr(T), S = cdr(S) ) {
! 1376: DetS *= car(S)^car(T);
! 1377: TrIS += car(T)/car(S);
! 1378: }
! 1379: C = 1/eval(DetS^(1/2*N)*2^(1/2*N*M));
! 1380: return C;
! 1381: }
! 1382:
! 1383: /*
! 1384: * PL=[P1,P2,...]: sizes of blocks
! 1385: * SL=[S1,S2,...]: [S1xP1,S2xP2,..]
! 1386: */
! 1387:
! 1388: def prob_by_hgm(M,N,PL,SL,X)
! 1389: "usage : n_wishartd.prob_by_hgm(M,N,[P1,P2,...],[S1,S2,...]|options) where M is the number of variables, N is the degrees of freedom, Pi is the size of i-th block and Si is the value of i-th (repeated) eigenvalue of Sigma. options : rk5=1 => use 5-th order Runge-Kutta method (default : rk4) double=1 => use machine double float (default : MPFR) step=k => set the number of steps=k (default : 10^4) init=x => set the maximal coordinate of the initial point=x (default : 1). For example, n_wishartd.prob_by_hgm(3,10,[2,1],[1/10,1],10|rk5=1,double=1) computes Pr[l1<10|diag(1/10,1/10,1)] by RK5 and machine double float."
! 1390: {
! 1391: A = (M+1)/2; C = (N+M+1)/2;
! 1392:
! 1393: B = []; V = []; Beta = []; SBeta = 0;
! 1394: Prev = 1;
! 1395: for ( T = PL, S = SL; T != []; T = cdr(T), S = cdr(S) ) {
! 1396: V = cons(strtov("y"+rtostr(Prev)),V);
! 1397: B = cons([Prev,Prev+car(T)-1],B);
! 1398: Prev += car(T);
! 1399: Beta = cons(1/(2*car(S)),Beta);
! 1400: SBeta += car(T)/(2*car(S));
! 1401: }
! 1402: if ( Prev != M+1 )
! 1403: error("invalid block specification");
! 1404: B = reverse(B); V = reverse(V);
! 1405: Beta = reverse(Beta);
! 1406:
! 1407: T0 = time();
! 1408: if ( type(Z=getopt(eq)) == -1 )
! 1409: Z = diagpf(M,B);
! 1410: T1 = time(); Tdiag = (T1[0]-T0[0])+(T1[1]-T0[1]);
! 1411:
! 1412: if ( type(Step=getopt(step)) == -1 )
! 1413: Step = 10000;
! 1414:
! 1415: if ( type(Double=getopt(double)) == -1 )
! 1416: Double = 0;
! 1417:
! 1418: Z = subst(Z,a,A,c,C);
! 1419: Init=getopt(init);
! 1420: Rk5 = getopt(rk5);
! 1421: if ( type(Rk5) == -1 ) Rk5 = 0;
! 1422: F = genrk45(Z,M,N,V,PL,SL,Beta,SBeta,X,Rk5,Step|double=Double,init=Init);
! 1423: if ( Print ) print(["diag",Tdiag,"pfaffian",Tp,"ps",Tps,"rk",Trk]);
! 1424: return F;
! 1425: }
! 1426:
! 1427: def prob_by_ps(M,N,PL,SL,X)
! 1428: {
! 1429: A = (M+1)/2; C = (N+M+1)/2;
! 1430:
! 1431: B = []; V = []; Beta = []; SBeta = 0;
! 1432: Prev = 1;
! 1433: for ( T = PL, S = SL; T != []; T = cdr(T), S = cdr(S) ) {
! 1434: V = cons(strtov("y"+rtostr(Prev)),V);
! 1435: B = cons([Prev,Prev+car(T)-1],B);
! 1436: Prev += car(T);
! 1437: Beta = cons(1/(2*car(S)),Beta);
! 1438: SBeta += car(T)/(2*car(S));
! 1439: }
! 1440: if ( Prev != M+1 )
! 1441: error("invalid block specification");
! 1442: B = reverse(B); V = reverse(V);
! 1443: Beta = reverse(Beta);
! 1444:
! 1445: if ( type(Z=getopt(eq)) == -1 )
! 1446: Z = diagpf(M,B);
! 1447:
! 1448: Z = subst(Z,a,A,c,C);
! 1449: GM = gammam(M,(M+1)/2)/gammam(M,(N+M+1)/2);
! 1450: C = GM*genprc(M,N,PL,SL,X);
! 1451:
! 1452: Beta = ltov(Beta)*eval(exp(0));
! 1453: X *= eval(exp(0));
! 1454: L = psvalue(Z,V,Beta*X); PS = L[0];
! 1455: MN2 = M*N/2;
! 1456: ExpF = eval(X^(MN2)*exp(-SBeta*X));
! 1457: return C*L[1]*ExpF;
! 1458: }
! 1459:
! 1460: def ps(Z,V,TD)
! 1461: {
! 1462: Tact = Tgr = Tactmul = 0;
! 1463: G = map(ptozp,map(pftop,Z));
! 1464: M = length(V);
! 1465: N = length(G);
! 1466: for ( I = 0, DV = []; I < M; I++ )
! 1467: DV = cons(strtov("d"+rtostr(V[I])),DV);
! 1468: DV = reverse(DV);
! 1469: VDV = append(V,DV);
! 1470: DG = map(decomp_op,G,V,DV);
! 1471: F = [1];
! 1472: R = 1;
! 1473: Den = 1;
! 1474: for ( I = 1 ; I <= TD; I++ ) {
! 1475: L = create_homo(V,I);
! 1476: FI = L[0]; CV = L[1];
! 1477: Eq = [];
! 1478: for ( K = 0; K < N; K++ ) {
! 1479: P = DG[K]; Len = length(P);
! 1480: /* P(0)F(I)*Den + (P(1)F(I-1)+...+P(I)F(0)) = 0 */
! 1481: T0=time();
! 1482: Lhs = dp_dtop(dp_weyl_act(P[0],dp_ptod(FI,V)),V)*Den;
! 1483: for ( T = F, J = 1; T != [] && J < Len; T = cdr(T), J++ )
! 1484: Lhs += dp_dtop(dp_weyl_act(P[J],dp_ptod(car(T),V)),V);
! 1485: T1=time(); Tact += T1[0]-T0[0];
! 1486: if ( Lhs ) {
! 1487: for ( T = dp_ptod(Lhs,V); T; T = dp_rest(T) ) {
! 1488: Eq = cons(dp_hc(T),Eq);
! 1489: }
! 1490: }
! 1491: }
! 1492: T0=time();
! 1493: #if 0
! 1494: Sol = solve_leq(Eq,CV);
! 1495: #else
! 1496: Sol = nd_f4(Eq,CV,0,0);
! 1497: #endif
! 1498: L = p_true_nf(FI,Sol,CV,0);
! 1499: Den1 = ilcm(Den,L[1]); MI = Den1/L[1];
! 1500: FI = L[0]*(Den1/L[1]);
! 1501: T1=time(); Tgr += T1[0]-T0[0];
! 1502:
! 1503: MJ = Den1/Den; Den = Den1;
! 1504: for ( S = [], T = F; T != []; T = cdr(T) )
! 1505: S = cons(MJ*car(T),S);
! 1506: F = cons(FI,reverse(S));
! 1507: R = R*MJ+car(F);
! 1508: }
! 1509: return R/Den;
! 1510: }
! 1511:
! 1512: def msubst(F,V,X)
! 1513: {
! 1514: M = length(V);
! 1515: for ( J = 0, F0 = F*eval(exp(0)); J < M; J++ )
! 1516: F0 = subst(F0,V[J],X[J]);
! 1517: return F0;
! 1518: }
! 1519:
! 1520: def psvalue(Z,V,X)
! 1521: {
! 1522: Tact = Tgr = Tactmul = 0;
! 1523: G = map(ptozp,map(pftop,Z));
! 1524: M = length(V);
! 1525: N = length(G);
! 1526: for ( I = 0, DV = []; I < M; I++ )
! 1527: DV = cons(strtov("d"+rtostr(V[I])),DV);
! 1528: DV = reverse(DV);
! 1529: VDV = append(V,DV);
! 1530: DG = map(decomp_op,G,V,DV);
! 1531: Prev = getopt(prev);
! 1532: if ( type(Prev) == -1 ) {
! 1533: F = [1];
! 1534: R = 1;
! 1535: Val = eval(exp(0));
! 1536: Den = 1;
! 1537: I = 1;
! 1538: } else {
! 1539: /* Prev = [R/Den,Val1,F,Den] */
! 1540: Den = Prev[3];
! 1541: F = Prev[2];
! 1542: R = Prev[0]*Den;
! 1543: I = length(F);
! 1544: Val = msubst(Prev[0],V,X);
! 1545: ValI = msubst(F[0],V,X)/Den;
! 1546: if ( Val-ValI == Val )
! 1547: return [Prev[0],Val,F,Den];
! 1548: }
! 1549: for ( ; ; I++ ) {
! 1550: L = create_homo(V,I);
! 1551: FI = L[0]; CV = L[1];
! 1552: Eq = [];
! 1553: for ( K = 0; K < N; K++ ) {
! 1554: P = DG[K]; Len = length(P);
! 1555: /* P(0)F(I)*Den + (P(1)F(I-1)+...+P(I)F(0)) = 0 */
! 1556: T0=time();
! 1557: Lhs = dp_dtop(dp_weyl_act(P[0],dp_ptod(FI,V)),V)*Den;
! 1558: for ( T = F, J = 1; T != [] && J < Len; T = cdr(T), J++ )
! 1559: Lhs += dp_dtop(dp_weyl_act(P[J],dp_ptod(car(T),V)),V);
! 1560: T1=time(); Tact += T1[0]-T0[0];
! 1561: if ( Lhs ) {
! 1562: for ( T = dp_ptod(Lhs,V); T; T = dp_rest(T) ) {
! 1563: Eq = cons(dp_hc(T),Eq);
! 1564: }
! 1565: }
! 1566: }
! 1567: T0=time();
! 1568: Sol = solve_leq(Eq,CV);
! 1569: L = p_true_nf(FI,Sol,CV,0);
! 1570: Den1 = ilcm(Den,L[1]); MI = Den1/L[1];
! 1571: FI = L[0]*(Den1/L[1]);
! 1572: T1=time(); Tgr += T1[0]-T0[0];
! 1573:
! 1574: MJ = Den1/Den; Den = Den1;
! 1575: for ( S = [], T = F; T != []; T = cdr(T) )
! 1576: S = cons(MJ*car(T),S);
! 1577: F = cons(FI,reverse(S));
! 1578: R = R*MJ+car(F);
! 1579: F0 = msubst(FI,V,X)/Den;
! 1580: Val1 = Val+F0;
! 1581: if ( Val1 == Val ) {
! 1582: if ( Print ) print(["I=",I,"act",Tact,"actmul",Tactmul,"gr",Tgr]);
! 1583: delete_uc();
! 1584: return [R/Den,Val1,F,Den];
! 1585: } else {
! 1586: if ( Print ) print([F0],2);
! 1587: Val = Val1;
! 1588: }
! 1589: }
! 1590: }
! 1591:
! 1592: /* P -> [P0,P1,...] Pi = y^a dy^b, |a|-|b|=i */
! 1593:
! 1594: def decomp_op(P,V,DV)
! 1595: {
! 1596: M = length(V);
! 1597: VDV = append(V,DV);
! 1598: D = dp_ptod(P,VDV);
! 1599: for ( I = 0, T = D; T; T = dp_rest(T), I++ );
! 1600: for ( T = D, NotSet = 1; T; T = dp_rest(T) ) {
! 1601: E = dp_etov(T);
! 1602: for ( K = 0, J = 0; J < M; J++ )
! 1603: K += E[J]-E[M+J];
! 1604: if ( NotSet ) {
! 1605: Max = Min = K; NotSet = 0;
! 1606: } else if ( K > Max ) Max = K;
! 1607: else if ( K < Min ) Min = K;
! 1608: }
! 1609: W = vector(Max-Min+1);
! 1610: for ( T = D; T; T = dp_rest(T) ) {
! 1611: E = dp_etov(T);
! 1612: for ( K = 0, J = 0; J < M; J++ )
! 1613: K += E[J]-E[M+J];
! 1614: W[K-Min] += dp_hm(T);
! 1615: }
! 1616: W = map(dp_ptod,map(dp_dtop,W,VDV),DV);
! 1617: return W;
! 1618: }
! 1619:
! 1620: def create_homo(V,TD)
! 1621: {
! 1622: M = length(V);
! 1623: for ( S = 0, I = 0; I < M; I++ ) S += V[I];
! 1624: Tmp = 0;
! 1625: Nc = 0;
! 1626: for ( RI = 0, D = dp_ptod(S^TD,V); D; D = dp_rest(D), Nc++ ) {
! 1627: E = dp_etov(D);
! 1628: T = uc();
! 1629: U = dp_dtop(dp_ht(D),V);
! 1630: RI += T*U;
! 1631: Tmp += T;
! 1632: }
! 1633: CV = vector(Nc);
! 1634: for ( I = 0; I < Nc; I++ ) {
! 1635: CV[I] = var(Tmp); Tmp -= CV[I];
! 1636: }
! 1637: return [RI,vtol(CV)];
! 1638: }
! 1639:
! 1640: def act_op(P,V,F)
! 1641: {
! 1642: N = length(V);
! 1643: for ( R = 0; P; P = dp_rest(P) ) {
! 1644: C = dp_hc(P); E = dp_etov(P);
! 1645: for ( I = 0, S = F; I < N; I++ )
! 1646: for ( J = 0; J < E[I]; J++ ) S = diff(S,V[I]);
! 1647: T0 = time();
! 1648: R += C*S;
! 1649: T1 = time(); Tactmul += T1[0]-T0[0];
! 1650: }
! 1651: return R;
! 1652: }
! 1653:
! 1654: def gen_basis(D,DV)
! 1655: {
! 1656: N = length(D);
! 1657: B = [];
! 1658: E = vector(N);
! 1659: while ( 1 ) {
! 1660: B = cons(dp_dtop(dp_vtoe(E),DV),B);
! 1661: for ( I = N-1; I >= 0; I-- )
! 1662: if ( E[I] < D[I]-1 ) break;
! 1663: if ( I < 0 ) return reverse(B);
! 1664: E[I]++;
! 1665: for ( J = I+1; J < N; J++ ) E[J] = 0;
! 1666: }
! 1667: }
! 1668:
! 1669: def pfaffian(Z)
! 1670: {
! 1671: N = length(Z);
! 1672: D = vector(N);
! 1673: Mat = vector(N);
! 1674: V = vars(Z);
! 1675: DV = vector(N);
! 1676: for ( I = 0; I < N; I++ )
! 1677: DV[I] = strtov("d"+rtostr(V[I]));
! 1678: DV = vtol(DV);
! 1679: for ( I = 0; I < N; I++ ) {
! 1680: ZI = Z[I];
! 1681: HI = ZI[0][1];
! 1682: DI = dp_dtop(HI,PF_DV);
! 1683: for ( J = 0; J < N; J++ )
! 1684: if ( var(DI) == DV[J] )
! 1685: break;
! 1686: D[J] = deg(DI,DV[J]);
! 1687: }
! 1688: B = gen_basis(D,DV);
! 1689: Dim = length(B);
! 1690: Hist = [];
! 1691: for ( I = 0; I < N; I++ ) {
! 1692: if ( Print ) print(["I=",I]);
! 1693: A = matrix(Dim,Dim);
! 1694: for ( K = 0; K < Dim; K++ )
! 1695: for ( L = 0; L < Dim; L++ )
! 1696: A[K][L] = [];
! 1697: for ( J = 0; J < Dim; J++ ) {
! 1698: if ( Print ) print(["J=",J],2);
! 1699: Mon = DV[I]*B[J];
! 1700: for ( T = Hist; T != []; T = cdr(T) )
! 1701: if ( (Q = tdiv(Mon,car(T)[0])) && Q==var(Q) ) break;
! 1702: if ( (T != []) ) {
! 1703: for ( L = 0; L < N; L++ )
! 1704: if ( Q == DV[L] ) break;
! 1705: Red1 = muldpf(V[L],car(T)[1]);
! 1706: Red = nfpf0(Red1,Z);
! 1707: } else
! 1708: Red = nfpf0([mtot(Mon,1)],Z);
! 1709: Hist = cons([Mon,Red],Hist);
! 1710: for ( T = Red; T != []; T = cdr(T) ) {
! 1711: T0 = car(T); C = T0[0]; Op = dp_dtop(T0[1],PF_DV);
! 1712: for ( K = 0; K < Dim; K++ )
! 1713: if ( B[K] == Op ) break;
! 1714: if ( K == Dim )
! 1715: error("afo");
! 1716: else
! 1717: A[J][K] = C;
! 1718: }
! 1719: }
! 1720: if ( Print ) print("");
! 1721: A = map(ctord,A);
! 1722: Lcm = [];
! 1723: for ( K = 0; K < Dim; K++ )
! 1724: for ( L = 0; L < Dim; L++ )
! 1725: Lcm = lcmf(Lcm,A[K][L][1]);
! 1726: for ( K = 0; K < Dim; K++ )
! 1727: for ( L = 0; L < Dim; L++ ) {
! 1728: Num = A[K][L][0]; Den = A[K][L][1];
! 1729: A[K][L] = Num*ftorat(divf(Lcm,Den));
! 1730: }
! 1731: Lcm = ftorat(Lcm);
! 1732: if ( !Lcm ) Lcm = 1;
! 1733: Mat[I] = [A,Lcm];
! 1734: }
! 1735: return [Mat,B];
! 1736: }
! 1737:
! 1738: def pfaffian2(Z,V,Beta,XV)
! 1739: {
! 1740: N = length(Z);
! 1741: D = vector(N);
! 1742: Mat = vector(N);
! 1743: DV = vector(N);
! 1744: for ( I = 0; I < N; I++ )
! 1745: DV[I] = strtov("d"+rtostr(V[I]));
! 1746: DV = vtol(DV);
! 1747: for ( I = 0; I < N; I++ ) {
! 1748: ZI = Z[I];
! 1749: HI = ZI[0][1];
! 1750: DI = dp_dtop(HI,PF_DV);
! 1751: for ( J = 0; J < N; J++ )
! 1752: if ( var(DI) == DV[J] )
! 1753: break;
! 1754: D[J] = deg(DI,DV[J]);
! 1755: }
! 1756: B = gen_basis(D,DV);
! 1757: Dim = length(B);
! 1758: Hist = [];
! 1759: for ( I = 0; I < N; I++ ) {
! 1760: if ( Print ) print(["I=",I]);
! 1761: A = matrix(Dim,Dim);
! 1762: for ( K = 0; K < Dim; K++ )
! 1763: for ( L = 0; L < Dim; L++ )
! 1764: A[K][L] = [];
! 1765: for ( J = 0; J < Dim; J++ ) {
! 1766: if ( Print ) print(["J=",J],2);
! 1767: Mon = DV[I]*B[J];
! 1768: for ( T = Hist; T != []; T = cdr(T) )
! 1769: if ( (Q = tdiv(Mon,car(T)[0])) && Q==var(Q) ) break;
! 1770: if ( (T != []) ) {
! 1771: for ( L = 0; L < N; L++ )
! 1772: if ( Q == DV[L] ) break;
! 1773: Red1 = muldpf(V[L],car(T)[1]);
! 1774: Red = nfpf0(Red1,Z);
! 1775: } else
! 1776: Red = nfpf0([mtot(Mon,1)],Z);
! 1777: Hist = cons([Mon,Red],Hist);
! 1778: for ( T = Red; T != []; T = cdr(T) ) {
! 1779: T0 = car(T); C = T0[0]; Op = dp_dtop(T0[1],PF_DV);
! 1780: for ( K = 0; K < Dim; K++ )
! 1781: if ( B[K] == Op ) break;
! 1782: if ( K == Dim )
! 1783: error("afo");
! 1784: else
! 1785: A[J][K] = C;
! 1786: }
! 1787: }
! 1788: if ( Print ) print("");
! 1789: A = map(ctord,A);
! 1790: Lcm = [];
! 1791: for ( K = 0; K < Dim; K++ )
! 1792: for ( L = 0; L < Dim; L++ )
! 1793: Lcm = lcmf(Lcm,A[K][L][1]);
! 1794: for ( K = 0; K < Dim; K++ )
! 1795: for ( L = 0; L < Dim; L++ ) {
! 1796: Num = A[K][L][0]; Den = A[K][L][1];
! 1797: A[K][L] = Num*ftorat(divf(Lcm,Den));
! 1798: }
! 1799: Lcm = ftorat(Lcm);
! 1800: if ( !Lcm ) Lcm = 1;
! 1801: for ( K = 0; K < N; K++ ) {
! 1802: A = subst(A,V[K],Beta[K]*XV);
! 1803: Lcm = subst(Lcm,V[K],Beta[K]*XV);
! 1804: }
! 1805: A = map(red,A/Lcm);
! 1806: Lcm = 1;
! 1807: for ( K = 0; K < Dim; K++ )
! 1808: for ( L = 0; L < Dim; L++ )
! 1809: Lcm = lcm(dn(A[K][L]),Lcm);
! 1810: for ( K = 0; K < Dim; K++ )
! 1811: for ( L = 0; L < Dim; L++ )
! 1812: A[K][L] = nm(A[K][L])*sdiv(Lcm,dn(A[K][L]));
! 1813: Mat[I] = [A,Lcm];
! 1814: }
! 1815: Lcm = 1;
! 1816: for ( I = 0; I < N; I++ )
! 1817: Lcm = lcm(Mat[I][1],Lcm);
! 1818: R = matrix(Dim,Dim);
! 1819: for ( I = 0; I < N; I++ ) {
! 1820: A = Mat[I][0];
! 1821: for ( K = 0; K < Dim; K++ )
! 1822: for ( L = 0; L < Dim; L++ )
! 1823: R[K][L] += Beta[I]*nm(A[K][L])*sdiv(Lcm,Mat[I][1]);
! 1824: }
! 1825: Deg = 0;
! 1826: for ( K = 0; K < Dim; K++ )
! 1827: for ( L = 0; L < Dim; L++ ) {
! 1828: DegT = deg(R[K][L],XV);
! 1829: if ( DegT > Deg ) Deg = DegT;
! 1830: }
! 1831: RA = vector(Deg+1);
! 1832: for ( I = 0; I <= Deg; I++ )
! 1833: RA[I] = map(coef,R,I);
! 1834: return [[RA,Lcm],B];
! 1835: }
! 1836:
! 1837: def evalpf(F,V,X)
! 1838: {
! 1839: if ( !F ) return 0;
! 1840: F0 = F;
! 1841: N = length(V);
! 1842: for ( I = 0; I < N; I++ )
! 1843: F0 = subst(F0,V[I],X[I]);
! 1844: return F0;
! 1845: }
! 1846:
! 1847: def eval_pfaffian(Mat,Beta,SBeta,MN2,V,XI,F)
! 1848: {
! 1849: N = length(V);
! 1850: Dim = size(Mat[0][0])[0];
! 1851: R = vector(Dim);
! 1852: Mul = vector(N);
! 1853: X = XI*Beta;
! 1854: X = vtol(X);
! 1855: for ( K = 0; K < N; K++ )
! 1856: Mul[K] = Beta[K]/evalpf(Mat[K][1],V,X);
! 1857: for ( K = 0; K < N; K++ ) {
! 1858: MatK = Mat[K][0];
! 1859: for ( I = 0; I < Dim; I++ ) {
! 1860: for ( U = J = 0; J < Dim; J++ )
! 1861: U += substr2np(MatK[I][J],V,X)*F[J];
! 1862: R[I] += Mul[K]*U;
! 1863: }
! 1864: }
! 1865: for ( I = 0; I < Dim; I++ ) R[I] -= (SBeta-MN2/XI)*F[I];
! 1866: return R;
! 1867: }
! 1868:
! 1869: def eval_pfaffian2(Mat,Beta,SBeta,MN2,V,T,XI,F)
! 1870: {
! 1871: N = length(V);
! 1872: MA = Mat[0]; Den = subst(Mat[1],T,XI);
! 1873: Len = length(MA);
! 1874: Dim = size(MA[0])[0];
! 1875: R = vector(Dim);
! 1876: XII = 1;
! 1877: for ( I = 0; I < Len; I++ ) {
! 1878: R += MA[I]*(XII*F);
! 1879: XII *= XI;
! 1880: }
! 1881: R = R/Den - (SBeta-MN2/XI)*F;
! 1882: return R;
! 1883: }
! 1884:
! 1885:
! 1886: def genrk45(Z,M,N,V,PL,SL,Beta,SBeta,X,RK5,Step)
! 1887: {
! 1888: if ( type(Double=getopt(double)) == -1 ) {
! 1889: ctrl("bigfloat",1);
! 1890: Double = 0;
! 1891: } else
! 1892: ctrl("bigfloat",0);
! 1893: if ( type(Init=getopt(init)) == -1 )
! 1894: Init = 1;
! 1895: Len = length(V);
! 1896: DV = vector(Len);
! 1897: for ( I = 0; I < Len; I++ )
! 1898: DV[I] = strtov("d"+rtostr(V[I]));
! 1899: DV = vtol(DV);
! 1900: A = (1+M)/2; C = (1+M+N)/2;
! 1901: Z0 = subst(Z,a,A,c,C);
! 1902: T0 = time();
! 1903: Q = pfaffian2(Z0,V,Beta,t); Mat = Q[0]; Base = Q[1];
! 1904: MatLen = length(Q[0][0]);
! 1905: for ( I = 0; I < MatLen; I++ )
! 1906: Mat[0][I] *= eval(exp(0));
! 1907: T1 = time(); Tp = (T1[0]-T0[0])+(T1[1]-T0[1]);
! 1908: T0 = time();
! 1909: Beta = ltov(Beta)*eval(exp(0));
! 1910: X *= eval(exp(0));
! 1911: X1 = Beta*X;
! 1912: X0 = 0;
! 1913: for ( I = 0; I < Len; I++ )
! 1914: if ( Beta[I] > X0 ) X0 = Beta[I];
! 1915: X0 /= Init;
! 1916: /* Max Beta[I] is set to Init */
! 1917: Beta0 = Beta/X0;
! 1918: X0 = 1/X0;
! 1919: /* now Beta0 = X0*Beta */
! 1920: #if 0
! 1921: Prec = setbprec();
! 1922: setbprec(2*Prec);
! 1923: #endif
! 1924: L = psvalue(Z0,V,Beta0); PS = L[0];
! 1925: T1 = time(); Tps = (T1[0]-T0[0])+(T1[1]-T0[1]);
! 1926: T0 = time();
! 1927: #if 0
! 1928: setbprec(Prec);
! 1929: #endif
! 1930: Dim = length(Base);
! 1931: MN2 = M*N/2;
! 1932: ExpF = eval(X0^(MN2)*exp(-SBeta*X0));
! 1933: F = vector(Dim);
! 1934: for ( I = 0; I < Dim; I++ ) {
! 1935: DPS = act_op(dp_ptod(Base[I],DV),V,PS);
! 1936: for ( J = 0; J < Len; J++ )
! 1937: DPS = subst(DPS,V[J],Beta0[J]);
! 1938: F[I] = DPS*ExpF;
! 1939: }
! 1940: F0 = F*1;
! 1941: while ( 1 ) {
! 1942: F = F0*1;
! 1943: H = eval((X-X0)/Step);
! 1944: if ( Double ) {
! 1945: ctrl("bigfloat",0);
! 1946: X0 = todouble(X0);
! 1947: H = todouble(H);
! 1948: Beta = map(todouble,Beta);
! 1949: SBeta = todouble(SBeta);
! 1950: F = map(todouble,F);
! 1951: }
! 1952: R = [];
! 1953: GM = gammam(M,(M+1)/2)/gammam(M,(N+M+1)/2);
! 1954: O = eval(exp(0));
! 1955: if ( RK5 ) {
! 1956: A2 = 1/5*O; B21 = 1/5*O;
! 1957: A3 = 3/10*O;B31 = 3/40*O; B32 = 9/40*O;
! 1958: A4 = 3/5*O; B41 = 3/10*O; B42 = -9/10*O; B43 = 6/5*O;
! 1959: A5 = 1*O; B51 = -11/54*O; B52 = 5/2*O; B53 = -70/27*O; B54 = 35/27*O;
! 1960: A6 = 7/8*O; B61 = 1631/55296*O; B62 = 175/512*O; B63 = 575/13824*O; B64 = 44275/110592*O; B65 = 253/4096*O;
! 1961: C1 = 37/378*O; C2 = 0*O; C3 = 250/621*O; C4 = 125/594*O; C5 = 0*O; C6 = 512/1771*O;
! 1962: D1 = 2825/27648*O; D2 = 0*O; D3 = 18575/48384*O; D4 = 13525/55296*O; D5 = 277/14336*O; D6 = 1/4*O;
! 1963: for ( I = 0; I < Step; I++ ) {
! 1964: if ( Print ) {
! 1965: if ( !(I%1000) ) print([I],2);
! 1966: if ( !(I%10000) ) print("");
! 1967: }
! 1968: XI = X0+I*H;
! 1969: K1 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI,F);
! 1970: K2 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+A2*H,F+B21*K1);
! 1971: K3 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+A3*H,F+B31*K1+B32*K2);
! 1972: K4 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+A4*H,F+B41*K1+B42*K2+B43*K3);
! 1973: K5 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+A5*H,F+B51*K1+B52*K2+B53*K3+B54*K4);
! 1974: K6 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+A6*H,F+B61*K1+B62*K2+B63*K3+B64*K4+B65*K5);
! 1975: F += C1*K1+C2*K2+C3*K3+C4*K4+C5*K5+C6*K6;
! 1976: T = GM*F[0]*genprc(M,N,PL,SL,XI+H);
! 1977: if ( T < 0 || T > 1 ) break;
! 1978: R = cons([XI+H,T],R);
! 1979: }
! 1980: } else {
! 1981: for ( I = 0; I < Step; I++ ) {
! 1982: if ( Print ) {
! 1983: if ( !(I%1000) ) print([I],2);
! 1984: if ( !(I%10000) ) print("");
! 1985: }
! 1986: XI = X0+I*H;
! 1987: K1 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI,F);
! 1988: K2 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+H/2,F+1/2*K1);
! 1989: K3 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+H/2,F+1/2*K2);
! 1990: K4 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+H,F+K3);
! 1991: F += 1/6*K1+1/3*K2+1/3*K3+1/6*K4;
! 1992: T = GM*F[0]*genprc(M,N,PL,SL,XI+H);
! 1993: if ( T < 0 || T > 1 ) break;
! 1994: R = cons([XI+H,T],R);
! 1995: }
! 1996: }
! 1997: if ( Print ) print("");
! 1998: if ( I == Step ) {
! 1999: T1 = time(); Trk = (T1[0]-T0[0])+(T1[1]-T0[1]);
! 2000: return reverse(R);
! 2001: } else {
! 2002: Step *= 2;
! 2003: if ( Print ) print("Step=",0); if ( Print ) print(Step);
! 2004: }
! 2005: }
! 2006: }
! 2007:
! 2008: def solve_leq(L,V)
! 2009: {
! 2010: N = length(V);
! 2011: B = [];
! 2012: for ( I = 0; I < N; I++, L = cdr(L) ) B = cons(car(L),B);
! 2013: while ( 1 ) {
! 2014: Sol = nd_f4(B,V,0,0);
! 2015: if ( length(Sol) != N ) {
! 2016: B = cons(car(L),Sol); L = cdr(L);
! 2017: } else {
! 2018: return Sol;
! 2019: }
! 2020: }
! 2021: }
! 2022:
! 2023:
! 2024: def ptom(P,V)
! 2025: {
! 2026: M = length(P); N = length(V);
! 2027: A = matrix(M,N+1);
! 2028: for ( I = 0; I < M; I++ ) {
! 2029: F = ptozp(P[I]);
! 2030: VT = V;
! 2031: J = 0;
! 2032: while ( F )
! 2033: if ( type(F) <= 1 ) {
! 2034: A[I][N] = F; break;
! 2035: } else {
! 2036: for ( FV = var(F); FV != car(VT); VT = cdr(VT), J++ );
! 2037: A[I][J] = coef(F,1);
! 2038: F -= A[I][J]*FV;
! 2039: }
! 2040: }
! 2041: return A;
! 2042: }
! 2043:
! 2044: def vton(V1,V)
! 2045: {
! 2046: N = length(V);
! 2047: for ( I = 0; I < N; I++ )
! 2048: if ( V1 == V[I] ) return I;
! 2049: error("vton : no such variable");
! 2050: }
! 2051:
! 2052: def ftotex(F,V)
! 2053: {
! 2054: if ( F == [] ) return "1";
! 2055: R = "";
! 2056: for ( T = F; T != []; T = cdr(T) ) {
! 2057: Fac = car(T)[0]; M = car(T)[1];
! 2058: V1 = var(Fac); NV1 = vton(V1,V);
! 2059: if ( Fac == V1 ) SFac = "y_{"+rtostr(NV1)+"}";
! 2060: else {
! 2061: V2 = var(Fac-V1);
! 2062: NV2 = vton(V2,V);
! 2063: SFac = "(y_{"+rtostr(NV1)+"}-y_{"+rtostr(NV2)+"})";
! 2064: }
! 2065: if ( M == 1 ) R += SFac;
! 2066: else R += SFac+"^{"+rtostr(M)+"}";
! 2067: }
! 2068: return R;
! 2069: }
! 2070:
! 2071: def ctotex(C,V)
! 2072: {
! 2073: R = "";
! 2074: for ( T = C; T != []; T = cdr(T) ) {
! 2075: if ( T != C ) R += "+";
! 2076: N = quotetotex(objtoquote(car(T)[0]));
! 2077: if ( car(T)[1] == [] ) {
! 2078: R += "("+N+")";
! 2079: } else {
! 2080: D = ftotex(car(T)[1],V);
! 2081: R += "\\frac{"+N+"}{"+D+"}";
! 2082: }
! 2083: }
! 2084: return R;
! 2085: }
! 2086:
! 2087: def optotex(Op,DV)
! 2088: {
! 2089: E = dp_etov(Op);
! 2090: N = length(E);
! 2091: R = "";
! 2092: for ( I = 0; I < N; I++ )
! 2093: if ( E[I] )
! 2094: if ( E[I] == 1 )
! 2095: R += "\\partial_{"+rtostr(I)+"}";
! 2096: else
! 2097: R += "\\partial_{"+rtostr(I)+"}^{"+rtostr(E[I])+"}";
! 2098: return R;
! 2099: }
! 2100:
! 2101: def pftotex(P,V,DV)
! 2102: {
! 2103: R = "";
! 2104: for ( T = P; T != []; T = cdr(T) ) {
! 2105: if ( T != P ) R += "+";
! 2106: C = ctotex(car(T)[0],V); Op = optotex(car(T)[1],DV);
! 2107: R += "("+C+")"+Op+"\\\\";
! 2108: }
! 2109: return R;
! 2110: }
! 2111:
! 2112: def eqtotex(Z)
! 2113: {
! 2114: shell("rm -f __tmp__.tex");
! 2115: output("__tmp__.tex");
! 2116: N = length(Z);
! 2117: print("\\documentclass[12pt]{jsarticle}");
! 2118: print("\\begin{document}");
! 2119: print("\\large\\noindent");
! 2120: for ( I = 0; I < N; I++ ) {
! 2121: print("$h_"+rtostr(I)+"=",0);
! 2122: T = mulcpf([[-1,[]]],Z[I]);
! 2123: print(pftotex(T,PF_V,PF_DV),0);
! 2124: print("$\\\\");
! 2125: }
! 2126: print("\\end{document}");
! 2127: output();
! 2128: shell("latex __tmp__");
! 2129: shell("xdvi __tmp__");
! 2130: }
! 2131:
! 2132: /*
! 2133: * for a partition L=[N1,...,Nl]
! 2134: * compute the coefficient of x1^N1*...*xl^Nl in phi((x1+...+xM)^N)
! 2135: * where phi(x(s(1))^n1*...*x(s(k))^nk)=x1^n1*...*xk^nk (n1>=...>=nk)
! 2136: * if count_dup(L)=[I1,...,Is], then the coefficient is
! 2137: * C=N!/(N1!*...*Nl!)*M!/((M-l)!*I1!*...*Is!)
! 2138: * return C*<<0,...0,N1,...,Nl,0,...,0>>
! 2139: position 0 Base
! 2140: */
! 2141:
! 2142: def coef_partition(L,N,M,Base)
! 2143: {
! 2144: R = fac(N);
! 2145: for ( T = L; T != []; T = cdr(T) )
! 2146: R = sdiv(R,fac(car(T)));
! 2147: S = count_dup(L);
! 2148: R *= fac(M)/fac(M-length(L));
! 2149: for ( T = S; T != []; T = cdr(T) )
! 2150: R = sdiv(R,fac(car(T)));
! 2151: E = vector(length(PF_DV));
! 2152: for ( I = Base, T = L; T != []; T = cdr(T), I++ )
! 2153: E[I] = car(T);
! 2154: return R*dp_vtoe(E);
! 2155: }
! 2156:
! 2157: /*
! 2158: * for a partition L=[N1,...,Nk], return [I1,...,Is]
! 2159: * s.t. N1=...=N(I1)>N(I1+1)=...=N(I1+I2)>...
! 2160: */
! 2161:
! 2162: def count_dup(L)
! 2163: {
! 2164: D = [];
! 2165: T = L;
! 2166: while ( T != [] ) {
! 2167: T0 = car(T); T = cdr(T);
! 2168: for ( I = 1; T != [] && car(T) == T0; T = cdr(T) ) I++;
! 2169: D = cons(I,D);
! 2170: }
! 2171: return D;
! 2172: }
! 2173:
! 2174: /* generate (N_1,...,N_L) s.t. N_1+...+N_L=N, N_1>=...>=N_L>=1, L<= K */
! 2175:
! 2176: def all_partition(N,K)
! 2177: {
! 2178: R = [];
! 2179: for ( I = 1; I <= K; I++ )
! 2180: R = append(partition(N,I),R);
! 2181: return map(reverse,R);
! 2182: }
! 2183:
! 2184: /* generate (N_1,...,N_K) s.t. N_1+...+N_K=N, 1<=N_1<=...<=N_K */
! 2185:
! 2186: def partition(N,K)
! 2187: {
! 2188: if ( N < K ) return [];
! 2189: else if ( N == K ) {
! 2190: S = [];
! 2191: for ( I = 0; I < K; I++ )
! 2192: S = cons(1,S);
! 2193: return [S];
! 2194: } else if ( K == 0 )
! 2195: return [];
! 2196: else {
! 2197: R = [];
! 2198: L = partition(N-K,K);
! 2199: for ( T = L; T != []; T = cdr(T) ) {
! 2200: S = [];
! 2201: for ( U = car(T); U != []; U = cdr(U) )
! 2202: S = cons(car(U)+1,S);
! 2203: R = cons(reverse(S),R);
! 2204: }
! 2205: L = partition(N-1,K-1);
! 2206: for ( T = L; T != []; T = cdr(T) )
! 2207: R = cons(cons(1,car(T)),R);
! 2208: return R;
! 2209: }
! 2210: }
! 2211:
! 2212: def mul_lopitalpf2(P,I,N,K,I0,I1,E)
! 2213: {
! 2214: YI = PF_V[I]; DYI = PF_DV[I];
! 2215: R = [mtot(DYI^2,1)];
! 2216: for ( J = I0; J <= I1; J++ ) {
! 2217: if ( J == I ) continue;
! 2218: YJ = PF_V[J]; DYJ = PF_DV[J];
! 2219: R = addpf([mtot(1/2*DYI,YI-YJ)],R);
! 2220: R = addpf([mtot(-1/2*DYJ,YI-YJ)],R);
! 2221: }
! 2222: S = subpf(P[I],R);
! 2223: R = loph(E,I,I0,I1);
! 2224: if ( type(getopt(high)) == -1 )
! 2225: S = mul_lopitalpf(S,N,K,I0,I1,E)[1];
! 2226: else
! 2227: S = 0;
! 2228: return [R,S];
! 2229: }
! 2230:
! 2231: /* the highest part of lim_{Y(I+1),...,Y(I+N-1)->YI} E(PI) */
! 2232: /* where E1 = E+(0,...,2,...,0) */
! 2233:
! 2234: def loph(E,I,I0,I1)
! 2235: {
! 2236: E1 = E*1;
! 2237: E1[I] += 2;
! 2238: R = [[[[1,[]]],dp_vtoe(E1)]];
! 2239: S = lopn(E,I,I0,I1);
! 2240: S = mulcpf([[1/2,[]]],S);
! 2241: R = addpf(R,S);
! 2242: R = simppf(R,I0,I1-I0+1);
! 2243: return R;
! 2244: }
! 2245:
! 2246: /* lim_{Y(I+1),...,Y(I+N-1)->YI} dy^E(1/(YI-Y(I+1))+...+1/(YI-Y(I+N-1))) */
! 2247:
! 2248: def lopn(E,I,I0,I1)
! 2249: {
! 2250: R = [];
! 2251: for ( K = I0; K <= I1; K++ ) {
! 2252: if ( K != I ) {
! 2253: L = lop1(E,I,K);
! 2254: R = addpf(R,L);
! 2255: }
! 2256: }
! 2257: return R;
! 2258: }
! 2259:
! 2260: /* lim_{YJ->YI} dy^E(1/(YI-YJ) */
! 2261: def lop1(E,I,J)
! 2262: {
! 2263: YI = PF_V[I]; YJ = PF_V[J];
! 2264: DYI = PF_DV[I]; DYJ = PF_DV[J];
! 2265: R = [];
! 2266: R = addpf([mtot(DYI,YI-YJ)],R);
! 2267: R = addpf([mtot(-DYJ,YI-YJ)],R);
! 2268: N = length(PF_V);
! 2269: BI = E[I]; BJ = E[J];
! 2270: for ( K = 1, D = 1; K < N; K++ )
! 2271: if ( K != I && K != J )
! 2272: D *= PF_DV[K]^E[K];
! 2273: S = [mtot(-1/(BJ+1)*DYI^(BI+1)*DYJ^(BJ+1)*D,1)];
! 2274: for ( K = 0; K <= BI; K++ )
! 2275: if ( -BI+BJ+2*K+2 )
! 2276: S = addpf([mtot(fac(BI)*fac(BJ)*(-BI+BJ+2*K+2)/fac(BI-K)/fac(BJ+K+2)*DYI^(BI-K)*DYJ^(BJ+K+2)*D,1)],S);
! 2277: return S;
! 2278: }
! 2279:
! 2280: def getchi(N)
! 2281: {
! 2282: CHI=[
! 2283: [5 ,11.07049769,1.145476226],
! 2284: [10 ,18.30703805,3.940299136],
! 2285: [20 ,31.41043284,10.85081139],
! 2286: [50 ,67.50480655,34.76425168],
! 2287: [100,124.3421134,77.92946517],
! 2288: [500,553.1268089,449.1467564],
! 2289: [1000, 1074.679449,927.594363]];
! 2290: for ( T = CHI; T != []; T = cdr(T) )
! 2291: if ( car(T)[0] == N ) return car(T);
! 2292: return [];
! 2293: }
! 2294:
! 2295: def message(S) {
! 2296: dp_gr_print(S);
! 2297: Print = S;
! 2298: }
! 2299: endmodule$
! 2300: end$
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to wishartd.rr CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / asir-contrib / testing / noro