/* $OpenXM: OpenXM/src/asir-contrib/testing/rewriting.rr,v 1.1 2005/03/30 05:10:40 takayama Exp $ */

/*  
  OpenXM版の Risa/Asir で実行のこと. OpenXM 版の関数を用いるため.
*/
/* $Id: quote2.rr,v 1.5 2005/03/30 05:03:44 taka Exp $
  このファイルは quotetolist でリストに変換したデータに対して
  パターンマッチおよびそれを応用した変形を行う.
  テストプログラムのため効率は無視.   (append の多用, 無駄な2重呼び出し, など))
*/

extern Debug$
Debug=0$
def dprint(X) {
  if (Debug) print(X);
}
def dprint0(X) {
  if (Debug) print(X,0);
}

/*
  変数パターンの書き方
  pn(name)  pattern の前と後ろをとり pn
  pn("x")

  Todo: pn(name,length,type) 
        pn("x","rest")?

  関数パターンの書き方. (Todo)
  fn(name,argv)
  fn("f",1,23) --> f(1,23) へ.
*/

/*
  Rule の例1:
    sin(3*@pi) 等を 0 に書き換える例:
      quote(sin(pn("n")*@pi))
  --> f(n)
  
  def f(X) { if (Xが整数) return 0; else sin(X*@pi); }

   Rule の左辺は quote 型のパターン. 右辺はかならず asir の関数.

   [function,sin,[b_op,*,[function,pn,[internal,x]],[function,@pi]]]
   [function,fn,[internal,f],[internal,x]]

   下の test0(), test1(), test2() を参照.  
*/

/*
 例: 不定積分.  test3() を参照.
*/

/*   Todo:
 例. Mathematica の N[ ] 相当の関数をユーザが書けるように.
    nn(sin(cos(@pi)+sqrt(2))) 
    --> nn(sin(nn(cos(nn(@pi)))+nn(sqrt(nn(2)))))

 例: 冪級数の計算を quote で実現.
        sort や expand は組み込みで.

 例: Mathematica の Expand[], Toghether[] 相当のもの.

 例: D の掛け算を パターンマッチで実現.  

 例: (x^(1/n))^n --> x 等.

*/

/*
 トップレベルの関数達.  (stylesheet の考えに似てる.)
  apply_rule1(Obj,rule). 
  apply_rule1 は iterator の一種. 右辺はつねに関数.

  Todo: rules はユーザ定義のものと default rule 名がある.
      たとえば sort とか展開, 0 の削除は組み込み rule として欲しい.
*/

def node(F) {
   return [F[0],F[1]];
}
/* Number of  child */
def nchild(F) {
   return length(F)-2;
}
def child(F,K) {
   return F[K+2];
}

/* 
   リスト F が リスト P に(先頭からの比較で)マッチしたら 1.
   そうでないから 0.  幅優先探索.
    Todo: P に任意関数を含む仕組みはまだ実装してない.  
         たとえば  quote(nn(fn("f"))) 
         この場合 quote(nn(sin(1.3))) に f=sin , 引数 1.3 で match.
         この場合 quote(nn(cos(1.3))) に f=cos , 引数 1.3 で match.
         nn(f(g(x)+h(x))) --> nn(f(nn(g(x))+nn(h(x)))) としたい.

*/
def match0(F,P)  {
  dprint0("F="); dprint(F);
  dprint0("P="); dprint(P);  

  if (type(F) != type(P)) return 0;
  if (type(F) != 4) {
    if (F == P) return 1;
    else return 0;
  }
  Node = node(F);
  Node2 = node(P);
  if (Node2 == ["function","pn"]) return 2;
  if (Node != Node2) return 0;
  N = nchild(F);
  if (N != nchild(P)) return 0;
  for (I=0; I<N; I++) {
     C = child(F,I);
     C2 = child(P,I);
     if (!match0(C,C2)) return 0;
  }
  return 1;
}

/* F と P が match0 するとき  bindingTable をもどす.
  [[変数の名前(文字列), 値(list)], ...]
*/
def makeBind(F,P) {
  Ans = [ ];
  if (F == P) return Ans;

  Node = node(F);
  Node2 = node(P);

  if (Node2 == ["function", "pn"]) {
     Ans = append(Ans,[[P[2][1],F]]);
     return Ans;
  }
  N = nchild(F);
  for (I=0; I<N; I++) {
     C = child(F,I);
     C2 = child(P,I);
     Ans = append(Ans,makeBind(C,C2));
  }
  return Ans;
}

/* 
   Tree の中を幅優先探索で検索して 置き換える.
   幅優先探索なので, 同じ rule にマッチするものが入れ子になった場合, 
   内側は置き換えられない.
   Todo: 深さ優先探索.
   Todo: 書き換えがおこったかのフラグ.
*/
def rp(F,P,Q) {
  dprint0("rp, F="); dprint(F);
  dprint0("rp, P="); dprint(P);
  dprint0("rp, Q="); dprint(P);
  if (match0(F,P)) {
     BindTable = makeBind(F,P);
     dprint0("BindTable="); dprint(BindTable);
     return applyfunction0(Q,BindTable);
  }
  if (type(F) != 4) return F;
  Node = node(F);
  N = nchild(F);
  Ans = Node;
  for (I=0; I<N; I++) {
    T = rp(child(F,I),P,Q);
    Ans = append(Ans,[T]); 
  }
  return Ans;
}

/* ["f","x"],[["x",[internal,3]]]  の時は
   f(3) を計算する.
*/
def applyfunction0(Q,BindTable) {
  B = [ ];
  N = length(BindTable);
  /* BindTable の右辺値を quote(...) なる文字列に */
  for (I=0; I<N; I++) {
    B = append(B,[[BindTable[I][0],"quote("+quote_input_form_quote_list(BindTable[I][1])+")"]]);
  } 
  dprint0("applyfunction0: "); dprint(B);
  N = length(Q)-1; /* 引数の数 */
  M = length(B);   /*  binding table のサイズ */
  R = Q[0]+"(";
  for (I=0; I<N; I++) {
    X = rtostr(Q[I+1]); /* 変数 */
    /* binding Table をサーチ */
    for (J=0; J<M; J++) {
       Y = rtostr(B[J][0]);
       if (X == Y) {
          R = R+B[J][1];
          if (I != N-1) R = R+",";
          break;
       }
       if (J == M-1) error("No binding data.");
    }
  }
  R = R+")";
  dprint0("R="); dprint(R);
  return eval_str(R);
}

/* 右規則関数.  sin(整数*@pi) を 0 に */
def r_sin_int(X) {
  /* いま X は quote 型 */
  Y = quotetolist(X);
  /* Todo: このようなものを作る機能は組み込みで欲しい. */
  R = "quote(sin("+quote_input_form_quote_list(Y)+"*@pi))";
  print(R);
  R = eval_str(R);
  /* Todo: X が 数字かどうか調べる機能も組み込みで欲しい.
  */
  if (Y[0] == "internal") {
     Z = eval_str(rtostr(Y[1]));
  }else{ 
    return quotetolist(R);
  }
  if (type(Z) == 0) return quotetolist(quote(0));
  if ((type(Z) == 1) &&   (ntype(Z) == 0)) return quotetolist(quote(0));
  return quotetolist(R);
}

/* L が左規則. R が右規則.  幅優先探索.
  例: 
    apply_rule1(quote(1+sin(3*@pi)*sin(@pi/2)),
                quote(sin(pn("x")*@pi)),
                ["r_sin_int","x"]);
*/
def apply_rule1(Obj,L,R) {
  dprint("--------  start of apply_rule1 ------------ ");
  Obj = quotetolist(Obj);
  L = quotetolist(L);
  R = rp(Obj,L,R);
  RR = "quote("+quote_input_form_quote_list(R)+")";
  dprint("--------  end of apply_rule1 ------------ ");
  return eval_str(RR);  
}

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

/* 右規則関数.   0 を戻す. */
def r_zero() {
  return quotetolist(quote(0));
}

/* 右規則関数.   恒等式 */
def r_id(X) {
  return quotetolist(X);
}

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

/* 幅優先探索の場合, これは simplify できず. */
def test2() {
  Rule1=[quote(sin(pn("x")*@pi)),["r_sin_int","x"]];
  R0 = quote(1+sin(sin(2*@pi)*@pi)*sin(@pi/2));
  print(print_input_form(R0));
  R=apply_rule1(R0,Rule1[0],Rule1[1]);
  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=apply_rule1(R,Rules[I][0],Rules[I][1]);
    }
  }
  return R;
}

end$