Local(FloatHiPtr) J;			/* current set of Jacobian vals */
Local(FloatHiPtr) H;                   	/* current set of Hessian vals */
Local(FloatHiPtr) D3;                   /* current set of third-order dervals */
Local(Int2) Der2num;			/* number of 2nd order derivatives */
Local(Int2) Der3num;			/* number of 3nd order derivatives */

Local(void) CopyToX(Vector x);
Local(void) CopyToP(Vector x);

/**************************************************************************/
/* Returns Jacobian matrix of RHS using symbolic derivatives */
Local(void) SymJac(Vector x, Matrix Jac) {
int SUP,i,j,k,nactive;
VAR s;
  /* x - initial point; F = F^i(x); DF = DF^i(x); DFF; */
  CopyToX(x);
  CopyToP(x);
  J=stagenData->Der1(X,P,&zero);
  nactive=staFunc->GetVectorDim(x)-nphase;
  if (itmap==1) {
    for (i=0; i<nphase; i++) {
      MemCopy(Jac[i],J+(nphase+npar)*i,sizeof(FloatHi)*nphase);
      for (j=0; j<nactive; j++) 
	Jac[i][nphase+j]=*(J+(nphase+npar)*i+nphase+active[j]);
    }
    return;
  }
  MemCopy(DF,J,sizeof(FloatHi)*nphase*(nphase+npar));
  MemCopy(F,x,sizeof(FloatHi)*nphase);

  for (SUP=1; SUP<itmap; SUP++) {
    CopyToX(F);
    stagenData->Rhs(X,P,&zero,F);
    CopyToX(F);
    J=stagenData->Der1(X,P,&zero);
    for (i=0; i<nphase; i++) {
      for (j=0; j<nphase; j++) {
        for (k=0, s=0; k<nphase; k++) s+=(*(J+(nphase+npar)*i+k))*(*(DF+(nphase+npar)*k+j));
        *(DFF+(nphase+npar)*i+j)=s;
      }
    }
    for (i=0; i<nphase; i++) {
      for (j=0; j<npar; j++) {
        for (k=0,s=(*(J+(nphase+npar)*i+nphase+j)); k<nphase; k++) 
	  s+=(*(J+(nphase+npar)*i+k))*(*(DF+(nphase+npar)*k+nphase+j));
        *(DFF+(nphase+npar)*i+nphase+j)=s;
      }
    }
    MemCopy(DF,DFF,sizeof(FloatHi)*nphase*(nphase+npar));
  }

  for (i=0; i<nphase; i++) {
    MemCopy(Jac[i],DF+(nphase+npar)*i,sizeof(FloatHi)*nphase);
    for (j=0; j<nactive; j++) 
         Jac[i][nphase+j]=*(DF+(nphase+npar)*i+nphase+active[j]);
  }
}

/**************************************************************************/
/* analytical second-order derivatives of RHS */

Local(int) Ind2_(int i, int j) {
  return (i<=j) ? ind2_(i,j) : ind2_(j,i);
}

Local(FloatHiPtr) SymHess(Vector x) {
int SUP,i,j,k,l,m,indjk;
VAR s,r;
  /* x - initial point; F = F^i(x); DF = DF^i(x); DFF; */
  CopyToX(x);
  CopyToP(x);
  if (itmap==1) return stagenData->Der2(X,P,&zero);
  J=stagenData->Der1(X,P,&zero);
  MemCopy(DF,J,sizeof(FloatHi)*nphase*(nphase+npar));
  H=stagenData->Der2(X,P,&zero);
  MemCopy(D2,H,sizeof(FloatHi)*nphase*Der2num);
  MemCopy(F,x,sizeof(FloatHi)*nphase);

  for (SUP=1; SUP<itmap; SUP++) {
    CopyToX(F);
    stagenData->Rhs(X,P,&zero,F);
    CopyToX(F);
    J=stagenData->Der1(X,P,&zero);
    H=stagenData->Der2(X,P,&zero);
    /* DxDx */
    for (i=0; i<nphase; i++) {
      for (j=0; j<nphase; j++) {
        for (k=0; k<nphase; k++) {
	  indjk=Ind2_(j,k);
	  s=0;
	  for (l=0; l<nphase; l++) {
	    s+=(*(J+(nphase+npar)*i+l))*(*(D2+Der2num*l+indjk));
	    for (m=0; m<nphase; m++) {
	      s+=(*(H+Der2num*(i)+Ind2_(l,m)))*
	         (*(DF+(nphase+npar)*l+j))*
		 (*(DF+(nphase+npar)*m+k));
	    }
	  }
	  *(DD2+Der2num*i+indjk)=s;
	}
      }
    }
    /* DxDp */
    for (i=0; i<nphase; i++) {
      for (j=0; j<nphase; j++) {
        for (k=0; k<npar; k++) {
	  s=0;
	  for (l=0; l<nphase; l++) {
	    s+=(*(J+(nphase+npar)*i+l))*(*(D2+Der2num*l+ind2_(j,nphase+k)));
	    r=0;
	    for (m=0; m<nphase; m++) {
	      r+=(*(H+Der2num*(i)+Ind2_(l,m)))*
		 (*(DF+(nphase+npar)*m+nphase+k));
	    }
	    r+=*(H+Der2num*(i)+ind2_(l,nphase+k));
	    s+=r*(*(DF+(nphase+npar)*l+j));
	  }
	  *(DD2+Der2num*i+ind2_(j,nphase+k))=s;
	}
      }
    }
    /* DpDp */
    for (i=0; i<nphase; i++) {
      for (j=0; j<npar; j++) {
        for (k=0; k<npar; k++) {
	  indjk=Ind2_(nphase+j,nphase+k);
	  s=(*(H+Der2num*(i)+indjk));
	  for (l=0; l<nphase; l++) {
	    s+=(*(J+(nphase+npar)*i+l))*(*(D2+Der2num*l+indjk));
	    r=0;
	    for (m=0; m<nphase; m++) {
	      r+=(*(H+Der2num*(i)+Ind2_(l,m)))*
		 (*(DF+(nphase+npar)*m+nphase+k));
	    }
	    r+=*(H+Der2num*(i)+ind2_(l,nphase+k));
	    s+=r*(*(DF+(nphase+npar)*l+nphase+j));
	    s+=(*(H+Der2num*(i)+ind2_(l,nphase+j)))*(*(DF+(nphase+npar)*l+nphase+k));
	  }
	  *(DD2+Der2num*i+indjk)=s;
	}
      }
    }
    /* Dx */
    MemCopy(D2,DD2,sizeof(FloatHi)*nphase*Der2num);
    for (i=0; i<nphase; i++) {
      for (j=0; j<nphase; j++) {
        for (k=0, s=0; k<nphase; k++) s+=(*(J+(nphase+npar)*i+k))*(*(DF+(nphase+npar)*k+j));
        *(DFF+(nphase+npar)*i+j)=s;
      }
    }
    /* Dp */
    for (i=0; i<nphase; i++) {
      for (j=0; j<npar; j++) {
        for (k=0,s=(*(J+(nphase+npar)*i+nphase+j)); k<nphase; k++) 
	  s+=(*(J+(nphase+npar)*i+k))*(*(DF+(nphase+npar)*k+nphase+j));
        *(DFF+(nphase+npar)*i+nphase+j)=s;
      }
    }
    MemCopy(DF,DFF,sizeof(FloatHi)*nphase*(nphase+npar));
  }
  return D2;
}

/* Returns second order derivatives of RHS */ 
Local(FloatHi) HessianElem(int i, int j, int k) {

#define HessElem(i,j,k) (*(H+Der2num*(i)+ind2_(j,k)))
#define MixedHessElem(i,j,k) HessElem(i,j,nphase+active[k])  
  
  if (j<=k) 
    if (k<nphase) 
      return HessElem(i,j,k);
    else 
      return MixedHessElem(i,j,k-nphase);
  else 
    if (j<nphase) 
      return HessElem(i,k,j);
    else
      return MixedHessElem(i,k,j-nphase);

}

/**************************************************************************/
/* analytical third-order derivatives of RHS */

#define Swap(a,b) {t=a; a=b; b=t;}

Local(int) Ind3_(int j, int k, int l) {
  int t;
  /* ensure j<=k<=l */
  if (j>k) Swap(j,k);
  if (j>l) Swap(j,l);
  if (k>l) Swap(k,l);
  return ind3_(j,k,l);
}

Local(FloatHiPtr) SymD3(Vector x) {
int SUP,a,i,j,k,l,m,n,indik,indjk,indij,indijk;
VAR s,r;
  /* x - initial point; F = F^i(x); DF = DF^i(x); DFF; */
  CopyToX(x);
  CopyToP(x);
  if (itmap==1) return stagenData->Der3(X,P,&zero);
  J=stagenData->Der1(X,P,&zero);
  MemCopy(DF,J,sizeof(FloatHi)*nphase*(nphase+npar));
  H=stagenData->Der2(X,P,&zero);
  MemCopy(D2,H,sizeof(FloatHi)*nphase*Der2num);
  MemCopy(F,x,sizeof(FloatHi)*nphase);
  D3=stagenData->Der3(X,P,&zero);
  MemCopy(D3F,D3,sizeof(FloatHi)*nphase*Der3num);

  for (SUP=1; SUP<itmap; SUP++) {
    CopyToX(F);
    stagenData->Rhs(X,P,&zero,F);
    CopyToX(F);
    J=stagenData->Der1(X,P,&zero);
    H=stagenData->Der2(X,P,&zero);
    D3=stagenData->Der3(X,P,&zero);
    /* DxDxDx */
    for (a=0; a<nphase; a++) {
      for (i=0; i<nphase; i++) {
        for (j=0; j<nphase; j++) {
	  indij=Ind2_(i,j);
          for (k=0; k<nphase; k++) {
	    indik=Ind2_(i,k);
	    indjk=Ind2_(j,k);
	    indijk=Ind3_(i,j,k);
	    s=0;
	    for (l=0; l<nphase; l++) {
	      s+=(*(J+(nphase+npar)*a+l))*(*(D3F+Der3num*l+indijk));
	      for (m=0; m<nphase; m++) {
		s+=(*(H+Der2num*a+Ind2_(l,m)))*(
		   (*(D2+Der2num*l+indik))*(*(DF+(nphase+npar)*m+j))+
		   (*(D2+Der2num*m+indjk))*(*(DF+(nphase+npar)*l+i))+
		   (*(D2+Der2num*l+indij))*(*(DF+(nphase+npar)*m+k))
		   );
		r=0;
		for (n=0; n<nphase; n++) {
		  r+=(*(D3+Der3num*a+Ind3_(l,m,n)))*(*(DF+(nphase+npar)*n+k));
		}
		r*=(*(DF+(nphase+npar)*l+i))*(*(DF+(nphase+npar)*m+j));
		s+=r;
	      }
	    }
	    *(DD3F+Der3num*a+indijk)=s;
	  }
        }
      }
    }
    MemCopy(D3F,DD3F,sizeof(FloatHi)*nphase*Der3num);
    /* DxDx */
    for (i=0; i<nphase; i++) {
      for (j=0; j<nphase; j++) {
        for (k=0; k<nphase; k++) {
	  indjk=Ind2_(j,k);
	  s=0;
	  for (l=0; l<nphase; l++) {
	    s+=(*(J+(nphase+npar)*i+l))*(*(D2+Der2num*l+indjk));
	    for (m=0; m<nphase; m++) {
	      s+=(*(H+Der2num*(i)+Ind2_(l,m)))*
	         (*(DF+(nphase+npar)*l+j))*
		 (*(DF+(nphase+npar)*m+k));
	    }
	  }
	  *(DD2+Der2num*i+indjk)=s;
	}
      }
    }
    MemCopy(D2,DD2,sizeof(FloatHi)*nphase*Der2num);
    /* Dx */
    for (i=0; i<nphase; i++) {
      for (j=0; j<nphase; j++) {
        for (k=0, s=0; k<nphase; k++) s+=(*(J+(nphase+npar)*i+k))*(*(DF+(nphase+npar)*k+j));
        *(DFF+(nphase+npar)*i+j)=s;
      }
    }
    MemCopy(DF,DFF,sizeof(FloatHi)*nphase*(nphase+npar));
  }
  return D3F;
}

#define Der3Elem(i,j,k,l) (*(D3+Der3num*(i)+ind3_(j,k,l)))

/* Returns third order derivatives */
Local(FloatHi) D3Elem(int i, int j, int k, int l) {
  int t;
  /* ensure j<=k<=l */
  if (j>k) Swap(j,k);
  if (j>l) Swap(j,l);
  if (k>l) Swap(k,l);
  return Der3Elem(i,j,k,l);
}
#undef Swap