/*
   Create functions for computation of derivatives
   of RHS and related functions.
   N.B. Content creates derivatives of orders 1-3 by itself;
        For orders 4-5 the External Differentiator (ExtDiff, currently Maple) is used.
*/

#include "common.h"

#if ED
#define MAXDER		5	/* ==DERMAX in function.h */
#else	/* ED */
#define MAXDER		3	/* max order of derivatives computed by Content */
#endif	/* ED */

#define PN		(4+MAXDER*4)

#define CD_OOPEN	1	/* cannot open output file */

#include "corefunc.h"

static FILE *out;

#if ED

static int tfn;			/* number of temp files with ExtDiff input */

#define MIF	"tmp%i.mi"	/* name of ExtDiff input file */
#define MOF	"tmp%i.mo"	/* name of ExtDiff output file */
static char mif[sizeof(MIF)+10],mof[sizeof(MOF)+10];
static FILE *edif;

#endif	/* ED */

/*==================*/
/* Common functions */

/* Map 2nd- and 3d-order derivatives to intervals [0...] */
/* See autodif.c */
#define ind2_(i,j) ((i)+(((j)*((j)+1))>>1))
#define ind3_(i,j,k) ((i)+(((j)*((j)+1))>>1)+(k)*((k)+1)*((k)+2)/6)
#if ED
#define ind4_(i,j,k,l) (ind3_(i,j,k)+(l)*((l)+1)*((l)+2)*((l)+3)/24)
#define ind5_(i,j,k,l,m) (ind4_(i,j,k,l)+(m)*((m)+1)*((m)+2)*((m)+3)*((m)+4)/120)
#endif	/* ED */

typedef void (PNTR ProcessName)(CharPtr name, int len, int dim, VoidPtr user);

#define DIM_F	"%n [%n %i ]%n"		/* format used for dimension parsing */
#define ID_F	"%*[a-zA-Z0-9_()]"DIM_F	/* id[dim] */
#define GRO_F	"{%*[^}]}"DIM_F		/* {id id ...}[dim] */

/* Enumerates all the names in a given list (the same as in function.c) */
Local(Int2) EnumNames(CharPtr plist, ProcessName process, VoidPtr user) {
  int pos_id_e,pos_num,pos_e,dim;
  Int2 rc=0;
  for (; *plist && !rc; ) {
    pos_e=-1;
    sscanf(plist," %n",&pos_e);
    plist+=pos_e;
    if (pos_e<0 || !*plist) break;		/* only trailing blanks. !*p for cc */
    pos_id_e=pos_num=pos_e=0;
    dim=1;
    if (IS_ALPHA(*plist) || *plist=='_' || *plist==*GRO_F) {
      sscanf(plist, *plist==*GRO_F ? GRO_F : ID_F ,&pos_id_e,&pos_num,&dim,&pos_e);
      if (pos_id_e==0) pos_e=pos_id_e=(int)StrLen(plist);	/* last/only name is right-justified */
      else
	if (pos_num==0) pos_e=pos_id_e;			/* no [ after name */
	else
	  if (pos_e==0 || dim<2) rc=1;
    } else rc=2;
    process(plist,pos_id_e,dim,user);
    plist+=pos_e;
    pos_e=-1;
    sscanf(plist," ,%n",&pos_e);
    if (pos_e>=0) plist+=pos_e;
  }
  return rc;
}

/* substitute definitions of local vartibles by TempVars */
Local(CharPtr) SubstituteLocals(CharPtr RhsPtr) {
  int dim,pos;
  Char c;
  Char b[200];
  while (sscanf(RhsPtr," %s%n",b,&pos)==1) {
    /* Only directives and int are allowed and they are skipped */
    if (*b=='#') {
      RhsPtr+=pos;
      pos=0;
      sscanf(RhsPtr,"%[^\n]\n%n",b+StrLen(b),&pos);
      fprintf(out,"%s\n",b);
      RhsPtr+=pos;
      continue;
    }
    if (StrCmp(b,"int")==0) {
      RhsPtr+=pos;
      sscanf(RhsPtr,"%[^;];%n",b+StrLen(b),&pos);
      fprintf(out,"%s;\n",b);
      RhsPtr+=pos;
      continue;
    }
    if (StrCmp(b,"double")) break;
    RhsPtr+=pos;
    do {
      sscanf(RhsPtr," %[a-zA-Z0-9] %n",b,&pos);
      RhsPtr+=pos;
      if (*RhsPtr=='[') {
	sscanf(RhsPtr,"[ %i ]%n",&dim,&pos);
	RhsPtr+=pos;
      } else dim=1;
      
      if (sscanf(RhsPtr," %c%n",&c,&pos)==1) RhsPtr+=pos;
      else c=';';
      fprintf(out,"  TempVar(%s,%i);\n",b,dim);
    } while (c==',');
  }
  return RhsPtr;
}


/* Substitutes all occurences of *what by *with in *source */
Local(CharPtr) StrSubst(CharPtr source, CharPtr what, CharPtr with) {
  CharPtr p;
  CharPtr n;
  for (; (p=StringStr(source,what),p); ) {
    n=(CharPtr)MemNew(StrLen(source)+1-StrLen(what)+StrLen(with));
    StringNCpy(n,source,(Int2)(p-source));
    n[(Int2)(p-source)]='\0';
    StringCat(n,with);
    StringCat(n,p+StrLen(what));
    MemFree(source);
    source=n;
  }
  return source;
}

typedef void (PNTR EnumArgsCallback)(CharPtr name, int last);

/* Enumerates all the name in argument list */
Local(void) EnumArgs(CharPtr arglist, EnumArgsCallback process) {
  CharPtr p;
  int pos;
  Char name[30],c;
  p=StrChr(arglist,'(');
  if (p) arglist=p+1;
  while (sscanf(arglist," %*[a-zA-Z0-9] %[a-zA-Z0-9] %c%n",name,&c,&pos)==2) {
    process(name,c==')');
    arglist+=pos;
  }
}


static CharPtr DepVarPtr;
static CharPtr IndepVarPtr;
static CharPtr FuncPtr;
static CharPtr RhsPtr;
static CharPtr OutFile;
static CharPtr FuncAttr;
Local(CharPtr) HeadPtr[MAXDER];
Local(CharPtr) DefsPtr[MAXDER];
Local(CharPtr) TailPtr[MAXDER];
Local(CharPtr) UndefsPtr[MAXDER];

/* See function.c. Must be the same */
#define DER_NUM		1	/* numerically */
#define DER_SYM		2	/* analitically */
#define DER_USER	3	/* user-supplied */


/*=====================*/
/* C++ header and tail */

#define TMPFUN		"tmp_fun.c"	/* C++ file for user-defined functions */
#define TMPFUN1		"tmp_fun1.h"	/* #define INDEX_ ... */

static Char Header[]=
  "\n"
  "#include <stdio.h>\n\n"
  "#include \"autodif.h\"\n\n"
  "FILE *MStream_;\n"
  "FILE *FStream_;\n\n"
  ;

static Char Tail[]=
  "\n"
  "int main(void) {\n"
  "#define ORDER__ %i\n"
  "  MStream_=fopen(\"%s\",\"wt\");\n"
  "  FStream_=fopen(\""TMPFUN"\",\"wt\");\n"
  "  fprintf(MStream_,\"#include \\\""TMPFUN1"\\\"\\n\");\n\n"
  "  fprintf(MStream_,\"static double v__[INDEX_];\\n\");\n"
  "#if ORDER__&1\n  fprintf(MStream_,\"static double d1__[INDEX_][D1N_];\\n\");\n#endif\n"
  "#if ORDER__&2\n  fprintf(MStream_,\"static double d2__[INDEX_][D2N_];\\n\");\n#endif\n"
  "#if ORDER__&4\n  fprintf(MStream_,\"static double d3__[INDEX_][D3N_];\\n\");\n#endif\n"
  "  fprintf(MStream_,\"double sv__[SFINDEX_], sd__[SFINDEX_];\\n\");\n"
  "  fprintf(MStream_,\"#include \\\""TMPFUN"\\\"\\n\");\n\n"
  "#if ORDER__&1\n  _Hod_=1; rtInit_(MStream_,FStream_);\n  %s\n  rtTerm_();\n#endif\n"
  "#if ORDER__&2\n  _Hod_=2; rtInit_(MStream_,FStream_);\n  %s\n  rtTerm_();\n#endif\n"
  "#if ORDER__&4\n  _Hod_=3; rtInit_(MStream_,FStream_);\n  %s\n  rtTerm_();\n#endif\n"
  "#if ORDER__&8\n  _Hod_=4; rtInit_(MStream_,FStream_);\n  %s\n  rtTerm_();\n#endif\n"
  "#if ORDER__&16\n  _Hod_=5; rtInit_(MStream_,FStream_);\n  %s\n  rtTerm_();\n#endif\n"
  "  fclose(MStream_);\n"
  "  fclose(FStream_);\n"
  "  FStream_=fopen(\""TMPFUN1"\",\"wt\");\n"
  "  fprintf(FStream_,\"#define INDEX_   %%i\\n\",rtGetMaxuf_());\n"
  "  fprintf(FStream_,\"#define SFINDEX_ %%i\\n\",rtGetMaxsf_());\n"
  "#if ORDER__&1\n  fprintf(FStream_,\"#define D1N_     %%i\\n\",_Niv_);\n#endif\n"
  "#if ORDER__&2\n  fprintf(FStream_,\"#define D2N_     %%i\\n\",ind2_(_Niv_-1,_Niv_-1)+1);\n#endif\n"
  "#if ORDER__&4\n  fprintf(FStream_,\"#define D3N_     %%i\\n\",ind3_(_Niv_-1,_Niv_-1,_Niv_-1)+1);\n#endif\n"
  "  fclose(FStream_);\n"
  "  return 0;\n"
  "}\n"
  ;

/*================================================*/
/* Create source C++ text for partial derivatives */

static Char DerHead[]=
  "\n"
  "void %.*s_Der%i(void) {\n"
  "  _Niv_=%i;\n"
  "  _Hod_=%i;\n"
  ;

static Char DerTail[]=
  "\n"
  "}\n\n"
  ;

static int ivn,dvn;	/* independent and dependent variable numbers */

Local(void) JIGVar(CharPtr name, int len, int dim, VoidPtr vn, Boolean indep) {
  int l;
  CharPtr p;
  Char c;
  c=name[len];
  name[len]='\0';
  p=StrChr(name,'(');
  name[len]=c;
  if (p) {
    l=(int)(p-name);
    fprintf(out,"#ifdef %.*s\n",l,name);
    fprintf(out,"  Gradient %.*s_%.*s(%.*s,%i,%i,%s);\n",
		l,name,(int)(name+len-p-2),p+1,len,name,*(int *)vn,dim,indep?"INDEP":"GRAD");
    fprintf(out,"#else\n");
    fprintf(out,"  Gradient %.*s_%.*s(\"%.*s\",%i,%i,%s);\n",
		l,name,(int)(name+len-p-2),p+1,len,name,*(int *)vn,dim,indep?"INDEP":"GRAD");
    fprintf(out,"#endif\n");
    fprintf(out,"#undef %.*s\n",l,name);
    fprintf(out,"#define %.*s(n) %.*s_##n\n",l,name,l,name);
  } else
    fprintf(out,"  %s(%.*s,%i,%i);\n",indep?"IndepVar":"GradVar",len,name,*(int *)vn,dim);
  *(int *)vn+=dim;
}

Local(void) JIndepVar(CharPtr name, int len, int dim, VoidPtr ivn) {
  JIGVar(name,len,dim,ivn,TRUE);
}

Local(void) JGradVar(CharPtr name, int len, int dim, VoidPtr dvn) {
  JIGVar(name,len,dim,dvn,FALSE);
}

#if _WIN
#pragma argsused
#endif
Local(void) JIndepGradUndef(CharPtr name, int len, int dim, VoidPtr dvn) {
  CharPtr p;
  Char c;
  c=name[len];
  name[len]='\0';
  p=StrChr(name,'(');
  name[len]=c;
  if (p) {
    fprintf(out,"#undef %.*s\n",(int)(p-name),name);
  }
}

#if _WIN
#pragma argsused
#endif
static void JVarCount(CharPtr name, int len, int dim, VoidPtr ivn) {
  *(int PNTR)ivn+=dim;
}

/* 
   Lines containing "/*==" are printed as literals,
   e.g. without being put into fprintf(MStream,...);
*/
Local(void) PrintFprintf(CharPtr p) {
  CharPtr q;
  Char c;
  while (*p) {
    q=StrPBrk(p,"\n\r");
    if (q) {c=*q; *q='\0';}
    else q=p+StrLen(p);
    if (StrStr(p,"/*=="))
      fprintf(out,"%.*s\n",(int)(q-p),p);	/*== comment */
    else
      fprintf(out,"  fprintf(MStream_,\"%.*s\\n\");\n",(int)(q-p),p);
    if (q) *q=c;
    while (*q=='\n' || *q=='\r') q++;
    p=q;
  }
}

static void Der(CharPtr name, int len, int lvl) {
  CharPtr p;
  CharPtr rhs;
  /* Compute the number of independent vars */
  ivn=0;
  EnumNames(IndepVarPtr,JVarCount,&ivn);
  fprintf(out,DerHead,len,name,lvl,ivn,lvl);
  /* Output the header */
  PrintFprintf(HeadPtr[lvl-1]);
  /* Output fprintf(...,#defines) that map user's names onto formal parameter */
  PrintFprintf(DefsPtr[lvl-1]);
  /* Output declarations of independent and dependent variables and
     run-time #defines that map user's names onto _x, _p, _d1 */
  ivn=1;
  EnumNames(IndepVarPtr,JIndepVar,&ivn);
  dvn=1;
  EnumNames(DepVarPtr,JGradVar,&dvn);
  /* substitute depvar' by depvar_ */
  for (p=RhsPtr; (p=StrChr(p,'\''))!=NULL; *p++='_');	/* substitute ' by _ */
  /* printf TempVars instead of locals */
  rhs=SubstituteLocals(RhsPtr);
#if ED
  if (lvl>3) {	/* specific for External differentiator */
    if (tfn==0) {
      p=GetParamString(SFS_COMPILER,"EXTDIFFB");
      edif=fopen(p,"wt");
      if (!edif) fprintf(out,"#error Cannot create file '%s'\n",p);
    }
    if (edif) {
      tfn++;
      sprintf(mif,MIF,tfn);
      sprintf(mof,MOF,tfn);
      /* Put to External Differentiator's batch file */
      p=GetParamString(SFS_COMPILER,"EXTDIFF");
      fprintf(edif,p,mif,mof);
      fprintf(edif,"\n");
      /* Pre-source text hook */
      fprintf(out,"  fprintf(MStream_,\"#include \\\"%s\\\"\\n\");\n",mof);
      fprintf(out,"  rtExtDiff_(\"%s\",1);\n",mif);
    }
  }
#endif	/* ED */
  /* output rhs */
  fprintf(out,rhs);	/* Body */
  fprintf(out,"\n");
#if ED
  if (lvl>3) {	/* specific for External differentiator */
    if (edif) {
      /* Post-source text hook */
      fprintf(out,"  rtExtDiff_(\"%s\",0);\n",mof);
    }
  }
#endif	/* ED */
  EnumNames(IndepVarPtr,JIndepGradUndef,NULL);
  EnumNames(DepVarPtr,JIndepGradUndef,NULL);
  /* Output the tail and #undefs */
  PrintFprintf(TailPtr[lvl-1]);
  PrintFprintf(UndefsPtr[lvl-1]);
  fprintf(out,DerTail,lvl);
}

/*==========================================*/
/* Modify source for user-defined functions */

Local(int) FuncArgN;

Local(CharPtr) head;

Local(void) CppArgs(CharPtr name, int last) {
  fprintf(out,"Gradient _%s%s",name,last ? "" : ", ");
  FuncArgN++;
}

Local(void) CppFormat(CharPtr name, int last) {
  fprintf(out,"  fprintf(MStream_,\"%%s%s\",GetExpr_(_%s));\n",last ? "" : ",",name);
}
#if _WIN
#pragma argsused
#endif
Local(void) CppIndep(CharPtr name, int last) {
  fprintf(out,"    IndepVar(%s,%i,1);\n",name,++FuncArgN);
}

Local(int) CppChain_i,CppChain_j,CppChain_k;
Local(CharPtr) CppChain_name;

#if _WIN
#pragma argsused
#endif
Local(void) CppChain1(CharPtr name, int last) {
  fprintf(out,"    _p_=Add__(_p_,Mult__(\"d1__%s[%i]\",GetDer1T_(_%s,_i_),0),'b');\n",
	  CppChain_name,CppChain_i++,name);	
}

#if _WIN
#pragma argsused
#endif
Local(void) CppChain2a(CharPtr name, int last) {
  fprintf(out,"      _q_=Add__(_q_,Mult__(\"d2__%s[%i]\",GetDer1T_(_%s,_j_),0),'b');\n",
	  CppChain_name,ind2_(CppChain_i,CppChain_j),name);
  CppChain_j++;	
}

#if _WIN
#pragma argsused
#endif
Local(void) CppChain2(CharPtr name, int last) {
  fprintf(out,"      _p_=Add__(_p_,Mult__(\"d1__%s[%i]\",GetDer2T_(_%s,_i_,_j_),0),'b');\n",
	  CppChain_name,CppChain_i,name);
  fprintf(out,"      _q_=0;\n");
  CppChain_j=0;
  EnumArgs(head,CppChain2a);
  fprintf(out,"      _p_=Add__(_p_,Mult__(_q_,GetDer1T_(_%s,_i_),0),'b');\n",
	  name);
  CppChain_i++;
}

Local(CharPtr) CppChain3a_name,CppChain3b_name;

#if _WIN
#pragma argsused
#endif
Local(void) CppChain3b(CharPtr name, int last) {
  fprintf(out,"      _w_=Mult__(Mult__(GetDer1T_(_%s,_i_),GetDer1T_(_%s,_j_),0),GetDer1T_(_%s,_k_),'l');\n",
		CppChain3a_name,CppChain3b_name,name);
  if (CppChain_i!=CppChain_k) {
    if (CppChain_i!=CppChain_j && CppChain_j!=CppChain_k) {
      fprintf(out,"      _w_=Add__(_w_,Mult__(Mult__(GetDer1T_(_%s,_i_),GetDer1T_(_%s,_k_),0),GetDer1T_(_%s,_j_),'l'),'b');\n",
		CppChain3a_name,CppChain3b_name,name);
      fprintf(out,"      _w_=Add__(_w_,Mult__(Mult__(GetDer1T_(_%s,_j_),GetDer1T_(_%s,_i_),0),GetDer1T_(_%s,_k_),'l'),'b');\n",
		CppChain3a_name,CppChain3b_name,name);
    }
    fprintf(out,"      _w_=Add__(_w_,Mult__(Mult__(GetDer1T_(_%s,_j_),GetDer1T_(_%s,_k_),0),GetDer1T_(_%s,_i_),'l'),'b');\n",
		CppChain3a_name,CppChain3b_name,name);
    fprintf(out,"      _w_=Add__(_w_,Mult__(Mult__(GetDer1T_(_%s,_k_),GetDer1T_(_%s,_i_),0),GetDer1T_(_%s,_j_),'l'),'b');\n",
		CppChain3a_name,CppChain3b_name,name);
    if (CppChain_i!=CppChain_j && CppChain_j!=CppChain_k)
      fprintf(out,"      _w_=Add__(_w_,Mult__(Mult__(GetDer1T_(_%s,_k_),GetDer1T_(_%s,_j_),0),GetDer1T_(_%s,_i_),'l'),'b');\n",
		CppChain3a_name,CppChain3b_name,name);
  }
  fprintf(out,"      _r_=Add__(_r_,Mult__(\"d3__%s[%i]\",_w_,0),'b');\n",
	      CppChain_name,ind3_(CppChain_i,CppChain_j,CppChain_k));
  CppChain_k++;
}

#if _WIN
#pragma argsused
#endif
Local(void) CppChain3a(CharPtr name, int last) {
  fprintf(out,"      _w_=Mult__(GetDer1T_(_%s,_i_),GetDer2T_(_%s,_j_,_k_),0);\n",
		CppChain3a_name,name);
  if (CppChain_i!=CppChain_j)
    fprintf(out,"      _w_=Add__(_w_,Mult__(GetDer1T_(_%s,_i_),GetDer2T_(_%s,_j_,_k_),0),'b');\n",
		name,CppChain3a_name);
  fprintf(out,"      _w_=Add__(_w_,Mult__(GetDer1T_(_%s,_j_),GetDer2T_(_%s,_i_,_k_),0),'b');\n",
		CppChain3a_name,name);
  if (CppChain_i!=CppChain_j)
    fprintf(out,"      _w_=Add__(_w_,Mult__(GetDer1T_(_%s,_j_),GetDer2T_(_%s,_i_,_k_),0),'b');\n",
		name,CppChain3a_name);
  fprintf(out,"      _w_=Add__(_w_,Mult__(GetDer1T_(_%s,_k_),GetDer2T_(_%s,_i_,_j_),0),'b');\n",
		CppChain3a_name,name);
  if (CppChain_i!=CppChain_j)
    fprintf(out,"      _w_=Add__(_w_,Mult__(GetDer1T_(_%s,_k_),GetDer2T_(_%s,_i_,_j_),0),'b');\n",
		name,CppChain3a_name);
  fprintf(out,"      _q_=Add__(_q_,Mult__(\"d2__%s[%i]\",_w_,0),'b');\n",
	      CppChain_name,ind2_(CppChain_i,CppChain_j));
  CppChain_k=0;
  CppChain3b_name=name;
  EnumArgs(head,CppChain3b);
  CppChain_j++;
}

#if _WIN
#pragma argsused
#endif
Local(void) CppChain3(CharPtr name, int last) {
  fprintf(out,"      _p_=Add__(_p_,Mult__(\"d1__%s[%i]\",GetDer3T_(_%s,_i_,_j_,_k_),0),'b');\n",
	  CppChain_name,CppChain_i,name);
  CppChain_j=0;
  CppChain3a_name=name;
  EnumArgs(head,CppChain3a);
  CppChain_i++;
}

#if _WIN
#pragma argsused
#endif
Local(void) CppCall(CharPtr name, int last) {
  fprintf(out,"    fprintf(MStream_,\"%s%s\");\n",name,last ? ")" : ",");
}

static ByteStorePtr undefs;

Local(CharPtr) prefix[MAXDER];

static void UserDefFuncs(int lvl) {
  CharPtr p;
  CharPtr q;
  CharPtr w;
  CharPtr rname;
  int pos,i,j,k,num,lvl1=lvl-1;
  Char name[30];
#if ED
  if (lvl>3 && FuncPtr && (StrTrim(FuncPtr),*FuncPtr)) {
    fprintf(out,"#error Derivatives of order %i cannot be build in the presence of local functions\n",lvl);
    return;
  }
#endif	/* ED */
  /* strip comments */
  for (p=FuncPtr; (p=StringStr(p,"//"))!=NULL; ) {
    q=StrChr(p,'\n');
    MemFill(p,' ',q ? (size_t)(q-p) : StrLen(p));
    p=q;
  }
  for (p=FuncPtr; (p=StringStr(p,"/*"))!=NULL; ) {
    q=StrStr(p,"*/");
    MemFill(p,' ',q ? (size_t)(q-p)+2 : StrLen(p));
    p=q;
  }
  /* process every function */
  StrCpy(name,prefix[lvl1]);	/* to make different order versions of functions different */
  rname=name+StrLen(prefix[lvl1]);	/* to real name without prefix */
  for (p=FuncPtr; sscanf(p," %s %n",rname,&pos)==1; ) {
    /* skip directives, if any */
    if (rname[0]=='#') {
      CharPtr defp=p;
      p+=pos;
      pos=-1;
      sscanf(p," %*[^\n\r] %n",&pos);
      if (pos<0) pos=StrLen(p);
      p+=pos;
      fprintf(out,"%.*s\n",(int)(p-defp),defp);
      continue;
    }
    /* find name and argument list */
    if (!StrCmp(rname,"static")) {
      p+=pos;
      sscanf(p," %s %n",rname,&pos);	/* read type */
    }
    p+=pos;
    sscanf(p," %[a-zA-Z0-9]",rname);	/* read name */
    /* don't try to differentiate user defined derivatives */
    if (sscanf(rname,"UserDer%i",&i)==1 && 1<=i && i<=3) {
      for (i=0; TRUE; ) {
	q=StrPBrk(p,"{}");
	if (q) {
	  p=q+1;
	  if (*q=='{') i++;
	  else {
	    i--;
	    if (i==0) break;
	  }
	} else {
	  p=StrChr(p,'\0');
	  break;
	}
      }
      continue;
    }
    fprintf(out,"#define %s %s\n",rname,name);
    BSWrite(undefs,"#undef ",7);
    BSWrite(undefs,rname,StrLen(rname));
    BSPutByte(undefs,'\n');
    q=StrChr(p,')')+1;
    head=(CharPtr)MemNew((size_t)(q-p)+1);
    StrNCpy(head,p,(size_t)(q-p));
    head[(size_t)(q-p)]='\0';
    p=StrChr(q,'{')+1;			/* p points to the next byte after opening { */
    /* int f?_name=1; */
    fprintf(out,"int f%i_%s=1;\n\n",lvl,rname);
    /* Gradient name(Gradient a1_,...) */
    fprintf(out,"static Gradient %s(",name);
    FuncArgN=0;
    EnumArgs(head,CppArgs);
    fprintf(out,") {\n  Gradient _t_;\n  char _b_[20],*_p_;\n  int _indx_,_i_;\n");
    switch (lvl) {
      case 3: fprintf(out,"  int _k_; char *_r_,*_w_;\n");
      case 2: fprintf(out,"  int _j_; char *_q_;\n");
    }
    fprintf(out,"  int s_Niv_=_Niv_;\n");
    /* if (first_name) {function's body} */
    fprintf(out,"  _Niv_=%i;\n",FuncArgN);	/* set _Niv_ before */
    fprintf(out,"  if (f%i_%s) {\n    Gradient _t_;\n    f%i_%s=0;\n",lvl,rname,lvl,rname);
    fprintf(out,"    rtFuncEntry_();\n");
    switch (lvl) {
      case 1: num=FuncArgN; break;
      case 2: num=ind2_(FuncArgN-1,FuncArgN-1)+1; break;
      case 3: num=ind3_(FuncArgN-1,FuncArgN-1,FuncArgN-1)+1; break;
    }
    fprintf(out,"    fprintf(MStream_,\"static double d%i__%s[%i];\\n\");\n",lvl,rname,num);
    fprintf(out,"    fprintf(MStream_,\"static double %s%s {\\n\");\n",prefix[lvl1],head);
    /* IndepVar for all arguments */
    FuncArgN=0;
    EnumArgs(head,CppIndep);
    /* TempVar for all temporaries */
    p=SubstituteLocals(p);
    /* function's body up to but not include return operator */
    for (q=p; (q=StrStr(q,"return"))!=NULL; q+=6)
      if (!IS_ALPHANUM(*(q-1)) && !IS_ALPHANUM(*(q+6))) break;
    if (!q) q=p+StrLen(p);
    fprintf(out,"%.*s\n\n",(int)(q-p),p);
    p=q;				/* p now points to the return operator */
    /* return */
    q=StrChr(p,'}')+1;			/* q now points to the next function */

    w=StrSubst(StrNCpy(w=(CharPtr)MemNew((int)(q-p)+1),p,(int)(q-p)),"\n","\\n");
    w=StrSubst(w,"\r","");
    fprintf(out,"    _t_=%.*s;\n",(int)(StrChr(w,';')-w-6),p+6);
    switch (lvl) {
      case 1:
	for (i=0; i<FuncArgN; i++)
	  fprintf(out,"    fprintf(MStream_,\"  d1__%s[%i]=%%s;\\n\",GetDer1_(_t_,%i));\n",
	  	rname,i,i);
	fprintf(out,"    fprintf(MStream_,\"  return %%s;\\n}\\n\\n\",GetExpr_(_t_));\n");
	break;
      case 2:
	for (i=0; i<FuncArgN; i++)
	  for (j=i; j<FuncArgN; j++)
	    fprintf(out,"    fprintf(MStream_,\"  d2__%s[%i]=%%s;\\n\",GetDer2_(_t_,%i,%i));\n",
	  	rname,ind2_(i,j),i,j);
       ret:
	fprintf(out,"    fprintf(MStream_,\"  return %s%s(\");\n",prefix[lvl1-1],rname);
	EnumArgs(head,CppCall);
	fprintf(out,"    fprintf(MStream_,\";\\n}\\n\\n\");\n");
	break;
      case 3:
 	for (i=0; i<FuncArgN; i++)
	  for (j=i; j<FuncArgN; j++)
	    for (k=j; k<FuncArgN; k++)
	      fprintf(out,"    fprintf(MStream_,\"  d3__%s[%i]=%%s;\\n\",GetDer3_(_t_,%i,%i,%i));\n",
		rname,ind3_(i,j,k),i,j,k);
       goto ret;
    }	
    MemFree(w);

    fprintf(out,"    rtFuncExit_();\n");
    fprintf(out,"  }\n");
    fprintf(out,"  _Niv_=s_Niv_;\n");

    /* fprintf(MStream,"v__[]=name(params)"); */
    fprintf(out,"  _indx_=rtGetIndex_();\n");
    fprintf(out,"  sprintf(_b_,\"v__[%%i]\",_indx_);\n");
    fprintf(out,"  fprintf(MStream_,\"  %%s=%s(\",_b_);\n",name);
    EnumArgs(head,CppFormat);
    fprintf(out,"  fprintf(MStream_,\");\\n\");\n");

    fprintf(out,"  _t_.SetExpr(_b_);\n");
    fprintf(out,"  for (_i_=0; _i_<_Niv_; _i_++) {\n");
    fprintf(out,"    _p_=0;\n");
    CppChain_name=rname;
    CppChain_i=0;
    EnumArgs(head,CppChain1);
    fprintf(out,"    if (_p_) {\n");
    fprintf(out,"      sprintf(_b_,\"d1__[%%i][%%i]\",_indx_,_i_);\n");
    fprintf(out,"      fprintf(MStream_,\"  %%s=%%s;\\n\",_b_,_p_);\n");
    fprintf(out,"      MemFree__(_p_);\n");
    fprintf(out,"      _t_.SetDer1(_i_,_b_);\n");
    fprintf(out,"    } else _t_.SetDer1(_i_,0);\n");
    if (lvl>=2) {
      fprintf(out,"    for (_j_=_i_; _j_<_Niv_; _j_++) {\n");
      fprintf(out,"      _p_=0;\n");
      CppChain_i=0;
      EnumArgs(head,CppChain2);
      fprintf(out,"      if (_p_) {\n");
      fprintf(out,"        sprintf(_b_,\"d2__[%%i][%%i]\",_indx_,ind2_(_i_,_j_));\n");
      fprintf(out,"        fprintf(MStream_,\"  %%s=%%s;\\n\",_b_,_p_);\n");
      fprintf(out,"        MemFree__(_p_);\n");
      fprintf(out,"        _t_.SetDer2(_i_,_j_,_b_);\n");
      fprintf(out,"      } else _t_.SetDer2(_i_,_j_,0);\n");
      if (lvl>=3) {
	fprintf(out,"    for (_k_=_j_; _k_<_Niv_; _k_++) {\n");
	fprintf(out,"      _p_=_q_=_r_=0;\n");
	CppChain_i=0;
	EnumArgs(head,CppChain3);
	fprintf(out,"      if (_p_||_q_||_r_) {\n");
	fprintf(out,"        if (!_q_) {_q_=_r_; _r_=0;}\n");
	fprintf(out,"        if (!_p_) {_p_=_q_; _q_=0;}\n");
	fprintf(out,"        sprintf(_b_,\"d3__[%%i][%%i]\",_indx_,ind3_(_i_,_j_,_k_));\n");
	fprintf(out,"        fprintf(MStream_,\"  %%s=%%s;\\n\",_b_,_p_);\n");
	fprintf(out,"        MemFree__(_p_);\n");
	fprintf(out,"        if (_q_) {\n");
	fprintf(out,"          fprintf(MStream_,\"  %%s+=%%s;\\n\",_b_,_q_);\n");
	fprintf(out,"          MemFree__(_q_);\n");
	fprintf(out,"        }\n");
	fprintf(out,"        if (_r_) {\n");
	fprintf(out,"          fprintf(MStream_,\"  %%s+=%%s;\\n\",_b_,_r_);\n");
	fprintf(out,"          MemFree__(_r_);\n");
	fprintf(out,"        }\n");
	fprintf(out,"        _t_.SetDer3(_i_,_j_,_k_,_b_);\n");
	fprintf(out,"      } else _t_.SetDer3(_i_,_j_,_k_,0);\n");
	fprintf(out,"    }\n");
      }      
      fprintf(out,"    }\n");
    }
    fprintf(out,"  }\n");

    fprintf(out,"  return _t_;\n");

    fprintf(out,"}\n");
    MemFree(head);
    p=q;
  }
#undef rname
}

/*==================*/
/* Main entry point */

typedef CharPtr PNTR Param;

Local(Int2) fn,level;

#if _WIN
#pragma argsused
#endif
Local(void) ProcessOneFunction(CharPtr name, int len, int dim, VoidPtr p) {
#define param ((Param)p)
  Int2 k;
  DepVarPtr=param[5+PN*fn];	/* list of dependent variable names */
  IndepVarPtr=param[6+PN*fn];	/* list of independent variable names */
  RhsPtr=param[8+PN*fn];	/* rhs */
  for (k=0; k<MAXDER; k++) {
    HeadPtr[k]=param[9+PN*fn+4*k];	/* header of der */
    DefsPtr[k]=param[10+PN*fn+4*k];	/* #defines for der */
    TailPtr[k]=param[11+PN*fn+4*k];	/* tail of der (return etc.) */
    UndefsPtr[k]=param[12+PN*fn+4*k];	/* #undefs for der */
  }
  Der(name,len,(int)level);
  fn++;
#undef param
}

#if ED
Local(CharPtr) l1,l2,l3,l4,l5;
#else	/* ED */
Local(CharPtr) l1,l2,l3;
#endif	/* ED */

#if _WIN
#pragma argsused
#endif
Local(void) AddOneName(CharPtr name, int len, int dim, VoidPtr dummy) {
  int l;
  l=StrLen(l1);
  sprintf(l1+l,"%.*s_Der1(); ",len,name);
  sprintf(l2+l,"%.*s_Der2(); ",len,name);
  sprintf(l3+l,"%.*s_Der3(); ",len,name);
#if ED
  sprintf(l4+l,"%.*s_Der4(); ",len,name);
  sprintf(l5+l,"%.*s_Der5(); ",len,name);
#endif	/* ED */
}

/*
   param contains the following information:
   0		name of output C file;
   1		name of work file for C++ text;
   2		functions attributes (type ATTRIBUTES name(params) {...};
   3		list of function names (functions to be differentiated);
   4		unused;
		Let n enumerates functions, n=0,1,..number_of_functions-1 (see 3).
#if ED
		Each groups of PN=4+5*4=24 pointers describes one function:
#else
		Each groups of PN=16 pointers describes one function:
#endif
   5+0+PN*n	list of dependent variables;
   5+1+PN*n	list of independent variables;
   5+2+PN*n	local functions (NULL for all except for the first function, i.e. Rhs);
   5+3+PN*n	function's body, e.g. Rhs;
#if ED
		The following 5 groups of 4 pointers related to
#else
		The following 3 groups of 4 pointers related to
#endif
		function for derivatives of order k+1 of the function (k=0,...);
   5+4+PN*n+4*k	header of function for derivatives of order k+1;
   5+5+PN*n+4*k	#defines for formal parameters;
   5+6+PN*n+4*k	tail of function for derivatives of order k+1;
   5+7+PN*n+4*k	#undefs for foral paramenters.
*/

Entry(Int2) autodif(Param param, Uint1Ptr d) {
  CharPtr u;
  int moad;		/* max order of analitical derivatives */
  int ads;
#if ED
  int i;
  char buf[10];
#endif	/* ED */
  OutFile=param[0];	/* output C file name (with C functions for derivatives) */
			/* param[1] is intermidiate C++ file name */
  FuncAttr=param[2];	/* functions attribute (type ATTRIBUTES name(params) {...} */
			/* param[3] is a list of functions to be differentiated */
			/* param[4] is reserved */
  FuncPtr=param[7];	/* local functions of Rhs */

#if ED
  tfn=0;
  unlink(GetParamString(SFS_COMPILER,"EXTDIFFB"));	/* erase ExtDiff batch/command file */
  
  for (moad=0, i=MAXDER-1; i>=0; i--)
    if (d[i]==DER_SYM) {moad=i+1; break;}
#else	/* ED */
  moad=d[2]==DER_SYM ? 3 :(d[1]==DER_SYM ? 2 : (d[0]==DER_SYM ? 1 : 0));
#endif	/* ED */

  out=fopen(param[1],"wt");
  if (!out) return CD_OOPEN;
  if (moad) {
    fprintf(out,Header);
#if ED
    for (i=0; i<MAXDER; i++) {
      sprintf(buf,"D%i_",i+1);
      prefix[i]=StringSave(buf);
    }
#else	/* ED */
    prefix[0]="D1_";
    prefix[1]="D2_";
    prefix[2]="D3_";
#endif	/* ED */
    ads=0;
    for (level=1; level<=moad; level++) {
      if (d[level-1]!=DER_SYM) continue;
      ads|=1<<(level-1);
      undefs=BSNew(1);
      UserDefFuncs(level);
      fn=0;
      EnumNames(param[3],ProcessOneFunction,param);
      u=(CharPtr)BSMerge(undefs,NULL);
      fputs(u,out);
      MemFree(u);
      BSFree(undefs);
    }

    l1=MemGet(1000,FALSE);
    l2=MemGet(1000,FALSE);
    l3=MemGet(1000,FALSE);
#if ED
    l4=MemGet(1000,FALSE);
    l5=MemGet(1000,FALSE);
    *l4=*l5='\0';
#endif	/* ED */
    *l1=*l2=*l3='\0';

    EnumNames(param[3],AddOneName,NULL);
    

#if ED
    fprintf(out,Tail,ads,OutFile,l1,l2,l3,l4,l5);
#else	/* ED */
    fprintf(out,Tail,ads,OutFile,l1,l2,l3);
#endif	/* ED */

    MemFree(l1);
    MemFree(l2);
    MemFree(l3);
#if ED
    MemFree(l4);
    MemFree(l5);
    for (i=0; i<MAXDER; i++)
      MemFree(prefix[i]);
#endif	/* ED */

#if ED
    if (tfn) {
      fclose(edif);
#if _UNIX
      chmod(GetParamString(SFS_COMPILER,"EXTDIFFB"),0700);  /* rwx --- --- (0700 is octal!) */
#endif	/* _UNIX */
    }
#endif	/* ED */

  }
  
  fclose(out);
  return 0;
}