/* -*- mode:C++ ; compile-command: "g++-3.4 -DHAVE_CONFIG_H -I. -I..  -DIN_GIAC -g -c derive.cc" -*- */
#include "giacPCH.h"
/*
 *  Copyright (C) 2000,14 B. Parisse, Institut Fourier, 38402 St Martin d'Heres
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
using namespace std;
#include <stdexcept>
#include <cmath>
#include "derive.h"
#include "usual.h"
#include "symbolic.h"
#include "unary.h"
#include "poly.h"
#include "sym2poly.h" // for equalposcomp
#include "tex.h"
#include "prog.h"
#include "intg.h"
#include "subst.h"
#include "plot.h"
#include "modpoly.h"
#include "moyal.h"
#include "alg_ext.h"
#include "giacintl.h"

#ifndef NO_NAMESPACE_GIAC
namespace giac {
#endif // ndef NO_NAMESPACE_GIAC

   gen eval_before_diff(const gen & expr,const gen & variable,GIAC_CONTEXT){
    identificateur tmp_x_id(" eval_before_diff_x");
    gen tmp_x(tmp_x_id);
    gen res=subst(expr,variable,tmp_x,false,contextptr); // replace variable by a non affected identifier
    gen save_vx_var=vx_var;
    if (variable==vx_var) vx_var=tmp_x;
    int m=calc_mode(contextptr);
    calc_mode(0,contextptr);
    res=eval(res,1,contextptr); // eval res (all identifiers except X will be replaced by their values)
    res=eval(res,1,contextptr); // eval res (all identifiers except X will be replaced by their values)
    calc_mode(m,contextptr);
    vx_var=save_vx_var;
    res=subst(res,tmp_x,variable,false,contextptr);
    return res;
  }

  bool depend(const gen & g,const identificateur & i){
    if (g.type==_IDNT)
      return *g._IDNTptr==i;
    if (g.type==_SYMB)
      return depend(g._SYMBptr->feuille,i);
    if (g.type!=_VECT)
      return false;
    const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end();
    for (;it!=itend;++it){
      if (depend(*it,i))
	return true;
    }
    return false;
  }

  static int count_noncst(const gen & g,const identificateur & i){
    if (g.type!=_VECT)
      return depend(g,i)?1:0;
    int res=0;
    for (unsigned j=0;j<g._VECTptr->size();++j){
      if (depend((*g._VECTptr)[j],i))
	++res;
    }
    return res;
  }

  static gen derive_SYMB(const gen &g_orig,const identificateur & i,GIAC_CONTEXT){
    const symbolic & s = *g_orig._SYMBptr;
    if (s.sommet==at_pnt){
      gen f=g_orig._SYMBptr->feuille;
      if (f.type==_VECT && !f._VECTptr->empty()){
	vecteur v=*f._VECTptr;
	v[0]=derive(v[0],i,contextptr);
	f=gen(v,f.subtype);
	return symbolic(at_pnt,f);
      }
    }
    // if s does not depend on i return 0
    if (!depend(g_orig,i))
      return zero;
    // rational operators are treated first for efficiency
    if (s.sommet==at_plus){
      bool do_step=step_infolevel(contextptr)>1 && count_noncst(s.feuille,i)>1;
      if (do_step)
	gprintf(gettext("Derivative of %gen apply linearity: (u+v+...)'=u'+v'+..."),makevecteur(s),contextptr);
      if (s.feuille.type!=_VECT)
	return derive(s.feuille,i,contextptr);
      vecteur::const_iterator iti=s.feuille._VECTptr->begin(),itend=s.feuille._VECTptr->end();
      int taille=int(itend-iti);
      if (taille==2)
	return derive(*iti,i,contextptr)+derive(*(iti+1),i,contextptr);
      vecteur v;
      v.reserve(taille);
      gen e;
      for (;iti!=itend;++iti){
	e=derive(*iti,i,contextptr);
	if (is_undef(e))
	  return e;
	if (!is_zero(e))
	  v.push_back(e);
      }
      if (v.size()==1)
	return v.front();
      if (v.empty())
	return zero;
      gen res=_plus(gen(v,_SEQ__VECT),contextptr); // symbolic(at_plus,v);
      if (do_step)
	gprintf(gettext("Hence derivative of %gen by linearity is %gen"),makevecteur(g_orig,res),contextptr);
      return res;
    }
    if (s.sommet==at_prod){
      bool do_step=step_infolevel(contextptr)>1 && count_noncst(s.feuille,i)>1;
      if (s.feuille.type==_VECT && s.feuille._VECTptr->size()==2 && s.feuille._VECTptr->back().is_symb_of_sommet(at_inv) && !is_constant_wrt(s.feuille._VECTptr->back()._SYMBptr->feuille,i,contextptr)){
	gen u=s.feuille._VECTptr->front(),v=s.feuille._VECTptr->back()._SYMBptr->feuille;
	if (do_step)
	  gprintf(gettext("Derivative of %gen/%gen, a quotient: (u/v)'=(u'*v-u*v')/v^2"),makevecteur(u,v),contextptr);
	return (derive(u,i,contextptr)*v-u*derive(v,i,contextptr))/pow(v,2,contextptr);
      }
      if (do_step)
	gprintf(gettext("Derivative of %gen, apply product rule: (u*v*...)'=u'*v*...+u*v'*...+..."),makevecteur(s),contextptr);
      if (s.feuille.type!=_VECT)
	return derive(s.feuille,i,contextptr);
      vecteur::const_iterator itbegin=s.feuille._VECTptr->begin(),itj,iti,itend=s.feuille._VECTptr->end();
      int taille=int(itend-itbegin);
      // does not work because of is_linear_wrt e.g. for cos(3*pi/4)
      // if (taille==2) return derive(*itbegin,i,contextptr)*(*(itbegin+1))+(*itbegin)*derive(*(itbegin+1),i,contextptr);
      vecteur v,w;
      v.reserve(taille);
      w.reserve(taille);
      gen e;
      for (iti=itbegin;iti!=itend;++iti){
	w.clear();
	e=derive(*iti,i,contextptr);
	if (is_undef(e))
	  return e;
	if (!is_zero(e)){
	  for (itj=itbegin;itj!=iti;++itj)
	    w.push_back(*itj);
	  w.push_back(e);
	  ++itj;
	  for (;itj!=itend;++itj)
	    w.push_back(*itj);
	  v.push_back(_prod(w,contextptr));
	}
      }
      if (v.size()==1)
	return v.front();
      if (v.empty())
	return zero;
      gen res=symbolic(at_plus,gen(v,_SEQ__VECT));
      if (do_step)
	gprintf(gettext("Hence derivative of %gen by product rule is %gen"),makevecteur(g_orig,res),contextptr);
      return res;
    }
    if (s.sommet==at_neg)
      return -derive(s.feuille,i,contextptr);
    if (s.sommet==at_pow){
      if (s.feuille.type!=_VECT || s.feuille._VECTptr->size()!=2)
	return gensizeerr(contextptr);
      gen base = s.feuille._VECTptr->front(),exponent=s.feuille._VECTptr->back();
      if (step_infolevel(contextptr)>1){
	if (is_constant_wrt(exponent,i,contextptr))
	  gprintf(gettext("Derivative of a power: (%gen)'=(%gen)*(%gen)'*%gen"),makevecteur(symb_pow(base,exponent),exponent,base,symb_pow(base,exponent-1)),contextptr);
	else
	  gprintf(gettext("Derivative of a power: (%gen)'=%gen*(%gen)'*ln(%gen)+(%gen)*(%gen)'*%gen"),makevecteur(symb_pow(base,exponent),symb_pow(base,exponent),exponent,base,exponent,base,symb_pow(base,exponent-1)),contextptr);
      }
      gen dbase=derive(base,i,contextptr),dexponent=derive(exponent,i,contextptr);
      // diff(base^exponent)=diff(exp(exponent*ln(base)))
      // =base^exponent*diff(exponent)*ln(base)+base^(exponent-1)*exponent*diff(base)
      gen expm1=exponent+gen(-1);
      if (is_zero(dexponent))
	return exponent*dbase*pow(base,expm1,contextptr);
      return dexponent*ln(base,contextptr)*s+exponent*dbase*pow(base,expm1,contextptr);
    }
    if (s.sommet==at_inv){
      if (step_infolevel(contextptr)>1)
	gprintf(gettext("Derivative of inv(u)=-u'/u^2 with u=%gen"),makevecteur(s.feuille),contextptr);
      if (s.feuille.is_symb_of_sommet(at_pow)){
	gen & f = s.feuille._SYMBptr->feuille;
	if (f.type==_VECT && f._VECTptr->size()==2)
	  return derive(symb_pow(f._VECTptr->front(),-f._VECTptr->back()),i,contextptr);
      }
      return rdiv(-derive(s.feuille,i,contextptr),pow(s.feuille,2),contextptr);
    }
    if (equalposcomp(inequality_tab,s.sommet))
      return 0;
    if (s.sommet==at_fsolve && s.feuille.type==_VECT && s.feuille._VECTptr->size()>=2){
      vecteur v=*s.feuille._VECTptr;
      if (v[1].is_symb_of_sommet(at_equal) && v[1]._SYMBptr->feuille.type==_VECT && !v[1]._SYMBptr->feuille._VECTptr->empty())
	v[1]=v[1]._SYMBptr->feuille._VECTptr->front();
      gen eq=remove_equal(v[0]),y=v[1],x=i; // fsolve(eq(x,y),y) -> y(x), dy/dx=-(deq/dx)/(deq/dy)
      gen res=-derive(eq,x,contextptr)/derive(eq,y,contextptr);
      res=subst(res,y,s,false,contextptr);
      return res;
    }
    if (s.sommet==at_rootof){
      gen f=s.feuille;
      if (f.type==_VECT && f._VECTptr->size()==2 && f._VECTptr->front().type==_VECT && f._VECTptr->back().type==_VECT){
	vecteur P=*f._VECTptr->front()._VECTptr;
	vecteur Q=*f._VECTptr->back()._VECTptr;
	// d/dx(P(y))=(dP/dx)(y) + (dP/dy) (dy/dx)
	//           =(dP/dx)(y) - (dP/dy) (dQ/dx)/(dQ/dy)
	// where y dependency is in the list polynomial P/Q
	gen Px=trim(*derive(P,i,contextptr)._VECTptr,0);
	Px=symb_rootof(Px,Q,contextptr);
	gen Qx=trim(*derive(Q,i,contextptr)._VECTptr,0);
	Qx=symb_rootof(Qx,Q,contextptr);
	gen Py=derivative(P);
	Py=symb_rootof(Py,Q,contextptr);
	gen Qy=derivative(Q);
	Qy=symb_rootof(Qy,Q,contextptr);
	gen res=Px-(Py*Qx)/Qy;
	res=normal(res,contextptr);
	return res;
      }
      return gensizeerr(gettext("Derivative of rootof currently not handled"));
    }
    if (step_infolevel(contextptr)>1 && s.feuille.type!=_VECT){
      if (s.feuille==i){
	int save_step=step_infolevel(contextptr);
	step_infolevel(contextptr)=0;
	gen der=derive_SYMB(g_orig,i,contextptr);
	step_infolevel(contextptr)=save_step;
	gprintf(gettext("Derivative of elementary function %gen is %gen"),makevecteur(g_orig,der),contextptr);
      }
      else
	gprintf(gettext("Derivative of a composition: (%gen)'=(%gen)'*f'(%gen) where f=%gen"),makevecteur(g_orig,s.feuille,s.feuille,s.sommet),contextptr);
    }
    if (s.sommet==at_UTPT){
      if (s.feuille.type!=_VECT || s.feuille._VECTptr->size()!=2)
	return gensizeerr(contextptr);
      gen & arg=s.feuille._VECTptr->back();
      return -derive(arg,i,contextptr)*_student(s.feuille,contextptr);
    }
    if (s.sommet==at_UTPC){
      if (s.feuille.type!=_VECT || s.feuille._VECTptr->size()!=2)
	return gensizeerr(contextptr);
      gen & arg=s.feuille._VECTptr->back();
      return -derive(arg,i,contextptr)*_chisquare(s.feuille,contextptr);
    }
    if (s.sommet==at_UTPF){
      if (s.feuille.type!=_VECT || s.feuille._VECTptr->size()!=3)
	return gensizeerr(contextptr);
      gen & arg=s.feuille._VECTptr->back();
      return -derive(arg,i,contextptr)*_snedecor(s.feuille,contextptr);
    }
    if (s.sommet==at_program){
      return gensizeerr(gettext("Expecting an expression, not a function"));
    }
    if (s.sommet==at_ln){ 
      if (s.feuille.is_symb_of_sommet(at_abs) )
	return rdiv(derive(s.feuille._SYMBptr->feuille,i,contextptr),s.feuille._SYMBptr->feuille,contextptr);
      if (s.feuille.is_symb_of_sommet(at_inv))
	return -derive(symbolic(at_ln,s.feuille._SYMBptr->feuille),i,contextptr);
      if (s.feuille.is_symb_of_sommet(at_prod)){
	gen res;
	const gen &f=s.feuille._SYMBptr->feuille;
	if (f.type==_VECT){
	  const_iterateur it=f._VECTptr->begin(),itend=f._VECTptr->end();
	  for (;it!=itend;++it)
	    res=res+derive(symbolic(at_ln,*it),i,contextptr);
	  return res;
	}
      }
    }
    if (s.feuille.type==_VECT){
      vecteur v=*s.feuille._VECTptr;
      int vs=int(v.size());
      if (vs>=3 && (s.sommet==at_ifte || s.sommet==at_when) ){
	for (int j=1;j<vs;++j){
	  gen & tmp=v[j];
	  tmp=derive(tmp,i,contextptr); // v[j]=derive(v[j],i,contextptr);
	  // if (is_undef(tmp)) return tmp; 
	  // commented otherwise diff(when(x<0,x^2+3,undef)) returns undef
	}
	return symbolic(s.sommet,gen(v,s.feuille.subtype));
      }
      if (s.sommet==at_piecewise){
	for (int j=0;j<vs/2;++j){
	  gen & tmp=v[2*j+1];
	  tmp=derive(tmp,i,contextptr); // v[2*j+1]=derive(v[2*j+1],i,contextptr);
	}
	if (vs%2){
	  gen & tmp=v[vs-1];
	  tmp=derive(tmp,i,contextptr); // v[vs-1]=derive(v[vs-1],i,contextptr);
	}
	return symbolic(s.sommet,gen(v,s.feuille.subtype));
      }
      if (vs==2 && s.sommet==at_NTHROOT){
	gen base = v[1],exponent=inv(v[0],contextptr);
	gen dbase=derive(base,i,contextptr),dexponent=derive(exponent,i,contextptr);
	// diff(base^exponent)=diff(exp(exponent*ln(base)))
	// =base^exponent*diff(exponent)*ln(base)+base^(exponent-1)*exponent*diff(base)
	if (is_zero(dexponent))
	  return exponent*dbase*s/v[1];
	return dexponent*ln(base,contextptr)*s+exponent*dbase*s/v[1];
      }
      if (vs>=3 && s.sommet==at_Beta){
	gen v0=v[0],v1=v[1],v2=v[2]; 
	if (!is_zero(derive(v0,i,contextptr)) || !is_zero(derive(v1,i,contextptr)) )
	  return gensizeerr("diff of incomplete beta with respect to non constant 1st or 2nd arg not implemented");
	// diff/v2 of int_0^v2 t^(v0-1)*(1-t)^(v1-1) dt
	gen tmp=pow(v2,v0-1,contextptr)*pow(1-v2,v1-1,contextptr)*derive(v2,i,contextptr);
	if (vs==4){
	  gen v3=v[3];
	  if (is_one(v3))
	    return tmp/Beta(v0,v1,contextptr);
	  return gensizeerr(contextptr);
	  gen v3p=derive(v3,i,contextptr);
	  if (!is_zero(v3p))
	    return tmp-pow(v3,v0-1,contextptr)*pow(1-v3,v1-1,contextptr)*v3p;
	}
	return tmp;
      }
      if (vs==4 && s.sommet==at_sum){
	gen v0=v[0],v1=v[1],v2=v[2],v3=v[3];
	if (!is_zero(derive(v1,i,contextptr)) || !is_zero(derive(v2,i,contextptr)) || ! is_zero(derive(v3,i,contextptr)) )
	  return gensizeerr("diff of sum with boundaries or mute variable depending on differentiation variable");
	if (is_inf(v2) || is_inf(v3))
	  *logptr(contextptr) << "Warning, assuming derivative commutes with infinite sum" << endl;
	return _sum(makesequence(derive(v0,i,contextptr),v1,v2,v3),contextptr);
      }
      if ( (vs==2 || (vs==3 && is_zero(v[2]))) && (s.sommet==at_upper_incomplete_gamma || s.sommet==at_lower_incomplete_gamma || s.sommet==at_Gamma)){
	gen v0=v[0],v1=v[1]; 
	if (!is_zero(derive(v0,i,contextptr)))
	  return gensizeerr("diff of incomplete gamma with respect to non constant 1st arg not implemented");
	// diff(int_v1^inf exp(-t)*t^(v0-1) dt)
	gen tmp1=exp(-v1,contextptr)*pow(v1,v0-1,contextptr)*derive(v1,i,contextptr);
	return (s.sommet==at_lower_incomplete_gamma)?tmp1:-tmp1;
      }
      if (vs==3 && (s.sommet==at_upper_incomplete_gamma || s.sommet==at_lower_incomplete_gamma || s.sommet==at_Gamma)){
	return derive(symbolic(s.sommet,makesequence(v[0],v[1]))/symbolic(at_Gamma,v[0]),i,contextptr);
      }
    }
    // now look at other operators, first onearg operator
    if (s.sommet.ptr()->D){
      if (s.feuille.type!=_VECT)
	return derive(s.feuille,i,contextptr)*(*s.sommet.ptr()->D)(1)(s.feuille,contextptr);
      // multiargs operators
      int taille=int(s.feuille._VECTptr->size());
      vecteur v;
      v.reserve(taille);
      vecteur::const_iterator iti=s.feuille._VECTptr->begin(),itend=s.feuille._VECTptr->end();
      gen e;
      for (int j=1;iti!=itend;++iti,++j){
	e=derive(*iti,i,contextptr);
	if (is_undef(e))
	  return e;
	if (!is_zero(e))
	  v.push_back(e*(*s.sommet.ptr()->D)(j)(s.feuille,contextptr));
      }
      if (v.size()==1)
	return v.front();
      if (v.empty())
	return zero;
      return symbolic(at_plus,gen(v,_SEQ__VECT));
    }
    // integrate
    if (s.sommet==at_integrate || s.sommet==at_HPINT){
      if (s.feuille.type!=_VECT)
	return s.feuille;
      vecteur v=*s.feuille._VECTptr;
      int nargs=int(v.size());
      if (nargs<=1)
	return s.feuille;
      if (nargs==2 && is_equal(v[1])){
	gen v1f=v[1]._SYMBptr->feuille;
	if (v1f.type==_VECT && v1f._VECTptr->size()==2){
	  gen v1f1=v1f._VECTptr->front();
	  gen v1f2=v1f._VECTptr->back();
	  v[1]=v1f1;
	  if (v1f2.is_symb_of_sommet(at_interval)){
	    v.push_back(v1f2._SYMBptr->feuille._VECTptr->front());
	    v.push_back(v1f2._SYMBptr->feuille._VECTptr->back());
	    nargs=4;
	  }
	}
      }
      gen res,newint;
      if (v[1]==i)
	res=v[0];
      else {
	res=subst(v[0],v[1],i,false,contextptr);
	newint=derive(v[0],i,contextptr);	 
	if (nargs<4)
	  newint=integrate_gen(newint,v[1],contextptr);
      }
      if (nargs==2)
	return res+newint;
      if (nargs==3)
	return derive(v[2],i,contextptr)*subst(res,i,v[2],false,contextptr);
      if (nargs==4){
	gen a3=derive(v[3],i,contextptr);
	gen b3=is_zero(a3)?zero:limit(res,i,v[3],-1,contextptr);
	gen a2=derive(v[2],i,contextptr);
	gen b2=is_zero(a2)?zero:limit(res,i,v[2],1,contextptr);
	return a3*b3-a2*b2+_integrate(gen(makevecteur(newint,v[1],v[2],v[3]),_SEQ__VECT),contextptr);
      }
      return gensizeerr(contextptr);
    }
    if (s.sommet==at_of && s.feuille.type==_VECT && s.feuille._VECTptr->size()==2){
      // assuming we do not have an index in a list or matrix!
      gen f=s.feuille._VECTptr->front();
      gen arg=s.feuille._VECTptr->back();
      gen darg=derive(arg,i,contextptr);
      if (!is_one(darg)){
	if (darg.type==_VECT){
	  gen res=0;
	  for (int i=0;i<int(darg._VECTptr->size());++i){
	    gen fprime=symbolic(at_derive,makesequence(f,i));
	    res += darg[i]*symbolic(at_of,makesequence(fprime,arg));
	  }
	  return res;
	}
	// f(arg)'=arg'*f'(arg)
	gen fprime=symbolic(at_derive,f);
	return darg*symbolic(at_of,makesequence(fprime,arg));
      }
    }
    // multi derivative and multi-indice derivatives
    if (s.sommet==at_derive){
      if (s.feuille.type!=_VECT)
	return symbolic(at_derive,gen(makevecteur(s.feuille,vx_var,2),_SEQ__VECT));
      if (s.feuille._VECTptr->size()==2){ // derive(f,x)
	gen othervar=(*s.feuille._VECTptr)[1];
	if (othervar.type!=_IDNT) return gensizeerr(gettext("derive.cc/derive_SYMB"));
	if (*othervar._IDNTptr==i){ // _FUNCnd derivative
	  vecteur res(*s.feuille._VECTptr);
	  symbolic sprime(s);
	  res.push_back(2);
	  return symbolic(at_derive,gen(res,_SEQ__VECT));
	}
	else {
	  vecteur var;
	  var.push_back(othervar);
	  var.push_back(i);
	  vecteur nderiv;
	  nderiv.push_back(1);
	  nderiv.push_back(1);
	  return symbolic(at_derive,gen(makevecteur((*s.feuille._VECTptr)[0],var,nderiv),_SEQ__VECT));
	}
      }
      else { // derive(f,x,n)
	if (s.feuille._VECTptr->size()!=3)  return gensizeerr(gettext("derive.cc/derive_SYMB"));
	gen othervar=(*s.feuille._VECTptr)[1];
	if (othervar.type==_IDNT){
	  if (*othervar._IDNTptr==i){ // n+1 derivative
	    vecteur vprime=(*s.feuille._VECTptr);
	    vprime[2] += 1;
	    return symbolic(s.sommet,gen(vprime,_SEQ__VECT));
	  }
	  else {
	    vecteur var;
	    var.push_back(othervar);
	    var.push_back(i);
	    vecteur nderiv;
	    nderiv.push_back((*s.feuille._VECTptr)[2]);
	    nderiv.push_back(1);
	    return symbolic(at_derive,gen(makevecteur((*s.feuille._VECTptr)[0],var,nderiv),_SEQ__VECT));
	  }
	} // end if othervar.type==_IDNT
	else { // othervar.type must be _VECT
	  if (othervar.type!=_VECT)  return gensizeerr(gettext("derive.cc/derive_SYMB"));
	  gen nder((*s.feuille._VECTptr)[2]);
	  if (nder.type!=_VECT ||
	      nder._VECTptr->size()!=othervar._VECTptr->size())  return gensizeerr(gettext("derive.cc/derive_SYMB"));
	  vecteur nderiv(*nder._VECTptr);
	  int pos=equalposcomp(*othervar._VECTptr,i);
	  if (pos){
	    nderiv[pos-1]=nderiv[pos-1]+1;
	  }
	  else {
	    othervar._VECTptr->push_back(i);
	    nderiv.push_back(1);
	  }
	  return symbolic(at_derive,gen(makevecteur((*s.feuille._VECTptr)[0],othervar,nderiv),_SEQ__VECT));
	}
      }
    }
    if (s.sommet==at_re || s.sommet==at_im || s.sommet==at_conj){
      return s.sommet(derive(s.feuille,i,contextptr),contextptr);
    }
    // no info about derivative
    return symbolic(at_derive,gen(makevecteur(s,i),_SEQ__VECT));
    //i.dbgprint();
    //s.dbgprint();
  }

  static gen derive_VECT(const vecteur & v,const identificateur & i,GIAC_CONTEXT){
    vecteur w;
    w.reserve(v.size());
    vecteur::const_iterator it=v.begin(),itend=v.end();
    for (;it!=itend;++it){
      gen tmp=derive(*it,i,contextptr);
      if (is_undef(tmp))
	return tmp;
      w.push_back(tmp);
    }
    return w;
  }

  gen derive(const gen & e,const identificateur & i,GIAC_CONTEXT){
    if (abs_calc_mode(contextptr)==38 && i.id_name[0]>='A' && i.id_name[0]<='Z'){
      identificateur tmp("xdiff");
      gen ee=subst(e,i,tmp,true,contextptr);
      ee=eval(ee,1,contextptr);
      ee=subst(ee,i,tmp,true,contextptr);
      ee=derive(ee,tmp,contextptr);
      ee=subst(ee,tmp,i,true,contextptr);
      return ee;
    }
    switch (e.type){
    case _INT_: case _DOUBLE_: case _ZINT: case _CPLX: case _MOD: case _REAL: case _USER: case _FLOAT_:
      return 0;
    case _IDNT:
      if (is_undef(e))
	return e;
      if (*e._IDNTptr==i)
	return 1;
      else
	return 0;
    case _SYMB:
      return derive_SYMB(e,i,contextptr);
    case _VECT: {
      gen res=derive_VECT(*e._VECTptr,i,contextptr);
      if (res.type==_VECT) res.subtype=e.subtype;
      return res;
    }
    case _FRAC:
      return fraction(derive(e._FRACptr->num,i,contextptr)*e._FRACptr->den-(e._FRACptr->num)*derive(e._FRACptr->den,i,contextptr),e._FRACptr->den);
    default:
      return gentypeerr(contextptr);
    }
    return 0;
  }

  static gen _VECTderive(const gen & e,const vecteur & v,GIAC_CONTEXT){
    vecteur w;
    w.reserve(v.size());
    vecteur::const_iterator it=v.begin(),itend=v.end();
    for (;it!=itend;++it){
      gen tmp=derive(e,*it,contextptr);
      if (is_undef(tmp))
	return tmp;
      w.push_back(tmp);
    }
    return w;    
  }

  static gen derivesymb(const gen& e,const gen & var,GIAC_CONTEXT){
    identificateur x(" x");
    gen xx(x);
    gen f=subst(e,var,xx,false,contextptr);
    f=derive(f,x,contextptr);
    f=subst(f,xx,var,false,contextptr);
    return f;
  }
  gen derive(const gen & e,const gen & vars,GIAC_CONTEXT){
    //  cout << e << " " << vars << endl;
    if (is_equal(e))
      return symb_equal(derive(e._SYMBptr->feuille[0],vars,contextptr),
			derive(e._SYMBptr->feuille[1],vars,contextptr));
    switch (vars.type){
    case _INT_:
      return symbolic(at_derive,makesequence(e,vars));
    case _IDNT:
      return derive(e,*vars._IDNTptr,contextptr);
    case _VECT:
      return _VECTderive(e,*vars._VECTptr,contextptr);
    case _SYMB:
      return derivesymb(e,vars,contextptr);
    default:
      return gensizeerr(contextptr);
    }
    return 0;
  }

  gen derive(const gen & e,const gen & vars,const gen & nderiv,GIAC_CONTEXT){
    if (is_equal(e))
      return symb_equal(derive(e._SYMBptr->feuille[0],vars,nderiv,contextptr),
			derive(e._SYMBptr->feuille[1],vars,nderiv,contextptr));
    if (nderiv.type==_INT_){
      int n=nderiv.val;
      gen ecopie(e),eprime(e);
      int j=1;
      for (;j<=n;++j){
	eprime=ratnormal(derive(ecopie,vars,contextptr),contextptr);
	if (is_undef(eprime))
	  return eprime;
	if ( (eprime.type==_SYMB) && (eprime._SYMBptr->sommet==at_derive))
	  break;
	ecopie=eprime;
      }
      if (j==n+1)
	return eprime;
      return symbolic(at_derive,gen(makevecteur(ecopie,vars,n+1-j),_SEQ__VECT));
    }
    // multi-index derivation
    if (nderiv.type!=_VECT ||
	vars.type!=_VECT)  return gensizeerr(gettext("derive.cc/derive"));
    int s=int(nderiv._VECTptr->size());
    if (s!=signed(vars._VECTptr->size()))  return gensizeerr(gettext("derive.cc/derive"));
    int j=0;
    gen ecopie(e);
    for (;j<s;++j){
      ecopie=derive(ecopie,(*vars._VECTptr)[j],(*nderiv._VECTptr)[j],contextptr);
    }
    return ecopie;
  }  

  symbolic symb_derive(const gen & a,const gen & b){
    return symbolic(at_derive,gen(makevecteur(a,b),_SEQ__VECT));
  }

  gen symb_derive(const gen & a,const gen & b,const gen &c){
    if (is_zero(c))
      return a;
    if (is_one(c))
      return symb_derive(a,b);
    return symbolic(at_derive,gen(makevecteur(a,b,c),_SEQ__VECT));
  }

  gen _derive(const gen & args,GIAC_CONTEXT){
    if (args.type==_STRNG && args.subtype==-1) return  args;
    if (is_equal(args))
      return apply_to_equal(args,_derive,contextptr);
#ifndef NSPIRE
    if (calc_mode(contextptr)==1 && args.type!=_VECT)
      return _derive(makesequence(args,ggb_var(args)),contextptr);
#endif
    vecteur v;
    if (args.type==_VECT && args.subtype==_POLY1__VECT)
      return gen(derivative(*args._VECTptr),_POLY1__VECT);
    if (args.type==_VECT)
      v=plotpreprocess(gen(*args._VECTptr,_SEQ__VECT),contextptr);
    else
      v=plotpreprocess(makesequence(args,vx_var),contextptr);
    if (v.size()>1 && v[1].is_symb_of_sommet(at_unquote))
      v[1]=eval(v[1],1,contextptr);
    if (is_undef(v))
      return v;
    if (step_infolevel(contextptr) && v.size()==2 && v[0].type==_SYMB)
      gprintf(step_derive_header,gettext("===== Derive %gen with respect to %gen ====="),makevecteur(v[0],v[1]),contextptr);
    gen var,res;
    if (args.type!=_VECT && is_algebraic_program(v[0],var,res)){
      if (var.type==_VECT && var.subtype==_SEQ__VECT && var._VECTptr->size()==1)
	var=var._VECTptr->front();
      res=derive(res,var,contextptr);
      return symbolic(at_program,makesequence(var,0,res));
    }
    int s=int(v.size());
    if (s==2){
      if (v[1].type==_VECT && v[1].subtype==_SEQ__VECT){
	vecteur & w=*v[1]._VECTptr;
	int ss=int(w.size());
	gen res=v[0];
	for (int i=0;i<ss;++i)
	  res=ratnormal(derive(res,w[i],contextptr),contextptr);
	return res;
      }
      if (v[0].type==_SPOL1){
	sparse_poly1 res=*v[0]._SPOL1ptr;
	sparse_poly1::iterator it=res.begin(),itend=res.end();
	for (;it!=itend;++it){
	  gen e=it->exponent;
	  it->coeff=it->coeff*e;
	  it->exponent=e-1;
	}
	return res;
      }
      if (args.type!=_VECT && v[0].type==_VECT && v[0].subtype==_POLY1__VECT)
	return gen(derivative(*v[0]._VECTptr),_POLY1__VECT);
      return derive(v[0],v[1],contextptr);
    }
    if (s==3 && (v[2].type==_INT_ || (v[2].type==_VECT && v[2].subtype!=_SEQ__VECT)) )
      return derive( v[0],v[1],v[2],contextptr);    
    if (s<3)
      return gensizeerr(contextptr);
    if (s>=3 && v.back().is_symb_of_sommet(at_equal)){
      gen v_=gen(vecteur(v.begin(),v.end()-1),_SEQ__VECT);
      v_=_derive(v_,contextptr);
      return _subst(makesequence(v_,v.back()),contextptr);
    }
    const_iterateur it=v.begin()+1,itend=v.end();
    res=v[0];
    for (;it!=itend;++it)
      res=ratnormal(_derive(gen(makevecteur(res,*it),_SEQ__VECT),contextptr),contextptr);
    return res;
  }
  // "unary" version
  gen step_derive(const gen & args,GIAC_CONTEXT){
    if (step_infolevel(contextptr))
      ++step_infolevel(contextptr);
    gen res;
#ifndef NO_STDEXCEPT
    try {
      res=_derive(args,contextptr);
    } catch (std::runtime_error & e){
      res=string2gen(e.what(),false);
      res.subtype=-1;
    }
#else
    res=_derive(args,contextptr);
#endif
    if (step_infolevel(contextptr))
      --step_infolevel(contextptr);
    return res;
  }
  gen _diff(const gen & g,GIAC_CONTEXT){
    return _derive(g,contextptr);
  }
  static const char _derive_s []="diff";
  static string printasderive(const gen & feuille,const char * sommetstr,GIAC_CONTEXT){
    if (feuille.type!=_VECT){
      if (feuille.type>=_POLY && feuille.type!=_IDNT)
	return "("+feuille.print()+")'";
      return feuille.print()+"'";
    }
    return sommetstr+("("+feuille.print(contextptr)+")");
  }
  static string texprintasderive(const gen & feuille,const char * sommetstr,GIAC_CONTEXT){
    if (feuille.type!=_VECT)
      return gen2tex(feuille,contextptr)+"'";
    return "\\frac{\\partial \\left("+gen2tex(feuille._VECTptr->front(),contextptr)+"\\right)}{\\partial "+gen2tex(feuille._VECTptr->back(),contextptr)+"}";
  }
  static define_unary_function_eval4_quoted (__derive,&step_derive,_derive_s,printasderive,texprintasderive);
  define_unary_function_ptr5( at_derive ,alias_at_derive,&__derive,_QUOTE_ARGUMENTS,true);

  gen _grad(const gen & args,GIAC_CONTEXT){
    if ( args.type==_STRNG && args.subtype==-1) return  args;
    if (args.type!=_VECT)
      return gensizeerr(contextptr);
    if (args._VECTptr->size()==3){
      gen opt=args._VECTptr->back();
      if (opt.is_symb_of_sommet(at_equal) && (opt._SYMBptr->feuille[0]==at_coordonnees || (opt._SYMBptr->feuille[0].type==_INT_ && opt._SYMBptr->feuille[0].val==_COORDS))){
	gen coord=(*args._VECTptr)[1];
	gen res=_derive(makesequence(args._VECTptr->front(),coord),contextptr);
	if (res.type==_VECT){
	  vecteur resv=*res._VECTptr;
	  if (opt._SYMBptr->feuille[1]==at_sphere && resv.size()==3){
	    resv[1]=resv[1]/coord[0];
	    resv[2]=resv[2]/(coord[0]*sin(coord[1],contextptr));
	    return resv;
	  }
	  if (opt._SYMBptr->feuille[1]==at_cylindre && resv.size()>=2){
	    resv[1]=resv[1]/coord[0];
	    return resv;
	  }
	}
      }
    }
    if (args._VECTptr->size()!=2)
      return gensizeerr(contextptr);
    return _derive(args,contextptr);
  }
  static const char _grad_s []="grad";
  static define_unary_function_eval_quoted (__grad,&_grad,_grad_s);
  define_unary_function_ptr5( at_grad ,alias_at_grad,&__grad,_QUOTE_ARGUMENTS,true);

  gen critical(const gen & g,bool extrema_only,GIAC_CONTEXT){
    gen arg,var;
    if (g.type!=_VECT){
      arg=g;
      var=ggb_var(arg);
    }
    else {
      if (g.subtype!=_SEQ__VECT || g._VECTptr->size()<2)
	return gensizeerr(contextptr);
      arg=g._VECTptr->front();
      var=(*g._VECTptr)[1];
    }
    int savestep=step_infolevel(contextptr);
    gprintf(gettext("===== Critical points for %gen ====="),makevecteur(arg),contextptr);
    step_infolevel(contextptr)=0;
    gen d=_derive(makesequence(arg,var),contextptr);
    gen deq=_equal(makesequence(d,0*var),contextptr);
    // *logptr(contextptr) << "Critical points for "<< arg <<": solving " << deq << " with respect to " << var << endl;
    int c=calc_mode(contextptr);
    calc_mode(0,contextptr);
    gen s=_solve(makesequence(deq,var),contextptr);
    step_infolevel(contextptr)=savestep;
    gprintf(step_extrema1,gettext("Derivative of %gen with respect to %gen is %gen\nSolving %gen with respect to %gen answer %gen"),makevecteur(arg,var,d,deq,var,s.type==_VECT?change_subtype(s,_SEQ__VECT):s),contextptr);
    calc_mode(c,contextptr);
    vecteur ls=lidnt(s);
    for (int i=0;i<int(ls.size());++i){
      if (ls[i]==var || (var.type==_VECT && equalposcomp(*var._VECTptr,ls[i])))
	return gensizeerr("solve error while finding critical points");
    }
    if (s.type==_VECT){
      vecteur res;
      step_infolevel(contextptr)=0;
      gen d2=_derive(makesequence(d,var),contextptr);
      step_infolevel(contextptr)=savestep;
      gprintf(step_extrema2,gettext("Hessian at %gen : %gen"),makevecteur(var,d2),contextptr);
      // *logptr(contextptr) << "Hessian " << d2 << endl;
      vecteur v=*s._VECTptr;
      int vs=int(v.size());
      for (int i=0;i<vs;++i){
	gen g=simplify(subst(d2,var,v[i],false,contextptr),contextptr);
	gprintf(step_extrema3,gettext("Hessian at %gen : %gen"),makevecteur(v[i],g),contextptr);
	// *logptr(contextptr) << "Hessian at " << v[i] << ": " << g << endl;
	if (ckmatrix(g)){
	  g=evalf(g,1,contextptr);
	  if (g.type!=_VECT || !is_numericm(*g._VECTptr,true)){
	    gprintf(step_extrema4,gettext("%gen critical point (unknown type)"),makevecteur(v[i]),contextptr);
	    // *logptr(contextptr) << v[i] << " critical point (unknown type)" << endl;
	    continue;
	  }
	  vecteur w=megvl(*g._VECTptr,contextptr);
	  if (!ckmatrix(w)){
	    gprintf(step_extrema4,gettext("%gen critical point (unknown type)"),makevecteur(v[i]),contextptr);
	    // *logptr(contextptr) << v[i] << " critical point (unknown type)" << endl;
	    continue;
	  }
	  int j=0,ws=int(w.size());
	  for (;j<ws;++j){
	    if (is_zero(w[0][0])){
	      gprintf(step_extrema5,gettext("%gen critical point (0 as eigenvalue) %gen"),makevecteur(v[i],_diag(w,contextptr)),contextptr);
	      // *logptr(contextptr) << v[i] << " critical point (0 as eigenvalue) " << _diag(w,contextptr) << endl;
	      break;
	    }
	    if (is_positive(-w[0][0]*w[j][j],contextptr)){
	      gprintf(step_extrema5,gettext("%gen is a saddle point (2 eigenvalues with opposite sign) %gen"),makevecteur(v[i],_diag(w,contextptr)),contextptr);
	      // *logptr(contextptr) << v[i] << " saddle point (2 eigenvalues with opposite sign) " << _diag(w,contextptr) << endl;
	      break;
	    }
	  }
	  if (j==ws){
	    res.push_back(v[i]);
	    if (is_positive(w[0][0],contextptr)){
	      gprintf(step_extrema6,gettext("%gen is a local minimum %gen"),makevecteur(v[i],_diag(w,contextptr)),contextptr);
	      // *logptr(contextptr) << v[i] << " local minimum " << _diag(w,contextptr) << endl;
	    }
	    else {
	      gprintf(step_extrema6,gettext("%gen is a local maximum %gen"),makevecteur(v[i],_diag(w,contextptr)),contextptr);
	      // *logptr(contextptr) << v[i] << " local maximum " << _diag(w,contextptr) << endl;
	    }
	  }
	  continue;
	} // end multi-dimension
	int d=2;
	gen curd=d2;
	if (is_zero(g)){
	  for (++d;d< NEWTON_DEFAULT_ITERATION;++d){
	    curd=ratnormal(_derive(makesequence(curd,var),contextptr),contextptr);
	    g=simplify(subst(curd,var,v[i],false,contextptr),contextptr);
	    if (!is_zero(g))
	      break;
	  }
	}
	if (d%2==0 && is_strictly_positive(g,contextptr)){
	  gprintf(step_extrema7,gettext("%gen is a local minimum"),makevecteur(v[i]),contextptr);
	  // *logptr(contextptr) << v[i] << " local minimum" << endl;
	  res.push_back(v[i]);
	  continue;
	}
	if (d%2==0 && is_strictly_positive(-g,contextptr)){
	  gprintf(step_extrema7,gettext("%gen is a local maximum"),makevecteur(v[i]),contextptr);
	  // *logptr(contextptr) << v[i] << " local maximum" << endl;
	  res.push_back(v[i]);
	  continue;
	}
	if (d==NEWTON_DEFAULT_ITERATION){
	  gprintf(step_extrema4,gettext("%gen is a critical point (unknown type)"),makevecteur(v[i]),contextptr);
	  //*logptr(contextptr) << v[i] << " critical point (unknown type)" << endl;
	}
	else {
	  gprintf(step_extrema8,gettext("%gen is an inflection point"),makevecteur(v[i]),contextptr);
	  // *logptr(contextptr) << v[i] << " inflection point" << endl;
	}
      }
      if (extrema_only)
	return res;
      else
	return s;
    }
    else
      return s;
  }

#if defined USE_GMP_REPLACEMENTS || defined GIAC_GGB
  gen _extrema(const gen & g,GIAC_CONTEXT){
    return critical(g,true,contextptr);
  }
  static const char _extrema_s []="extrema";
  static define_unary_function_eval_quoted (__extrema,&_extrema,_extrema_s);
  define_unary_function_ptr5( at_extrema ,alias_at_extrema,&__extrema,_QUOTE_ARGUMENTS,true);
#endif

  gen _critical(const gen & g,GIAC_CONTEXT){
    return critical(g,false,contextptr);
  }
  static const char _critical_s []="critical";
  static define_unary_function_eval_quoted (__critical,&_critical,_critical_s);
  define_unary_function_ptr5( at_critical ,alias_at_critical,&__critical,_QUOTE_ARGUMENTS,true);

  // FIXME: This should not use any temporary identifier
  // Should define the identity operator and write again all rules here
  // NB: It requires all D operators for functions to be functions!
  gen _function_diff(const gen & g,GIAC_CONTEXT){
    if ( g.type==_STRNG && g.subtype==-1) return  g;
    if (g.is_symb_of_sommet(at_function_diff)){
      gen & f = g._SYMBptr->feuille;
      return symbolic(at_of,makesequence(gen(symbolic(at_composepow,makesequence(at_function_diff,2))),f));
    }
    if (g.is_symb_of_sommet(at_of)){
      gen & f = g._SYMBptr->feuille;
      if (f.type==_VECT && f._VECTptr->size()==2){
	gen & f1=f._VECTptr->front();
	gen & f2=f._VECTptr->back();
	if (f1.is_symb_of_sommet(at_composepow)){
	  gen & f1f=f1._SYMBptr->feuille;
	  if (f1f.type==_VECT && f1f._VECTptr->size()==2 && f1f._VECTptr->front()==at_function_diff){
	    return symbolic(at_of,makesequence(gen(symbolic(at_composepow,makesequence(at_function_diff,f1f._VECTptr->back()+1))),f2));
	  }
	}
      }
    }
    identificateur _tmpi(" _x");
    gen _tmp(_tmpi);    
    gen dg(derive(g(_tmp,contextptr),_tmp,contextptr));
    if (lop(dg,at_derive).empty()){
      identificateur tmpi(" x");
      gen tmp(tmpi);
      gen res=symb_program(tmp,zero,quotesubst(dg,_tmp,tmp,contextptr),contextptr);
      return res;
    }
    return symbolic(at_function_diff,g);
  }
  static const char _function_diff_s []="function_diff";
  static define_unary_function_eval (__function_diff,&_function_diff,_function_diff_s);
  define_unary_function_ptr5( at_function_diff ,alias_at_function_diff,&__function_diff,0,true);

  static const char _fonction_derivee_s []="fonction_derivee";
  static define_unary_function_eval (__fonction_derivee,&_function_diff,_fonction_derivee_s);
  define_unary_function_ptr5( at_fonction_derivee ,alias_at_fonction_derivee,&__fonction_derivee,0,true);

  gen _implicit_diff(const gen & args,GIAC_CONTEXT){
    if (is_undef(args)) return args;
    if (args.type!=_VECT || (args._VECTptr->size()!=3 && args._VECTptr->size()!=4))
      return gensizeerr(contextptr);
    int ndiff=1;
    if (args._VECTptr->size()==4){
      gen g=args._VECTptr->back();
      if (!is_integral(g) || g.type!=_INT_ || g.val<1)
	return gensizeerr(contextptr);
      ndiff=g.val;
    }
    gen eq(remove_equal(args._VECTptr->front())),x((*args._VECTptr)[1]),y((*args._VECTptr)[2]);
    gen yprime=-derive(eq,x,contextptr)/derive(eq,y,contextptr);
    if (ndiff==1)
      return yprime;
    gen yn=yprime;
    for (int n=2;n<=ndiff;++n){
      yn=ratnormal(derive(yn,x,contextptr)+derive(yn,y,contextptr)*yprime,contextptr);
    }
    return yn;
  }
  static const char _implicit_diff_s []="implicit_diff";
  static define_unary_function_eval (__implicit_diff,&_implicit_diff,_implicit_diff_s);
  define_unary_function_ptr5( at_implicit_diff ,alias_at_implicit_diff,&__implicit_diff,0,true);

  // mode==0 for domain, ==1 for singular values
  void domain(const gen & f,const gen & x,vecteur & eqs,vecteur &excluded,int mode,GIAC_CONTEXT){
    vecteur v=lvarxwithinv(f,x,contextptr);
    lvar(f,v);
    for (int i=0;i<int(v.size());++i){
      gen g=v[i];
      if (g.is_symb_of_sommet(at_nop))
	g=g._SYMBptr->feuille;
      if (is_constant_wrt(g,x,contextptr))
	continue;
      if (g.type!=_SYMB)
	continue;
      gen gf=g._SYMBptr->feuille;
      domain(gf,x,eqs,excluded,mode,contextptr);
      unary_function_ptr & u=g._SYMBptr->sommet;
      if (u==at_inv || u==at_Ei || (mode==1 && (u==at_ln || u==at_log10 || u==at_Ci))){
	excluded=mergevecteur(excluded,gen2vecteur(_solve(makesequence(symb_equal(gf,0),x),contextptr)));
	continue;
      }
      if (u==at_pow){
	if (mode==1){
	  excluded=mergevecteur(excluded,gen2vecteur(_solve(makesequence(symb_equal(gf[0],0),x),contextptr)));
	  continue;
	}
	if (is_greater(gf[1],0,contextptr))
	  eqs.push_back(symb_superieur_egal(gf[0],0));
	else
	  eqs.push_back(symb_superieur_strict(gf[0],0));
	continue;
      }
      if (u==at_ln || u==at_log10 || u==at_Ci){
	eqs.push_back(symb_superieur_strict(gf,0));
	continue;
      }
      if (u==at_acosh){
	if (mode==1)
	  excluded=mergevecteur(excluded,gen2vecteur(_solve(makesequence(symb_equal(gf,0),x),contextptr)));
	else
	  eqs.push_back(symb_superieur_egal(gf,1));
	continue;
      }
      if (u==at_asin || u==at_acos || u==at_atanh){
	if (mode==1)
	  excluded=mergevecteur(excluded,gen2vecteur(_solve(makesequence(symb_equal(pow(gf,2,contextptr),0),x),contextptr)));
	else
	  eqs.push_back(symb_inferieur_egal(pow(gf,2,contextptr),1));
	continue;
      }
      if (u==at_tan){
	excluded=mergevecteur(excluded,gen2vecteur(_solve(makesequence(symb_equal(symb_cos(gf),0),x),contextptr)));
	continue;
      }
      if (u==at_sin || u==at_cos || u==at_exp || u==at_atan)
	continue;
      if (u==at_sinh || u==at_cosh || u==at_tanh)
	continue;
      if (u==at_floor || u==at_ceil || u==at_round || u==at_abs || u==at_sign || u==at_max || u==at_min)
	continue;
      *logptr(contextptr) << g << " function not supported, doing like if it was defined" << endl;
    }
  }
  gen domain(const gen & f,const gen & x,int mode,GIAC_CONTEXT){
    // domain of expression f with respect to variable x
    if (x.type!=_IDNT){
      gen domainx(identificateur("domainx"));
      return domain(subst(f,x,domainx,false,contextptr),domainx,mode,contextptr);
    }
    vecteur eqs,excluded,res;
    bool b=complex_mode(contextptr);
    complex_mode(false,contextptr);
#ifndef NO_STDEXCEPT
    try {
#endif
      domain(f,x,eqs,excluded,mode,contextptr);
      res=gen2vecteur(_solve(makesequence(eqs,x),contextptr));
#ifndef NO_STDEXCEPT
    } catch (std::runtime_error & e ) { *logptr(contextptr) << e.what() << endl;}
#endif
    complex_mode(b,contextptr);
    comprim(excluded);
    if (mode==1)
      return excluded;
    if (excluded.empty())
      return res.size()==1?res.front():res;
    vecteur tmp;
    for (int i=0;i<int(excluded.size());++i){
      tmp.push_back(symbolic(at_different,makesequence(x,excluded[i])));
    }
    if (res.size()==1 && res.front()==x)
      return tmp.size()==1?tmp.front():symbolic(at_and,gen(tmp,_SEQ__VECT));
    else {
      // check if excluded values are solutions inside res
      for (int i=0;i<int(res.size());++i){
	for (int j=0;j<int(excluded.size());++j){
	  gen resi=subst(res[i],x,excluded[j],false,contextptr);
	  resi=eval(resi,1,contextptr);
	  if (is_zero(resi))
	    continue;
	  if (!res[i].is_symb_of_sommet(at_and)){
	    res[i]=symbolic(at_and,makesequence(res[i],tmp));
	    continue;
	  }
	  vecteur v=gen2vecteur(res[i]._SYMBptr->feuille);
	  v.push_back(tmp[j]);
	  res[i]=symbolic(at_and,gen(v,_SEQ__VECT));
	}
      }
      return res;
    }
    // not reached
    if (res.size()==1){
      tmp.insert(tmp.begin(),res.front());
      return symbolic(at_and,gen(tmp,_SEQ__VECT));
    }
    tmp.insert(tmp.begin(),symbolic(at_ou,gen(res,_SEQ__VECT)));
    return symbolic(at_and,gen(tmp,_SEQ__VECT));
  }
  gen _domain(const gen & args,GIAC_CONTEXT){
    if (is_undef(args)) return args;
    if (args.type!=_VECT || args.subtype!=_SEQ__VECT)
      return domain(args,vx_var,0,contextptr);
    vecteur v=*args._VECTptr;
    if (v.size()<2)
      return gensizeerr(contextptr);
    if (is_integral(v[1]))
      v.insert(v.begin()+1,vx_var);
    if (v.size()==2)
      v.push_back(0);
    if (v[2].type!=_INT_)
      return gensizeerr(contextptr);
    return domain(v[0],v[1],v[2].val,contextptr);
  }
  static const char _domain_s []="domain";
  static define_unary_function_eval (__domain,&_domain,_domain_s);
  define_unary_function_ptr5( at_domain ,alias_at_domain,&__domain,0,true);

#ifndef NO_NAMESPACE_GIAC
} // namespace giac
#endif // ndef NO_NAMESPACE_GIAC