Return to r-fd.rr CVS log

Up to [local] / OpenXM / src / asir-port / cgi

fd_hessian2(A,B,C,Xval) and r_fd_hessian2(A,B,C,Xval)
are added to cgi-asir-r-fd2.sh

/* $OpenXM: OpenXM/src/asir-port/cgi/r-fd.rr,v 1.11 2015/08/16 07:32:21 takayama Exp $ */
load("tk_fd.rr")$
import("tk_r.rr")$
import("oh_number.rr")$
import("test_hook.rr")$  /* To put overriden functions */

/* r_d2rat(0.3)  --> precision loss in truncation if not ctrl("bigfloat",1) */
ctrl("bigfloat",1)$
def r_d2rat(Y,Prec) {
  if ((type(Y) ==4)||(type(Y)==5)||(type(Y)==6)) return map(r_d2rat,Y,Prec);
  if ((type(Y)  == 1) && (ntype(Y) >= 1)) {
    S = rtostr(Y);  Y = "eval(("+S+")*exp(0));";
    /* print(Y); */
    Y = eval_str(Y);
    /* printf("Y=%a\n",Y); */
    /* return oh_number.rats(Y); */
    return rats2(Y | prec=Prec);  /* temporary  */
  }else return Y;
}
def r_ahvec(A,B,C,Y) { 
  if (type(getopt(prec))<0) Prec=20;
  else Prec=getopt(prec);
  Y = r_d2rat(Y,Prec);
  Ans=a_ahvec(A,B,C,Y);
  /*  Fans=map(rtostr,map(tk_fd.tk_number_rattofloat,Ans)); */
  Fans=map(deval,Ans);
  Fans = tk_r.asir2r_c(Fans);
  return Fans;
}

def a_ahvec(A,B,C,Y) { 
  R=tk_fd.ahvec_abc(A,B,C,Y|all=1);
  Gamma=R[1];
  Der=R[0];
  Z=R[2]*Gamma;
  Der2 = newvect(length(Der));
  for (I=0; I<length(Der); I++) Der2[I] = Der[I]*Gamma;
  Der2 = vtol(Der2);
  Ans=cons(Z,Der2);
  return(Ans);
}

/* temporary */
def rats2(X) {
  if (type(getopt(prec))<0) Prec=20;
  else Prec=getopt(prec);
  if (X == 0) return 0;
  Sign=1;
  if (X <0) {Sign=-1  ; X = -X;}
  Digit = number_floor(eval(log(X)/log(10)));  
  Num = number_floor((X/(10^Digit))*10^Prec);
  return Sign*(Num/(10^Prec))*(10^Digit);
}

def checkrats2() {
  for (I=0; I<10; I++) {  
     Sign=(-1)^(random()%2);
     X = eval(exp(0)*(random()/random())*10^(Sign*(random()%300))); /* 308 */
     printf("X=%a\n",X); 
     Y = rats2(X);
     printf("Y=%a\n",Y);
     if (number_abs(eval(Y*exp(0))/X-1) > 0.0000001) {
        printf("error: X = %a, Y=%a\n",X,Y);
     }
  }
}

def a_expect(A,B,C,Y) {
  E=tk_fd.expectation_abc(A,B,C,Y);
  return(E);
}
def r_expect(A,B,C,Y) {
  if (type(getopt(prec))<0) Prec=20;
  else Prec=getopt(prec);

  Y = r_d2rat(Y,Prec);
  E=a_expect(A,B,C,Y);
  Fans=map_deval(E);
  Fans = tk_r.asir2r_c(Fans);
  return Fans;
}

def r_ahmat(A,B,C,Y) { 
  if (type(getopt(prec))<0) Prec=20;
  else Prec=getopt(prec);
  Y = r_d2rat(Y,Prec);
  Ans=a_ahmat(A,B,C,Y);
  Fans=map_deval(Ans);
  Fans = tk_r.asir2r_c(Fans);
  return Fans;
}

def a_ahmat(A,B,C,Y) { 
  return(tk_fd.ahmat_abc(A,B,C,Y));
}

def r_log_ahmat(A,B,C,Y) { 
  if (type(getopt(prec))<0) Prec=20;
  else Prec=getopt(prec);
  Y = r_d2rat(Y,Prec);
  Ans=a_log_ahmat(A,B,C,Y);
  Fans=map_deval(Ans);
  Fans = tk_r.asir2r_c(Fans);
  return Fans;
}

def a_log_ahmat(A,B,C,Y) { 
  Ans=tk_fd.log_ahmat_abc(A,B,C,Y);
  return Ans;
}

def map_deval(L) {
  if (type(L) >=4) return(map(map_deval,L));
  return(deval(L));
}

def fd_hessian2(A,B,C,Xval) {
  H = tk_fd.fd_hessian2(A,B,C,Xval);
  return([H[0],matrix_matrix_to_list(H[1]),matrix_matrix_to_list(H[2])]);
}
def r_fd_hessian2(A,B,C,Xval) {
  H = tk_fd.fd_hessian2(A,B,C,Xval);
  return(tk_r.asir2r_c(H));
}
end$