Return to taka_adaptive.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

Fixed a bug on the "if" statement at the toplevel.

/* $OpenXM: OpenXM/src/asir-contrib/packages/src/taka_adaptive.rr,v 1.6 2004/04/24 06:35:57 takayama Exp $  */

/* adaptive_int(P,Q,Eps)
   adaptive_barns(A,B,Z)
   are defined.
*/
/* see foo2(), foo3(), foo4() as for examples */

/* Original from misc-2002/01/barns3.rr   
   It uses 9th order Newton-Cotes formula.
   Reference,  Mori, p.187, Iwanami */
/* #define DEBUG */
/* #define TEST_FOO */

/* -------------- name replacement for asir-contrib ---------- */
/* While debugging, comment out the following definitions */
#define A_LevelMax Taka_A_LevelMax
#define A_P Taka_A_P
#define A_Q Taka_A_Q
#define Fvalue Taka_A_Fvalue
#define Fvalue2 Taka_A_Fvalue2
#define A_N Taka_A_N
#define A_A Taka_A_A
#define A_B Taka_A_B
#define A_C Taka_A_C
#define A_S1 Taka_A_S1
#define A_Z Taka_A_Z
#define NC_weight Taka_A_NC_weight
#define NC_n Taka_A_NC_n
#define A_C1 Taka_A_C1
#define A_aaa Taka_A_aaa
#define A_bbb Taka_A_bbb
#define A_ppp Taka_A_ppp
#define A_qqq Taka_A_qqq
#define kf(a)  adaptive_kf(a)
#define kf2(a,b,c,z,s1,s2) taka_adaptive_kf2(a,b,c,z,s1,s2)
#define newton_cotes8(a,b) taka_adaptive_newton_cotes8(a,b)
#define adaptive_int0(a,b,eps,qold,steps,level) taka_adaptive_int0(a,b,eps,qold,steps,level)
#define barnsPFQ0(a,b,z)         taka_adaptive_barnsPFQ0(a,b,z)

#define graph(p,q,step)          adaptive_graph(p,q,step)
#define barnsPFQ(a,b,z)          adaptive_barnsPFQ(a,b,z)


/* 
 #define adaptive_int(p,q,eps) taka_adaptive_int(p,q,eps)
 #define graph(p,q,step) taka_adaptive_graph(p,q,step)
*/

/* BUG!!! NOTE!!! Add space before $, otherwise cpp cannot handle. */
/* ------------------------ end of name replacement -----------------*/

setprec(100) $

#define INIT_N       8 
#define BASE         2 

if (!Taka_A_Not_loaded)  {
  load("taka_pfp.rr") $
  load("taka_plot.rr") $
  pari(allocatemem, 10^7) $   
  Taka_A_Not_loaded =1 $
} else { } $

extern A_P,A_Q,Fvalue,Fvalue2,A_N,A_LevelMax $
extern A_A,A_B,A_C, A_S1,A_Z $
extern NC_weight,NC_n $
extern A_aaa, A_bbb $    /* Store parameters for p F q */
extern A_ppp, A_qqq $    /* A_ppp = p, A_qqq = q */

/* Don't change these values during the execution of the program */
A_LevelMax = 10 $  /*&ja $B$I$3$^$G:YJ,$r$f$k$9$+(B? */
A_N = INIT_N*(BASE^A_LevelMax) $
Fvalue = newvect(A_N+1) $
Fvalue2 = newvect(A_N+1) $
NC_n = 14175 $
NC_weight = newvect(9) $
NC_weight[0] = 3956/NC_n $
NC_weight[1] = 23552/NC_n $
NC_weight[2] = -3712/NC_n $
NC_weight[3] = 41984/NC_n $
NC_weight[4] = -18160/NC_n $
NC_weight[5] = 41984/NC_n $
NC_weight[6] = -3712/NC_n $
NC_weight[7] = 23552/NC_n $
NC_weight[8] = 3956/NC_n $
/* -------------------------------------------------------------- */

A_P = -3 $
A_Q =  20 $

/*&ja $B%F%9%HMQ$NCM(B */
A_A=7/2 $
A_B=7/2 $
A_C=31/5 $
A_Z=1.3+1.8*@i $
A_S1= -0.3 $
/* expects  -0.376954 - 0.282286 I */
extern A_C1 $
A_C1 = eval(1/(2*@pi))*pari(gamma,A_C)/(pari(gamma,A_A)*pari(gamma,A_B)) $
A_C1 *= eval(exp(A_S1*log(-A_Z))) $


/*&ja  Step $B$O%0%i%U$r=q$/$?$a$NI}(B.
       adaptive_kf() $B$N%0%i%U$rI=<((B.
*/
def graph(P,Q,Step) {
  T = [];
  for (S2 = P; S2 < Q ; S2 += Step) {
    K = kf(S2);
    T = cons([S2,number_abs(K)],T);
  }
  T = reverse(T);
  taka_plot_auto(T);
}

/*&ja  $B@QJ,$9$Y$-4X?t$NDj5A(B.
       $B%5%s%W%k(B.  2F1
*/
def kf(S2) {
  return kf2(A_A,A_B,A_C,A_Z,A_S1,S2);
}

/*&ja  -re(A), -re(B), -re(C) < S1 $B$G$J$$$H$$$1$J$$(B.  */
def kf2(A,B,C,Z,S1,S2) {
     K = pari(gamma,A+S1+@i*S2)*pari(gamma,B+S1+@i*S2)*
         pari(gamma,-S1-@i*S2)*eval(exp(@i*S2*log(-Z)));
     K = K/pari(gamma,C+S1+@i*S2);
}

/*&ja A0,B0,S0 $B$O@0?t(B
     Hmin = (Q-P)/INIT_N*2^LevelMax;
     P + A0*Hmin $B$+$i(B  P + B0*Hmin $B$^$G@QJ,(B.
     8 $BEyJ,$7$F7W;;(B.
*/
def newton_cotes8(A0,B0) {
  S = 0.0;
  S0 = (B0-A0)/8;
  H = eval(S0*(A_Q-A_P)/A_N);  
  for (I=A0,J=0; I<=B0; I += S0, J++) {
    if (! Fvalue2[I]) {
      Fvalue[I] = kf(A_P+I*(A_Q-A_P)/A_N);
      Fvalue2[I] = 1;
    }
    S += H*Fvalue[I]*NC_weight[J];
  }
  return S;
}

def adaptive_int0(A,B,Eps,Qold,Steps, Level) {
#ifdef DEBUG
  print("adaptive_int0 (A,B,Eps,Qold,Steps,Level):",0);
  print([A,B,Eps,Qold,Steps,Level]);
#endif
  if (Level > A_LevelMax) {
    print("Warning: Max depth is reached at ",0); print([A,B]);
    print("Check the shape of the graph by adaptive_graph(P,Q,Step)");
    return Qold;
  }
  NewSteps = Steps/BASE;

  NewEps = eval(Eps/BASE);

  S1 = newton_cotes8(A,A+(B-A)/2);
  S2 = newton_cotes8(A+(B-A)/2,B);

  /* Mori, p.192,  (12.13) */
  if (eval(number_abs(S1+S2-Qold)/1023) < eval(Eps/(A_Q-A_P))) {
    return S1+S2;
  }
  S1 = adaptive_int0(A,A+(B-A)/2,NewEps,S1,NewSteps,Level+1);
  S2 = adaptive_int0(A+(B-A)/2,B,NewEps,S2,NewSteps,Level+1);
  return S1+S2;
}


/*&usage begin: adaptive_int(P,Q,Eps)
  It numerically evaluates the integral of adaptive_kf(x) over 
  the interval [{P},{Q}].
  It uses an adaptive integration scheme based on 8th order Newton-Cotes 
  formula.
  ref: adaptive_graph
*/
def adaptive_int(P,Q,Eps) {
  A_P = P; A_Q = Q;  /* set the global variable */
  Steps = A_N/INIT_N;
  for (I=0; I<=A_N; I++) Fvalue2[I]=0;
  S0 = newton_cotes8(0,A_N);
  S = adaptive_int0(0,A_N,Eps,S0,Steps/BASE, 1);
  return S;
}

def foo() {
  /* S = adaptive_int(A_P,A_Q,0.0000001); */
  tstart;
  S = adaptive_int(A_P,A_Q,0.0001);
  tstop;
  return A_C1*S;  
}


/* ------------- barns integral representation for hypergeometric function --- */

/*&ja  barns $B@QJ,$N3K$r7W;;$9$k$?$a$NDj?t$N@_Dj(B, $B$*$h$SDj?tIt$r@_Dj$7La$9(B.  */
def barnsPFQ0(A,B,Z) {
  A_Z = Z;
  A_ppp = length(A);
  A_qqq = length(B);
  A_aaa = newvect(A_ppp);
  A_bbb = newvect(A_qqq);
  for (I=0; I<A_ppp; I++) {
    A_aaa[I] = A[I];
    if (-number_real_part(A[I]) >= A_S1) {
       print("I="+rtostr(I)+","+rtostr(A[I]));
       error("Out of valid parameters of A for barns integral.");
    }
  }
  for (I=0; I<A_qqq; I++) {
    A_bbb[I] = B[I];
    if (-number_real_part(B[I]) >= A_S1) {
       print("I="+rtostr(I)+","+rtostr(A[I]));
       error("Out of valid parameters of B for barns integral.");
    }
  }

  A_C1 = eval(exp(A_S1*log(-A_Z))/(2*@pi));
  for (I=0; I<A_qqq; I++) {
    A_C1 *= pari(gamma,A_bbb[I]);
  }
  T = 1;
  for (I=0; I<A_ppp; I++) {
    T *= pari(gamma,A_aaa[I]);
  }
  A_C1 = A_C1/T;
  return A_C1;
}

#ifndef TEST_FOO
def kf(S2) {
  K1 = 1;
  for (I=0; I<A_ppp; I++) {
      K1 *= pari(gamma,A_aaa[I]+A_S1+@i*S2);
  }
  K1 *= pari(gamma,-A_S1-@i*S2)*eval(exp(@i*S2*log(-A_Z)));

  K2 = 1;
  for (I=0; I<A_qqq; I++) {
    K2 *= pari(gamma,A_bbb[I]+A_S1+@i*S2);
  }
  return K1/K2;
}
#endif


def barnsPFQ(A,B,Z) {
  Opt = getopt();
  Eps = 0.0001;
  Graph = 0;
  for (I=0; I<length(Opt); I++) {
    if (Opt[I][0] == "interval") { 
       A_P = Opt[I][1][0];
       A_Q = Opt[I][1][1];
    }else if (Opt[I][0] == "error") {
       Eps = Opt[I][1];
    }else if (Opt[I][0] == "graph") {
       Graph = Opt[I][1];
    }else if (Opt[I][0] == "real_part_of_path") {
       A_S1 = Opt[I][1];
       if (A_S1 >= 0) error("Real part of the path must be negative.");
    }else{
      print(Opt[I]); 
      error("Unknown option name.");
    }
  }
  barnsPFQ0(A,B,Z);
  /* test.
     def kf(S) { return S^2; }
     graph(A_P,A_Q,0.5); */

#ifdef DEBUG
  graph(A_P,A_Q,0.5);
#else
  if (Graph) graph(A_P,A_Q,0.5);
#endif
  if ((number_abs(kf(A_P)) > Eps) || (number_abs(kf(A_Q)) > Eps)) {
    print("Warning: integrant does not decay to zero at the boundary of ",0);
    print([A_P,A_Q]);
    print("kf(P)="+rtostr(kf(A_P)));
    print("kf(Q)="+rtostr(kf(A_Q)));
  }
  return A_C1*adaptive_int(A_P,A_Q,Eps);

}

/* -------------------  test suits ---------------- */
def foo2() {
 /*  return barnsPFQ([7/2,7/2],[31/5],1.3+@i*1.8);
     It causes macro error.
  */
  A = [7/2,7/2]; B=[31/5];
  return barnsPFQ(A,B,1.3+@i*1.8);
}

def foo3() {
  A = [7/2,7/2]; B=[31/5];
  for (Z = 1.3+@i*1.8; number_real_part(Z) < 30; Z += 5) {
    tstart;
    print("Z="+rtostr(Z)+", F=",0);
    print(barnsPFQ(A,B,Z));
    tstop;
  }
/*
  For[z = 1.3+I*1.8, Re[z] < 30, z += 5,
    f = Timing[N[Hypergeometric2F1[7/2,7/2,31/5,z]]];
    Print["z=",z,", F=",f]];

  z=6.3 + 1.8 I, F={0. Second, 0.0354872 + 0.0544005 I}
    change interval to interval=[-3,30]
    adaptive_barnsPFQ([7/2,7/2],[31/5],6.3+1.8*@i | graph=1, interval=[-3,30]);
    (0.03597211652913092260+0.05415889075951361488*@i)
  z=11.3 + 1.8 I, F={0.01 Second, -0.0041952 + 0.0113084 I}
  z=16.3 + 1.8 I, F={0. Second, -0.00273279 + 0.00272846 I}
  z=21.3 + 1.8 I, F={0. Second, -0.00142342 + 0.000830411 I}
  z=26.3 + 1.8 I, F={0.01 Second, -0.000784199 + 0.000290081 I}
*/
}

def foo4() {
  A = [7/2,7/2,8/3]; B=[31/5,36/7]; Z = 1.3+@i*1.8;
  print(barnsPFQ(A,B,Z | graph=1));
  /* (-0.4896050291590236027+0.7354763340691728624*@i 
     N[HypergeometricPFQ[{7/2,7/2,8/3},{31/5,36/7},1.3+I*1.8]]
     -0.491559 + 0.737438 I
  */
  A = [7/2,7/2,8/3]; B=[31/5,36/7]; Z = 6.3+@i*1.8;
  print(barnsPFQ(A,B,Z | graph=1));
  /* (0.03557770119281511709-0.3947142140534400708*@i)
     N[HypergeometricPFQ[{7/2,7/2,8/3},{31/5,36/7},6.3+I*1.8]]
     0.0356712 - 0.394754 I
  */
}

end $