/* $Id: test1-tr.rr,v 1.16 2005/04/21 08:16:26 taka Exp $ */
/* $OpenXM: OpenXM/src/asir-contrib/testing/test1-tr.rr,v 1.6 2005/04/21 10:54:49 takayama Exp $ */

load("tr.rr")$




def test0() {
  A = quotetolist(quote(1+sin(x)+sin(3*@pi)*sin(0)));
  P = quotetolist(quote(sin(pn("x")*@pi)));
  Q = ["qt_sin_int","x"];
  print(A);
  print(P);
  print(Q);
  print("----------------");
  print(tr_match0(A,P));
  A2 = quotetolist(quote(sin(2*@pi)));
  print(tr_match0(A2,P));
  print("----------------");
  print("---- tr_make_binding --------");
  print(tr_make_binding(A2,P));
  print("-----tr_rp -------------");
  R=tr_rp(A,P,Q);
  print("--------------------");
  print(R);
  print("--------------------");
  return quote_input_form_quote_list(R);
}

def test1()  {
  Rule1=[quote(sin(pn("x")*@pi)),["qt_sin_int","x"]]; /* sin(整数*@pi) --> 0 */
  Rule2=[quote(0*pn("y")),       ["qt_zero"]];       /* 0*any --> 0 */
  Rule3=[quote(pn("y")*0),       ["qt_zero"]];       /* any*0 --> 0 */
  Rule4=[quote(pn("y")+0),       ["qt_id","y"]];       /* any+0 --> any */
  Rule5=[quote(0+pn("y")),       ["qt_id","y"]];       /* 0+any --> any */
  Rule6=[quote(sin(0)),          ["qt_zero"]];       /* sin(0) --> 0 */
  R0 = quote(1+sin(sin(2*@pi)*sin(@pi/2))+sin(5*@pi));
  print(print_input_form(R0));
  R=tr_apply_rule1(R0,Rule1[0],Rule1[1]);
  print(print_input_form(R));
  R=tr_apply_rule1(R,Rule2[0],Rule2[1]);
  print(print_input_form(R));
  R=tr_apply_rule1(R,Rule4[0],Rule4[1]);
  print(print_input_form(R));
  R=tr_apply_rule1(R,Rule6[0],Rule6[1]);
  print(print_input_form(R));
  R=tr_apply_rule1(R,Rule4[0],Rule4[1]);
  print(print_input_form(R));
  return R;
}


/* 不定積分計算の例  
    c x^n の和の不定積分 (c は x に依存せず)
   いろいろ 問題点あり:  たとえば c が 無いときの処理できず.
*/

/* 右辺関数.  c x^n の不定積分 (c は x に依存せず)
   Todo: 右辺関数を容易に書く方法.
*/
def r_integral0(C,N) {
  NN = eval_str(quote_input_form_quote_list(quotetolist(N)));
  CC = quote_input_form_quote_list(quotetolist(C));
  if (NN == -1) {
     R = "quote("+CC+"*log(x))";
  }else{
     R = "quote("+CC+"/"+rtostr(NN+1)+"*x^"+rtostr(NN+1)+")";
  }
  print("r_integral0:",0);print(R);
  R = eval_str(R);
  return quotetolist(R);
}
/* 右辺関数 積分の線型性 */
def r_int_linear(F,G) {
  FF = quote_input_form_quote_list(quotetolist(F));
  GG = quote_input_form_quote_list(quotetolist(G));
  R = "quote(integral("+FF+")+integral("+GG+"))";
  print("r_int_linear:",0);print(R);
  R = eval_str(R);
  return quotetolist(R);
}
def test3() {
  R0 = quote(1+integral(2*x^(-1)+2*x^2));
  return test3a(R0);
}
def test3a(R0)  {
  Rules=[
     /* c*x^n --> (c/(n+1))*x^(n+1) or c*log(x) */
     [quote(integral(pn("c")*x^pn("n"))),["r_integral0","c","n"]],
     [quote(integral(pn("f")+pn("g"))),  ["r_int_linear","f","g"]]
  ];
  print("Input=",0); print(print_input_form(R0));
  N = length(Rules);
  R = R0;
  for (J=0; J<3; J++) {  /* Todo: フラグがないので, とりあえず 3 回 */
    for (I=0; I<N; I++) {
      print(print_input_form(R));
      R=tr_apply_rule1(R,Rules[I][0],Rules[I][1]);
    }
  }
  return R;
}


/* 関数のマッチ.  N[] 相当.  test4(). */
/*  quote(nn(pn(f),qt_is_function(f))); は不要. qt_map_arg が処理 */
def test4() {
  Rule=[quote(nn(pn(f))),[qt_map_arg,nn,f]];
  R0 = quote(nn(sin(1/2)*cos(1/3)));
  print(print_input_form(R0));
  R=tr_apply_rule1(R0,Rule[0],Rule[1]);
  return R;
}

/* tr_apply_or_rule の試作 */

/* Flag 付き の tr_rp. 属性がないのでこれでやる. */
def tr_rp_flag(F,P,Q) {
  Flag = 0; 
  dprint0("tr_rp, F="); dprint(F);
  dprint0("tr_rp, P="); dprint(P);
  dprint0("tr_rp, Q="); dprint(Q);
  if (tr_match0(F,P)) {
     BindTable = tr_make_binding(F,P);
     dprint0("BindTable="); dprint(BindTable);
     return [1,tr_apply_function0(Q,BindTable)];
  }
  if (type(F) != 4) return F;
  Node = qt_node(F);
  N = qt_nchild(F);
  Ans = Node;
  for (I=0; I<N; I++) {
    T = tr_rp_flag(qt_child(F,I),P,Q);
    if (T[0] == 1) Flag = 1;
    Ans = append(Ans,[T[1]]); 
  }
  return [Flag,Ans];
}

extern Debug2$
Debug2=0$
/* 書き換え flag 付きの tr_apply_rule_flag */
def tr_apply_rule1_flag(Obj,L,R) {
  Flag = 0;
  if (Debug2) 
   print("--------  start of tr_apply_rule1_flag ------------ ");
  if (Debug2) print(print_input_form(Obj));
  Obj = quotetolist(Obj);
  L = quotetolist(L);
  R = tr_rp_flag(Obj,L,R);
  Flag=R[0]; R=R[1];
  if (type(R) == 17) R=quotetolist(R);
  RR = "quote("+listtoquote_str(R)+")";
  if (Debug2) {print("==> "+RR+"  by  "); print(listtoquote_str(L));}
  if (Debug2) print("--------  end of tr_apply_rule1_flag ------------ ");
  return [Flag,eval_str(RR)];  
}

def tr_apply_or_rules(Q,R) {
  Flag = 1;
  N = length(R);
  while (Flag) {
   Flag = 0;
   for (I=0; I<N; I++) {  
     Q = tr_apply_rule1_flag(Q,R[I][0],R[I][1]);
     if (Q[0]) {
       Flag = 1;
       dprint("Applied the rule "+rtostr(I));
     }
     Q = Q[1];
   }
  }
  return Q;
}
def test5() {
  Rule1=[quote(sin(pn(x)*@pi)),[qt_sin_int,x]]; /* sin(整数*@pi) --> 0 */
  Rule2=[quote(0*pn(y)),       [qt_zero]];       /* 0*any --> 0 */
  Rule3=[quote(pn(y)*0),       [qt_zero]];       /* any*0 --> 0 */
  Rule4=[quote(pn(y)+0),       [qt_id,y]];       /* any+0 --> any */
  Rule5=[quote(0+pn(y)),       [qt_id,y]];       /* 0+any --> any */
  Rule6=[quote(sin(0)),          [qt_zero]];       /* sin(0) --> 0 */
  R0 = quote(1+sin(sin(2*@pi)*sin(@pi/2))+sin(5*@pi));
  print(print_input_form(R0));
  R=tr_apply_rule1_flag(R0,Rule1[0],Rule1[1]);
  print([R[0],print_input_form(R[1])]);
  R=tr_apply_or_rules(R0,[Rule1,Rule2,Rule3,Rule4,Rule5,Rule6]);
  return R;
}

def qt_one() {
  return quote(1);
}
def tr_simp_sin(R0) {
  Rule1=[quote(sin(pn(x)*@pi)),[qt_sin_int,x]]; /* sin(整数*@pi) --> 0 */
  Rule2=[quote(0*pn(y)),       [qt_zero]];       /* 0*any --> 0 */
  Rule3=[quote(pn(y)*0),       [qt_zero]];       /* any*0 --> 0 */
  Rule4=[quote(pn(y)+0),       [qt_id,y]];       /* any+0 --> any */
  Rule5=[quote(0+pn(y)),       [qt_id,y]];       /* 0+any --> any */
  Rule6=[quote(sin(0)),        [qt_zero]];       /* sin(0) --> 0 */
  Rule7=[quote(cos(0)),        [qt_one]];         /* cos(0) --> 1 */
  /* print(print_input_form(R0)); */
  R=tr_apply_rule1_flag(R0,Rule1[0],Rule1[1]);
  /* print([R[0],print_input_form(R[1])]); */
  R=tr_apply_or_rules(R0,[Rule1,Rule2,Rule3,Rule4,Rule5,Rule6,Rule7]);
  return R;
}

/* 0+any, 0*any になる quote を 0 にする. 必須. cf. taka_series.expand1 */
def tr_simp_zero(R0) {
  Rule1=[quote(0*pn(y)),       [qt_zero]];       /* 0*any --> 0 */
  Rule2=[quote(pn(y)*0),       [qt_zero]];       /* any*0 --> 0 */
  Rule3=[quote(0/pn(y)),       [qt_zero]];       /* 0/any --> 0 */
  Rule4=[quote(pn(y)+0),       [qt_id,y]];       /* any+0 --> any */
  Rule5=[quote(0+pn(y)),       [qt_id,y]];       /* 0+any --> any */
  Rule6=[quote(-0),            [qt_zero,y]];       /* -0 --> 0 */
  R=tr_apply_or_rules(R0,[Rule1,Rule2,Rule3,Rule4,Rule5, Rule6]);
  return R;
}

/* 微分環の計算 */
/* x に依存してるか?  u, u_0, u_1, u_2, ... は x に依存してる.*/
def to_quote(L) {
  return eval_str("quote("+listtoquote_str(L)+")");
}
def dep6(Q) {
  if (type(Q) == 4) {
    Q = to_quote(Q);
  }  
  if (qt_is_dependent(Q,x)) return 1;
  if (qt_is_dependent(Q,u)) return 1;
  /* とりあえず 10 次までの f. --> なんとかせよ. */
  for (I=0; I<10; I++) {  
    if (qt_is_dependent(Q,idxtov(u,I))) return 1;
  }
  return 0;
}
def diff_lin(F,G) {
  if (type(F) == 4) F=to_quote(F);
  if (type(G) == 4) G=to_quote(G);
  return qt_replace(quote(diff(f)+diff(g)),[[f,F],[g,G]]);
}
def diff_mul(F,G) {
  F1 = dep6(F); G1 = dep6(G);
  if (type(F) == 4) F=to_quote(F);
  if (type(G) == 4) G=to_quote(G);
  if (F1 && G1)
    return qt_replace(quote(diff(f)*g+f*diff(g)),[[f,F],[g,G]]);
  if ((F1 == 1) &&  (G1 == 0)) 
    return qt_replace(quote(diff(f)*g),[[f,F],[g,G]]);
  if ((F1 == 0) &&  (G1 == 1)) 
    return qt_replace(quote(f*diff(g)),[[f,F],[g,G]]);
  if ((F1 == 0) && (G1 == 0))
    return qt_zero();
}
def qt_one() {
  return quote(1);
}
def diff_x_n(N) {
  N = eval_quote(N);
  N1=N-1;
  if (N1 == 0)  return qt_one();
  if (N1 == 1)  return quote(2*x);
  if (N1 > 1) return eval_str("quote("+rtostr(N)+"*x^"+rtostr(N1)+")");
}
/* F が u とか u_0, u_1, ... なら 1 を戻す. */
/* debug 用の入力.
  tr_check_pn(quote(u_1),quote(pn(x,is_u_variable(x))));
*/
def is_u_variable(F) {
  /* 述語の前の check point も debugger に欲しい. */
  /* print("is_u_variable: ",0); print(print_input_form(F)); */
  if (type(F) == 17) F=quotetolist(F);
  if (rtostr(F[0]) == "internal") {
    V = eval_str(rtostr(F[1]));
    if (vtoidx(V)[0] == "u") return 1;
  }
  return 0;
}
/*  u_i^n の微分をする.  n*u_{i+1}*u_i^{n-1} 
   Todo: もっと簡潔に quote を書けないか?
*/
def diff_u_n(F,N) {
  F = eval_quote(F);
  I = vtoidx(F); 
  if (length(I) == 1) I = 0; else I=I[1];
  N = eval_quote(N);
  N1=N-1;
  NextU = "u_"+rtostr(I+1);
  if (I == 0) U = "u"; else U = "u_"+rtostr(I);

  NN = objtoquote(N);
  NN1 = objtoquote(N1);
  NextU = objtoquote(eval_str(NextU));
  U = objtoquote(eval_str(U));

  if (N1 == 0)  return NextU;
  if (N1 == 1)  return qt_replace(quote(2*up*uu),[[up,NextU],[uu,U]]);
  if (N1 > 1) return qt_replace(quote(n*up*uu^m),[[up,NextU],[uu,U],
     [n,NN],[m,NN1]]);
}

def test6b() {
  T1=[quote(diff(x)),[qt_one]];
  T2=[quote(diff(x^pn(n))),[diff_x_n,n]];  /* is_poly? が欲しい. */
  R1=[quote(diff(pn(f)+pn(g))),[diff_lin,f,g]];
  R2=[quote(diff(pn(f)*pn(g))),[diff_mul,f,g]];

  A = quote(diff(2*4*x^3+x));
  print(print_input_form(A));
  R=tr_apply_or_rules(A,[R1,R2,T1,T2]);
  return R;
}

/* Use Debug2=1; は debug にとても有益. */
def test6() {
  T1=[quote(diff(x)),[qt_one]];
  T2=[quote(diff(x^pn(n))),[diff_x_n,n]];  /* is_poly? が欲しい. */
  T3=[quote(diff(pn(f,is_u_variable(f))^pn(n))),[diff_u_n,f,n]]; 
  R1=[quote(diff(pn(f)+pn(g))),[diff_lin,f,g]];
  R2=[quote(diff(pn(f)*pn(g))),[diff_mul,f,g]];

  /* A = quote(diff(2*x^3+x));*/
  A = quote(diff(2*u^3+x)); 
  print(print_input_form(A));
  R=tr_apply_or_rules(A,[R1,R2,T1,T2,T3]);
  return R;
}
end$