Annotation of OpenXM/src/asir-contrib/packages/src/taka_diffop.rr, Revision 1.5
1.5 ! takayama 1: /* $OpenXM: OpenXM/src/asir-contrib/packages/src/taka_diffop.rr,v 1.4 2004/08/02 05:35:07 takayama Exp $ */
! 2: #define DEBUG
! 3:
1.1 takayama 4: #define ANK_SIZE 20
1.4 takayama 5:
6: #define XX x
7: #define LOG log(x)
8: /* #define LOG log(-x) */
9:
10: module taka_dr1;
11:
12: localf diff2 $
13: localf mult$
14: localf in$
15: localf division$
16: localf gcd$
17: localf act$
18: localf initANK$
19: localf ank$
20: localf thetapower2diff$
21: localf euler2diff$
22: localf diff2euler$
23: localf dp_poch $
24: localf int_roots $
25: localf root_multiplicity $
26: localf indicial $
27: localf kcoef0 $
28: localf test1 $
1.5 ! takayama 29: localf ord0 $
! 30: localf pfq_infty $
1.1 takayama 31:
32: extern ANK $
33: ANK = newmat(ANK_SIZE+1,ANK_SIZE+1) $
34:
35: /* --------------------------------------------------------------------
36: Ring of differential operators with rational functions coefficients
37: in one variables.
38: R_1 = K < x, dx >
39: ---------------------------------------------------------------------
40: */
41:
42: /* cf. d_gcd.rr */
1.4 takayama 43: def mult(F,G){
1.2 takayama 44: N=deg(nm(F),dx);
1.1 takayama 45: A=0;
46: for(K=0;K<=N;K++){
47: A=A+(1/fac(K))*diff2(F,dx,K)*diff2(G,x,K);
48: }
49: return red(A);
50: }
51: def diff2(F,G,K){
52: for(I=0;I<K;I++){
53: F=diff(F,G);
54: }
55: return F;
56: }
57:
58:
1.4 takayama 59: def in(F) {
1.1 takayama 60: Dn = dn(F); F=nm(F);
61: D = deg(F,dx);
62: C = coef(F,D,dx);
63: return( red(C*dx^D/Dn) );
64: }
65:
1.4 takayama 66: def division(F,G) {
1.1 takayama 67: Q = 0; R = F;
68: while ((R != 0) && (deg(nm(R),dx) >= deg(nm(G),dx))) {
1.4 takayama 69: D = red(in(R)/in(G));
1.1 takayama 70: Q = red(Q+D);
1.4 takayama 71: R = red(R-mult(D,G));
1.1 takayama 72: }
73: return([Q,R]);
74: }
75:
76:
1.4 takayama 77: def gcd(F,G) {
1.1 takayama 78: S = 1; T = 0; SS=0; TT=1;
1.2 takayama 79: if (deg(nm(F),dx) > deg(nm(G),dx)) {
1.1 takayama 80: Changed = 0;
81: }else {
82: T=F; F=G; G=T; Changed = 1;
83: }
84: do {
1.4 takayama 85: D=division(F,G);
1.1 takayama 86: Q = D[0]; R = D[1];
87: SSS = S-Q*SS; TTT = T-Q*TT;
88: F = G; G = R;
89: S = SS; T = TT; SS = SSS; TT = TTT;
90: } while ( R != 0);
91: if (Changed) {
92: return [F,T,S];
93: }else {
94: return [F,S,T];
95: }
1.3 takayama 96: }
97:
1.4 takayama 98: def act(L,F) {
1.3 takayama 99: D = dn(L);
100: L = nm(L);
101: N = deg(L,dx);
102: R = 0;
103: for (I=0; I<=N; I++) {
104: R += coef(L,I,dx)*diff2(F,x,I);
105: }
106: return red(R/D);
1.1 takayama 107: }
108:
109: /* --------------------------------------------------------
110: euler operators
111: --------------------------------------------------------
112: */
113:
114: def initANK() { /* Note, 2/22, 2003 */
115: extern ANK ;
116: ANK[1][1] = 1;
117: for (NN=2; NN<ANK_SIZE; NN++) {
118: for (KK=1; KK<=NN; KK++) {
119: if ((NN-1 < 1) || (KK > NN-1)) A = 0;
120: else A = (ANK[NN-1][KK])*KK;
121: if (KK-1 < 1) B = 0;
122: else B = ANK[NN-1][KK-1];
123: ANK[NN][KK] = A+B;
124: }
125: }
126: }
127:
128: def ank(N,K) {
129: extern ANK ;
130: if (!ANK[1][1]) {
131: initANK();
132: }
133: if ( N < 1 ) error("ank() out of index.");
134: if ( K > N ) return 0;
135: if ( K <= 0 ) return 0;
136: if (N < ANK_SIZE) return ANK[N][K];
137: return ank(N-1,K)*K+ank(N-1,K-1);
138: }
139:
140: /* Expand (X D)^J in X and D */
141: def thetapower2diff(J,X,D) {
142: if (J < 0) error("thetapower2diff() out of index.");
143: if (J == 0) return 1;
144: R = 0;
145: for (I = 1; I <= J; I++) {
146: R += ank(J,I)*X^I*D^I;
147: }
148: return R;
149: }
150:
151: def euler2diff(F,Xlist,Tlist,Dlist) {
152: DN =dn(F);
153: F = nm(F);
154: for (I=0; I < length(Xlist); I++) {
155: T = Tlist[I]; X = Xlist[I]; D = Dlist[I];
156: Deg = deg(F,T); R = 0;
157: for (J=0; J<=Deg; J++) {
158: C = coef(F,J,T);
159: R += C*thetapower2diff(J,X,D);
160: }
161: F = R;
162: }
163: return red(F/DN);
164: }
1.4 takayama 165:
166: /* -----------------------------------
167: diff2euler(x*dx+c,t) ==> [t+c,1]
168: diff2euler(dx+c,t) ==> [t+cx,1/x]
169: ----------------------------------
170: */
171: def diff2euler(F,T) {
172: DN = dn(F);
173: F = nm(F);
174: F = dp_ptod(F,[x,dx]);
175: R = 0;
176: while (F != 0) {
177: P=dp_hm(F);
178: F=dp_rest(F);
179: E = dp_etov(P);
180: R += dp_hc(P)*dp_vtoe(newvect(2,[E[0]-E[1],0]))*dp_poch(E[1]);
181: }
182: R = dp_dtop(R,[x,T]);
183: R = eval_str(rtostr(R));
184: DN2= dn(R);
185: R = nm(R);
186: /* factor out x */
187: M = root_multiplicity(R,x,0);
188: return [red(R/x^M),x^M/(DN*DN2)];
189: }
190:
191: def dp_poch(K) {
192: One = <<0,0>>;
193: T = <<0,1>>;
194: R = One;
195: for (I=0; I<K; I++) {
196: R = R*(T-I*One);
197: }
198: return R;
199: }
200:
201: /* int_roots(x^2*(x-1)^3*(x^2+1),x); --> [[1,3], [0,2]]
202: non-negative integral roots.
203: */
204: def int_roots(F,S) {
205: Root = [];
206: F = fctr(F);
207: F = cdr(F);
208: N = length(F);
209: for (I=0; I<N; I++) {
210: P = F[I][0]; /* factor */
211: Pm = F[I][1]; /* Multiplicity */
212: if (deg(P,S) == 1) {
213: P = -subst(P,S,0); /* Find root. */
214: if (type(P) <=1) {
215: if (P >= 0) Root = cons([P,Pm],Root);
216: }
217: }
218: }
219: return Root;
220: }
221:
222: /* root_multiplicity(x^2*(x-1),x,0); --> 2
223: multiplicity at x = 0.
224: */
225: def root_multiplicity(F,S,P) {
226: F = fctr(F);
227: F = cdr(F);
228: N = length(F);
229: for (I=0; I<N; I++) {
230: Pf = F[I][0]; /* factor */
231: Pm = F[I][1]; /* Multiplicity */
232: if (deg(Pf,S) == 1) {
233: Pf = -subst(Pf,S,P); /* Find root. */
234: if (Pf == 0) {
235: return Pm;
236: }
237: }
238: }
239: return 0;
240: }
241:
242: /* L is a differential operator in x and dx */
243: def indicial(L) {
244: /* indicial polynomial is a polynomial in x and s */
245: L = diff2euler(L,s)[0];
246: L0 = subst(L,x,0);
247: return [L0,-(L-L0)]; /* inidicial and the rest */
248: }
249:
250:
251: /* ---------------------------------------------------
252: Determine a_{k,i}
253: a_{k,i} x^k (log(x))^i, i=0, 1,..., m
254: L = [L0,LR] (euler op in s, output of indicial). L0-LR
255: Prev : series solution in x upto degree k-1
256: 0 ==> slow version.
1.5 ! takayama 257: DONT USE THE VARIABLE atmp_
1.4 takayama 258: ----------------------------------------------------- */
259: def kcoef0(K,M,L,Prev) {
260: L0 = L[0]; D_L0 = euler2diff(L0,[x],[s],[dx]);
261: LR = L[1]; D_LR = euler2diff(LR,[x],[s],[dx]);
262: V = [ ];
263: F = 0;
264: M0 = root_multiplicity(L0,s,K); /* exclude starting terms */
265: for (I=M0; I<=M; I++) {
1.5 ! takayama 266: V = cons(util_v("atmp",[I]),V);
! 267: F = F + util_v("atmp",[I])*(XX)^K*(LOG)^I;
1.4 takayama 268: }
269: V = reverse(V);
270: /* print([M0,M,V,F]); */
271: P0 = act(D_L0,F);
272: P1 = act(D_LR,Prev);
273: P0=red(subst(P0,LOG,y)/x^K);
1.5 ! takayama 274: /*
1.4 takayama 275: P1=coef(subst(P1,LOG,y),K,x);
1.5 ! takayama 276: */
! 277: Ptmp = subst(P1,LOG,y);
! 278: P1=red(coef(nm(Ptmp),K,x)/dn(Ptmp));
! 279:
! 280: #ifdef DEBUG
1.4 takayama 281: print([P0,P1]);
1.5 ! takayama 282: #endif
1.4 takayama 283:
284: /* compare coefficients of y=LOG */
285: Eq = [ ];
1.5 ! takayama 286: PPP = nm(P0-P1);
1.4 takayama 287: for (I=0; I<=M; I++) {
1.5 ! takayama 288: /* G = coef(P0,I,y)-coef(P1,I,y); */
! 289: G = coef(PPP,I,y);
1.4 takayama 290: if (G != 0) {
291: Eq = cons(G,Eq);
292: }
293: }
294: EqAns = poly_solve_linear(Eq,V);
1.5 ! takayama 295: #ifdef DEBUG
! 296: print(EqAns);
! 297: #endif
1.4 takayama 298: return base_replace(F,EqAns);
299: }
300:
301: /* testing kcoef0 */
302: def test1() {
303: L = euler2diff(s^2-x*(s+1/2)^2,[x],[s],[dx]);
304: LL = indicial(L);
305: F0 = 1;
306: print(sprintf("Indicial = %a",LL));
307: for (I=1; I<4; I++) {
308: Fnew = kcoef0(I,1,LL,F0);
309: print(sprintf("I=%a, %a",I,Fnew));
310: F0 = F0 + Fnew;
311: }
312: print("------------------------------");
313: F1 = LOG;
314: for (I=1; I<4; I++) {
315: Fnew = kcoef0(I,1,LL,F1);
316: print(sprintf("I=%a, %a",I,Fnew));
317: F1 = F1 + Fnew;
318: }
319: print("------------------------------");
320:
321: return [F0,F1];
322: /* cf. SST, p.23
323: [1+(1/4)*x+(9/64)*x,
324: log(x)+(-2*(1/2)^2+1 + 1/4*log(x))*x]
325: */
326: }
327:
1.5 ! takayama 328: /* order at X=0 of F*/
! 329: def ord0(F,X) {
! 330: F = nm(F);
! 331: N = deg(F,X);
! 332: for (I=0; I<=N; I++) {
! 333: if (coef(F,I,X) != 0) return I;
! 334: }
! 335: return -1;
! 336: }
! 337:
! 338: /*
! 339: Construct series solutions of p F q (A,B) up to O(1/x^N)
! 340: Use "CHECK" to check if a given series is an approximate sol.
! 341: order must be equal to the last arg (5 or 3).
! 342: taka_dr1.pfq_infty([1/2,1/2],[0],5);
! 343: taka_dr1.pfq_infty([a,b],[c],5);
! 344: taka_dr1.pfq_infty([a,a],[c],5);
! 345: taka_dr1.pfq_infty([a,a,b],[c1,c2],3);
! 346: taka_dr1.pfq_infty([a,a,b],[2,3],3);
! 347: taka_dr1.pfq_infty([1/2,1/2,-1/3],[2,3],5);
! 348: taka_dr1.pfq_infty([1/2,1/2,1/2,-3/2],[2,3,4],5);
! 349: */
! 350: /*&usage begin: taka_dr1.pfq_infty(A,B,N)
! 351: It computes the set of the canonical series solutions up to degree {N}
! 352: at z = infty for the generalized hypergeometric function p F q ({A},{B};z).
! 353: The variable 1/z is put x.
! 354: example: taka_dr1.pfq_infty([1/2,1/2,-3/2],[2,3],5);
! 355: */
! 356: def pfq_infty(A,B,N) {
! 357: P = length(A);
! 358: Q = length(B);
! 359: L0 = -s;
! 360: for (I=0; I<Q; I++) {
! 361: L0 = (-s+B[I]-1)*L0;
! 362: }
! 363: L1 = 1;
! 364: for (I=0; I<P; I++) {
! 365: L1 = (-s+A[I])*L1;
! 366: }
! 367: L = L1-x*L0; /* operator at infty */
! 368:
! 369: L1 = cdr(fctr(L1));
! 370: M = length(L1);
! 371: Roots = [ ];
! 372: LogMult = 0; /* Max exponent of log(x) */
! 373: for (I=0; I<M; I++) {
! 374: Roots = cons(-subst(L1[I][0],s,0)/coef(L1[I][0],1,s),Roots);
! 375: if (L1[I][1]-1 > LogMult) LogMult = L1[I][1];
! 376: }
! 377: Roots = reverse(Roots);
! 378: M = length(Roots);
! 379: Ans = [ ];
! 380: for (Ri=0; Ri<M; Ri++) {
! 381: /* x^E F(x). Le is the operator for F(x) */
! 382: E = Roots[Ri];
! 383: Le = subst(L,s,s+E);
! 384: LL = indicial(Le);
! 385: MM = root_multiplicity(LL[0],s,0);
! 386:
! 387: for (K=0; K<MM; K++) {
! 388: F0 = (LOG)^K;
! 389:
! 390: #ifdef DEBUG
! 391: print(sprintf("Indicial = %a",LL));
! 392: #endif
! 393:
! 394: for (I=1; I<N; I++) {
! 395: Fnew = kcoef0(I,LogMult,LL,F0);
! 396:
! 397: #ifdef DEBUG
! 398: print(sprintf("I=%a, %a",I,Fnew));
! 399: #endif
! 400:
! 401: F0 = F0 + Fnew;
! 402: }
! 403:
! 404: #ifdef DEBUG
! 405: LLL = euler2diff(Le,[x],[s],[dx]);
! 406: FFF = act(LLL,F0);
! 407: print(sprintf("order in x is %a",ord0(FFF,x)));
! 408: #endif
! 409:
! 410: print("------------------------------");
! 411:
! 412: #ifdef DEBUG
! 413: /* debug; */
! 414: #endif
! 415:
! 416: Ans = cons([x^E,F0],Ans);
! 417: }
! 418: }
! 419: return reverse(Ans);
1.4 takayama 420: }
1.5 ! takayama 421:
1.4 takayama 422:
423: endmodule;
1.1 takayama 424:
425: end$
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to taka_diffop.rr CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / asir-contrib / packages / src