Annotation of OpenXM/src/asir-contrib/packages/src/taka_adaptive.rr, Revision 1.4
1.4 ! takayama 1: /* $OpenXM: OpenXM/src/asir-contrib/packages/src/taka_adaptive.rr,v 1.3 2003/01/17 02:08:40 takayama Exp $ */
1.1 takayama 2:
3: /* adaptive_int(P,Q,Eps)
4: adaptive_barns(A,B,Z)
5: are defined.
6: */
7: /* see foo2(), foo3(), foo4() as for examples */
8:
9: /* Original from misc-2002/01/barns3.rr
10: It uses 9th order Newton-Cotes formula.
11: Reference, Mori, p.187, Iwanami */
12: /* #define DEBUG */
13: /* #define TEST_FOO */
14:
15: /* -------------- name replacement for asir-contrib ---------- */
16: /* While debugging, comment out the following definitions */
17: #define A_LevelMax Taka_A_LevelMax
18: #define A_P Taka_A_P
19: #define A_Q Taka_A_Q
20: #define Fvalue Taka_A_Fvalue
21: #define Fvalue2 Taka_A_Fvalue2
22: #define A_N Taka_A_N
23: #define A_A Taka_A_A
24: #define A_B Taka_A_B
25: #define A_C Taka_A_C
26: #define A_S1 Taka_A_S1
27: #define A_Z Taka_A_Z
28: #define NC_weight Taka_A_NC_weight
29: #define NC_n Taka_A_NC_n
30: #define A_C1 Taka_A_C1
31: #define A_aaa Taka_A_aaa
32: #define A_bbb Taka_A_bbb
33: #define A_ppp Taka_A_ppp
34: #define A_qqq Taka_A_qqq
35: #define kf(a) adaptive_kf(a)
36: #define kf2(a,b,c,z,s1,s2) taka_adaptive_kf2(a,b,c,z,s1,s2)
37: #define newton_cotes8(a,b) taka_adaptive_newton_cotes8(a,b)
38: #define adaptive_int0(a,b,eps,qold,steps,level) taka_adaptive_int0(a,b,eps,qold,steps,level)
39: #define barnsPFQ0(a,b,z) taka_adaptive_barnsPFQ0(a,b,z)
40:
41: #define graph(p,q,step) adaptive_graph(p,q,step)
42: #define barnsPFQ(a,b,z) adaptive_barnsPFQ(a,b,z)
43:
44:
45: /*
46: #define adaptive_int(p,q,eps) taka_adaptive_int(p,q,eps)
47: #define graph(p,q,step) taka_adaptive_graph(p,q,step)
48: */
49:
50: /* BUG!!! NOTE!!! Add space before $, otherwise cpp cannot handle. */
51: /* ------------------------ end of name replacement -----------------*/
52:
53: setprec(100) $
54:
55: #define INIT_N 8
56: #define BASE 2
57:
58: if (!Taka_A_Not_loaded) {
59: load("taka_pfp.rr") $
60: load("taka_plot.rr") $
61: pari(allocatemem, 10^7) $
62: Taka_A_Not_loaded =1 $
63: } $
64:
65: extern A_P,A_Q,Fvalue,Fvalue2,A_N,A_LevelMax $
66: extern A_A,A_B,A_C, A_S1,A_Z $
67: extern NC_weight,NC_n $
68: extern A_aaa, A_bbb $ /* Store parameters for p F q */
69: extern A_ppp, A_qqq $ /* A_ppp = p, A_qqq = q */
70:
71: /* Don't change these values during the execution of the program */
72: A_LevelMax = 10 $ /*&ja $B$I$3$^$G:YJ,$r$f$k$9$+(B? */
73: A_N = INIT_N*(BASE^A_LevelMax) $
74: Fvalue = newvect(A_N+1) $
75: Fvalue2 = newvect(A_N+1) $
76: NC_n = 14175 $
77: NC_weight = newvect(9) $
78: NC_weight[0] = 3956/NC_n $
79: NC_weight[1] = 23552/NC_n $
80: NC_weight[2] = -3712/NC_n $
81: NC_weight[3] = 41984/NC_n $
82: NC_weight[4] = -18160/NC_n $
83: NC_weight[5] = 41984/NC_n $
84: NC_weight[6] = -3712/NC_n $
85: NC_weight[7] = 23552/NC_n $
86: NC_weight[8] = 3956/NC_n $
87: /* -------------------------------------------------------------- */
88:
89: A_P = -3 $
90: A_Q = 20 $
91:
92: /*&ja $B%F%9%HMQ$NCM(B */
93: A_A=7/2 $
94: A_B=7/2 $
95: A_C=31/5 $
96: A_Z=1.3+1.8*@i $
97: A_S1= -0.3 $
98: /* expects -0.376954 - 0.282286 I */
99: extern A_C1 $
100: A_C1 = eval(1/(2*@pi))*pari(gamma,A_C)/(pari(gamma,A_A)*pari(gamma,A_B)) $
101: A_C1 *= eval(exp(A_S1*log(-A_Z))) $
102:
103:
104: /*&ja Step $B$O%0%i%U$r=q$/$?$a$NI}(B.
105: adaptive_kf() $B$N%0%i%U$rI=<((B.
106: */
107: def graph(P,Q,Step) {
108: T = [];
109: for (S2 = P; S2 < Q ; S2 += Step) {
110: K = kf(S2);
111: T = cons([S2,number_abs(K)],T);
112: }
113: T = reverse(T);
114: taka_plot_auto(T);
115: }
116:
117: /*&ja $B@QJ,$9$Y$-4X?t$NDj5A(B.
118: $B%5%s%W%k(B. 2F1
119: */
120: def kf(S2) {
121: return kf2(A_A,A_B,A_C,A_Z,A_S1,S2);
122: }
123:
124: /*&ja -re(A), -re(B), -re(C) < S1 $B$G$J$$$H$$$1$J$$(B. */
125: def kf2(A,B,C,Z,S1,S2) {
126: K = pari(gamma,A+S1+@i*S2)*pari(gamma,B+S1+@i*S2)*
127: pari(gamma,-S1-@i*S2)*eval(exp(@i*S2*log(-Z)));
128: K = K/pari(gamma,C+S1+@i*S2);
129: }
130:
131: /*&ja A0,B0,S0 $B$O@0?t(B
132: Hmin = (Q-P)/INIT_N*2^LevelMax;
133: P + A0*Hmin $B$+$i(B P + B0*Hmin $B$^$G@QJ,(B.
134: 8 $BEyJ,$7$F7W;;(B.
135: */
136: def newton_cotes8(A0,B0) {
137: S = 0.0;
138: S0 = (B0-A0)/8;
139: H = eval(S0*(A_Q-A_P)/A_N);
140: for (I=A0,J=0; I<=B0; I += S0, J++) {
141: if (! Fvalue2[I]) {
142: Fvalue[I] = kf(A_P+I*(A_Q-A_P)/A_N);
143: Fvalue2[I] = 1;
144: }
145: S += H*Fvalue[I]*NC_weight[J];
146: }
147: return S;
148: }
149:
150: def adaptive_int0(A,B,Eps,Qold,Steps, Level) {
151: #ifdef DEBUG
152: print("adaptive_int0 (A,B,Eps,Qold,Steps,Level):",0);
153: print([A,B,Eps,Qold,Steps,Level]);
154: #endif
155: if (Level > A_LevelMax) {
156: print("Warning: Max depth is reached at ",0); print([A,B]);
157: print("Check the shape of the graph by adaptive_graph(P,Q,Step)");
158: return Qold;
159: }
160: NewSteps = Steps/BASE;
161:
162: NewEps = eval(Eps/BASE);
163:
164: S1 = newton_cotes8(A,A+(B-A)/2);
165: S2 = newton_cotes8(A+(B-A)/2,B);
166:
167: /* Mori, p.192, (12.13) */
168: if (eval(number_abs(S1+S2-Qold)/1023) < eval(Eps/(A_Q-A_P))) {
169: return S1+S2;
170: }
171: S1 = adaptive_int0(A,A+(B-A)/2,NewEps,S1,NewSteps,Level+1);
172: S2 = adaptive_int0(A+(B-A)/2,B,NewEps,S2,NewSteps,Level+1);
173: return S1+S2;
174: }
175:
176:
177: /*&usage begin: adaptive_int(P,Q,Eps)
178: It numerically evaluates the integral of adaptive_kf(x) over
179: the interval [{P},{Q}].
180: It uses an adaptive integration scheme based on 8th order Newton-Cotes
181: formula.
182: ref: adaptive_graph
183: */
184: def adaptive_int(P,Q,Eps) {
185: Steps = A_N/INIT_N;
186: for (I=0; I<=A_N; I++) Fvalue2[I]=0;
187: S0 = newton_cotes8(0,A_N);
188: S = adaptive_int0(0,A_N,Eps,S0,Steps/BASE, 1);
189: return S;
190: }
191:
192: def foo() {
193: /* S = adaptive_int(A_P,A_Q,0.0000001); */
194: tstart;
195: S = adaptive_int(A_P,A_Q,0.0001);
196: tstop;
197: return A_C1*S;
198: }
199:
200:
201: /* ------------- barns integral representation for hypergeometric function --- */
202:
203: /*&ja barns $B@QJ,$N3K$r7W;;$9$k$?$a$NDj?t$N@_Dj(B, $B$*$h$SDj?tIt$r@_Dj$7La$9(B. */
204: def barnsPFQ0(A,B,Z) {
205: A_Z = Z;
206: A_ppp = length(A);
207: A_qqq = length(B);
208: A_aaa = newvect(A_ppp);
209: A_bbb = newvect(A_qqq);
210: for (I=0; I<A_ppp; I++) {
211: A_aaa[I] = A[I];
212: if (-number_real_part(A[I]) >= A_S1) {
213: print("I="+rtostr(I)+","+rtostr(A[I]));
214: error("Out of valid parameters of A for barns integral.");
215: }
216: }
217: for (I=0; I<A_qqq; I++) {
218: A_bbb[I] = B[I];
219: if (-number_real_part(B[I]) >= A_S1) {
220: print("I="+rtostr(I)+","+rtostr(A[I]));
221: error("Out of valid parameters of B for barns integral.");
222: }
223: }
224:
225: A_C1 = eval(exp(A_S1*log(-A_Z))/(2*@pi));
226: for (I=0; I<A_qqq; I++) {
227: A_C1 *= pari(gamma,A_bbb[I]);
228: }
229: T = 1;
230: for (I=0; I<A_ppp; I++) {
231: T *= pari(gamma,A_aaa[I]);
232: }
233: A_C1 = A_C1/T;
234: return A_C1;
235: }
236:
237: #ifndef TEST_FOO
238: def kf(S2) {
239: K1 = 1;
240: for (I=0; I<A_ppp; I++) {
241: K1 *= pari(gamma,A_aaa[I]+A_S1+@i*S2);
242: }
243: K1 *= pari(gamma,-A_S1-@i*S2)*eval(exp(@i*S2*log(-A_Z)));
244:
245: K2 = 1;
246: for (I=0; I<A_qqq; I++) {
247: K2 *= pari(gamma,A_bbb[I]+A_S1+@i*S2);
248: }
249: return K1/K2;
250: }
251: #endif
252:
253:
254: def barnsPFQ(A,B,Z) {
255: Opt = getopt();
256: Eps = 0.0001;
257: Graph = 0;
258: for (I=0; I<length(Opt); I++) {
259: if (Opt[I][0] == "interval") {
260: A_P = Opt[I][1][0];
261: A_Q = Opt[I][1][1];
262: }else if (Opt[I][0] == "error") {
263: Eps = Opt[I][1];
264: }else if (Opt[I][0] == "graph") {
265: Graph = Opt[I][1];
1.4 ! takayama 266: }else if (Opt[I][0] == "real_part_of_path") {
! 267: A_S1 = Opt[I][1];
! 268: if (A_S1 >= 0) error("Real part of the path must be negative.");
1.1 takayama 269: }else{
270: print(Opt[I]);
271: error("Unknown option name.");
272: }
273: }
274: barnsPFQ0(A,B,Z);
275: /* test.
276: def kf(S) { return S^2; }
277: graph(A_P,A_Q,0.5); */
278:
279: #ifdef DEBUG
280: graph(A_P,A_Q,0.5);
281: #else
282: if (Graph) graph(A_P,A_Q,0.5);
283: #endif
284: if ((number_abs(kf(A_P)) > Eps) || (number_abs(kf(A_Q)) > Eps)) {
285: print("Warning: integrant does not decay to zero at the boundary of ",0);
286: print([A_P,A_Q]);
1.2 takayama 287: print("kf(P)="+rtostr(kf(A_P)));
288: print("kf(Q)="+rtostr(kf(A_Q)));
1.1 takayama 289: }
290: return A_C1*adaptive_int(A_P,A_Q,Eps);
291:
292: }
293:
294: /* ------------------- test suits ---------------- */
295: def foo2() {
296: /* return barnsPFQ([7/2,7/2],[31/5],1.3+@i*1.8);
297: It causes macro error.
298: */
299: A = [7/2,7/2]; B=[31/5];
300: return barnsPFQ(A,B,1.3+@i*1.8);
301: }
302:
303: def foo3() {
304: A = [7/2,7/2]; B=[31/5];
305: for (Z = 1.3+@i*1.8; number_real_part(Z) < 30; Z += 5) {
306: tstart;
307: print("Z="+rtostr(Z)+", F=",0);
308: print(barnsPFQ(A,B,Z));
309: tstop;
310: }
311: /*
312: For[z = 1.3+I*1.8, Re[z] < 30, z += 5,
313: f = Timing[N[Hypergeometric2F1[7/2,7/2,31/5,z]]];
314: Print["z=",z,", F=",f]];
315:
316: z=6.3 + 1.8 I, F={0. Second, 0.0354872 + 0.0544005 I}
317: change interval to interval=[-3,30]
318: adaptive_barnsPFQ([7/2,7/2],[31/5],6.3+1.8*@i | graph=1, interval=[-3,30]);
319: (0.03597211652913092260+0.05415889075951361488*@i)
320: z=11.3 + 1.8 I, F={0.01 Second, -0.0041952 + 0.0113084 I}
321: z=16.3 + 1.8 I, F={0. Second, -0.00273279 + 0.00272846 I}
322: z=21.3 + 1.8 I, F={0. Second, -0.00142342 + 0.000830411 I}
323: z=26.3 + 1.8 I, F={0.01 Second, -0.000784199 + 0.000290081 I}
324: */
325: }
326:
327: def foo4() {
328: A = [7/2,7/2,8/3]; B=[31/5,36/7]; Z = 1.3+@i*1.8;
329: print(barnsPFQ(A,B,Z | graph=1));
330: /* (-0.4896050291590236027+0.7354763340691728624*@i
331: N[HypergeometricPFQ[{7/2,7/2,8/3},{31/5,36/7},1.3+I*1.8]]
332: -0.491559 + 0.737438 I
333: */
334: A = [7/2,7/2,8/3]; B=[31/5,36/7]; Z = 6.3+@i*1.8;
335: print(barnsPFQ(A,B,Z | graph=1));
336: /* (0.03557770119281511709-0.3947142140534400708*@i)
337: N[HypergeometricPFQ[{7/2,7/2,8/3},{31/5,36/7},6.3+I*1.8]]
338: 0.0356712 - 0.394754 I
339: */
340: }
341:
342: end $
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to taka_adaptive.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