/* $OpenXM: OpenXM/src/asir-contrib/testing/taka_series.rr,v 1.1 2005/05/04 05:47:03 takayama Exp $ */

module tkseries;

localf quote_input_form0$
localf diff1$
localf diffn$
localf expand1$
localf myeval$

/* Todo */
/* qt.dtoq, qtod;  quote<-->dp   */
/* qt.ltoq, qtol;  quote<-->list */

#define O_N 1
#define O_P 2
#define O_R 3
#define O_LIST 4
#define O_VECT 5
#define O_MAT 6
#define O_STR 7
#define O_COMP 8
#define O_DP 9
#define O_USINT 10
#define O_ERR 11
#define O_GF2MAT 12
#define O_MATHCAP 13
#define O_F 14
#define O_GFMMAT 15
#define O_BYTEARRAY 16
#define O_QUOTE 17
#define O_OPTLIST 18
#define O_SYMBOL 19
#define O_RANGE 20
#define O_TB 21
#define O_DPV 22
#define O_QUOTEARG 23
#define O_VOID -1

def quote_input_form0(Q) {
  Q = quotetolist(Q);
  return quote_input_form_quote_list(Q);
}
def diffn(F,X,N) {
  if ((type(F) == O_LIST) || (type(F) == O_VECT) || (type(F) == O_MAT)) {
    return map(diffn,F,X,N); /* It will not work for O_MAT */
  }
  /* We assume that F is O_QUOTE, O_P or O_R */
  if (N == 0) return F;
   for (I=0; I<N; I++) {
     F = diff1(F,X);
  }
  return F;
}

def diff1(F,X) {
  if (type(F) == O_QUOTE) {
     /* It has not yet been implemented. */
     R = "quote(calc_diff("+quote_input_form0(F)+","+rtostr(X)+"))";
     return eval_str(R);
  }
  /* We assume O_P or O_R */
  F=diff(F,x);
  return F;
}

def myeval(C) {
  C = objtoquote(C);
  C = tr.simp_sin(C);  /* It is defined in test1-tr.rr */

  /* Simplifier for quote should be called instead. */
  if (qt.is_rational(C)) C = eval_quote(C);
  return objtoquote(C);  
}

/* Series expansion of F at X=C to the degree N */
def expand1(F,X,C,N) {
  R = 0;
  Fn = F;
  Cn = base_replace(Fn,[[X,C]]);
  R = tr.simp_zero(R+myeval(Cn));
  for (I=1; I<=N; I++) {
     Fn = diff1(Fn,X);
     Cn = base_replace(Fn,[[X,C]]);
     Cn = myeval(Cn)/objtoquote(fac(I));
     Cn = qt.cancel_number(Cn); /* simplify, e.g.,  2/4 --> 1/4 */
     if (qt.is_minus(Cn)) Cn = qt.add_paren0(Cn); /* -1/2 --> (-1/2) */
     Po = objtoquote(X-C)^objtoquote(I);
     R = R + tr.simp_zero(Cn*Po);
     R = tr.simp_zero(R);
  }
  Rem = eval_str("quote(oO("+rtostr(X)+","+rtostr(N+1)+"))");
  return R+Rem;
}


endmodule;

/* Screenshot 
[1798] load("tr.rr"); load("taka_series.rr");
1

[1874] tkseries.expand1(cos(x),x,0,7);
quote(1+-1/2*x^2+1/24*x^4+-1/720*x^6+oO(x, 8))

By tr.simp_unary(R),  +- changes to - .

[1875] print_tex_form(@@);
     1 -  \frac{ 1} { 2}   {x}^{ 2} +  \frac{ 1} { 24}   {x}^{ 4}  -  \frac{ 1} { 720}   {x}^{ 6} + {oO}( {x},  8)

[1876] tkseries.expand1(sin(x),x,0,7);
quote(1/1*x^1+-1/6*x^3+1/120*x^5+-1/5040*x^7+oO(x, 8))

[1877] print_tex_form(@@);
      \frac{ 1} { 1}   {x}^{ 1}  -  \frac{ 1} { 6}   {x}^{ 3} +  \frac{ 1} { 120}   {x}^{ 5}  -  \frac{ 1} { 5040}   {x}^{ 7} + {oO}( {x},  8)

*/

end$