/* -*- mode:C++ ; compile-command: "g++-3.4 -I.. -g -c pari.cc" -*- */
#include "giacPCH.h"

/*  PARI interface
 *  Copyright (C) 2001,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/>.
 */
#ifndef NSPIRE
using namespace std;
#endif
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <stdexcept>
#include "gen.h" 
#include "identificateur.h"
#include "sym2poly.h"
#include "plot.h"
#include "prog.h"
#include "usual.h"
#include "input_lexer.h"
#include "modpoly.h"
#include "giacintl.h"
#ifdef USE_GMP_REPLACEMENTS
#undef HAVE_LIBPARI
#endif

#ifdef HAVE_LIBPARI
#ifdef HAVE_PTHREAD_H
#include <pthread.h>
#endif

static long int abs(long int & l){
  if (l<0)
    return -l;
  else
    return l;
}
#include "pari.h"
extern "C" {
#include <pari/pari.h>
#include <pari/paripriv.h>
#ifdef ENABLE_TLS
  extern THREAD void *PARI_stack_limit;
#else
  extern void *PARI_stack_limit;
#endif
  extern entree functions_basic[];
}
#if PARI_VERSION_CODE<PARI_VERSION(2,4,0) // 132096
#define PARI23
#else
jmp_buf env; 
static void
gp_err_recover(long numerr)
{
  longjmp(env, numerr);
}

#endif
#include <cstdlib>

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

#ifdef HAVE_LIBPTHREAD
  static pthread_mutex_t * pari_mutex_ptr = 0;
#endif

#if PARI_VERSION_CODE >= PARI_VERSION(2,8,0)
#define PARI_DYNAMIC_STACK
#endif

  static map<string,entree *> pari_function_table;
  static void do_giac_pari_init(long maxprime){
#ifdef PARI_DYNAMIC_STACK
    long pari_mem_size=512000;
#else
    long pari_mem_size=64000000;
#endif
    if (getenv("PARI_SIZE")){
      string pari_size_s(getenv("PARI_SIZE"));
      pari_mem_size= atoi(pari_size_s.c_str());
    }
    // do not initialize INIT_JMP so that PARI error do not exit
    pari_init_opts(pari_mem_size,maxprime,INIT_SIGm | INIT_DFTm);
#ifdef PARI_DYNAMIC_STACK
    paristack_setsize(pari_mem_size, (1<<30)); // pari 2.8
    //paristack_alloc(pari_mem_size, (1<<30)); // pari 2.8
#endif
    entree * ptr=functions_basic;
    for (;ptr->name;++ptr){
      pari_function_table[ptr->name]=ptr;
    }
    // initialize variable ordering x,y,z
    gp_read_str("[x,y,z,t]");
  }

#if 0
  struct giac_pari_init {
    giac_pari_init(long maxprime) { 
      if(!avma){ 
	do_giac_pari_init(maxprime); 
      }
    }
  };
  static long pari_maxprime=100000;
  static giac_pari_init bidon(pari_maxprime);
#else
  static long pari_maxprime=100000;
  long get_pari_avma() {
    if(!avma){
      do_giac_pari_init(pari_maxprime);
    }
    return avma;
  }
#endif

  static gen pow2sizeof_long(pow(256,sizeof(long)));
  // Conversion of a GEN integer to a gen, using Horner method
  static gen t_INT2gen(const GEN & G){
    long Gs=signe(G);
    if (!Gs)
      return 0;
    //setsigne(G,1);
    long Gpl=lgefint(G)-2;
    // use one of the following code depending how pari codes long integers
    // Code 1
#if !defined(__APPLE__) && !defined(WIN32) || defined(WIN64)
    ref_mpz_t * mz = new ref_mpz_t;
    mpz_realloc2(mz->z,32*Gpl);
    mpz_import(mz->z,Gpl,-1,sizeof(GEN),0,0,&G[2]);
    //setsigne(G,Gs);
    if (Gs>0) 
      return mz;
    else
      return -gen(mz);
#else
    // Code 2
    --Gpl;
    long * Gp=int_MSW(G);
    gen res;
    for (int i=0;i<=Gpl;++i){
#ifdef INT128
      res=res*pow2sizeof_long+int128_t((ulonglong)*Gp);
#else
      res=res*pow2sizeof_long+longlong(unsigned(*Gp));
#endif
      Gp=int_precW(Gp);
    }
    return Gs<0?-res:res;
#endif
  }

  static gen t_REAL2gen(const GEN & G){
    long Gs=signe(G);
    if (!Gs)
      return 0.0;
    long n=lg(G);
    gen res;
    for (int i=2;i<n;++i){
#ifdef x86_64 // FIXME make a gen constructor from ulonglong
      unsigned long u=G[i];
      res=res*pow2sizeof_long;
      if (u%2)
	res += 1;
      res += 2*gen(longlong(u>>1));
#else
      res=res*pow2sizeof_long+longlong(unsigned(G[i]));
#endif
    }
    res=res*pow(plus_two,int(expo(G)+1-bit_accuracy(n)));
    res=_evalf(makesequence(res,int(n/3.3)),context0);
    return Gs<0?-res:res;
  }

  static gen t_POL2gen(const GEN & G,const vecteur & vars){
    if (!signe(G))
      return 0;
    long n=lg(G);
    vecteur res;
    for (long i=2;i<n;++i){
      res.push_back(GEN2gen((GEN)G[i],vars));
    }
    reverse(res.begin(),res.end());
    long vn=varn(G);
    gen x;
    if (vn<long(vars.size())){
      x=vars[vn];
      return symb_horner(res,x);
    }
    if (!vars.empty() && vn==long(vars.size())){
      x=vars[vn-1];
      return symb_horner(res,x);
    }
    else
      return gen(res,_POLY1__VECT);
  }

  static gen t_VEC2gen(const GEN & G,const vecteur & vars,long debut,long fin){
    vecteur res;
    for (long i=debut;i<fin;++i){
      res.push_back(GEN2gen((GEN)G[i],vars));
    }
    return gen(res,typ(G)==t_VEC?0:99);
  }

  static gen t_VECSMALL2gen(const GEN & G){
    long fin=lg(G);
    vecteur res;
    for (long i=1;i<fin;++i){
      res.push_back((int)G[i]);
    }
    return gen(res,98);
  }

  static gen t_MOD2gen(const GEN & G,const vecteur & vars){
    return makemod(GEN2gen((GEN) G[2],vars),GEN2gen((GEN) G[1],vars));
  }

  static gen t_POLMOD2gen(const GEN & G,const vecteur & vars){
    gen tmp=_fxnd(GEN2gen((GEN) G[2],vars),context0);
    gen n(tmp),d(1);
    if (tmp.type==_VECT && tmp._VECTptr->size()==2){
      n=tmp._VECTptr->front();
      d=tmp._VECTptr->back();
    }
    return eval(symbolic(at_rootof,makesequence(n,GEN2gen((GEN) G[1],vars))),1,context0)/d;
    find_or_make_symbol("Mod",tmp,0,false,context0);
    return symbolic(at_of,makesequence(tmp,gen(makevecteur(GEN2gen((GEN) G[2],vars),GEN2gen((GEN) G[1],vars)),_SEQ__VECT)));
  }

  static gen t_COMPLEX2gen(const GEN & G,const vecteur & vars){
    return gen(GEN2gen((GEN) G[1],vars),GEN2gen((GEN) G[2],vars));
  }

  static gen t_FRAC2gen(const GEN & G,const vecteur & vars){
    return fraction(GEN2gen((GEN)G[1],vars),GEN2gen((GEN)G[2],vars));
  }

  static gen t_QUAD2gen(const GEN & G,const vecteur & vars){
    // use w__IDNT_e like pari for all quadratics
    return GEN2gen((GEN) G[1],vars)+w__IDNT_e*GEN2gen((GEN)G[2],vars);
  }

  static gen t_PADIC2gen(const GEN & G,const vecteur & vars){
    gen O;
    find_or_make_symbol("O",O,0,false,context0);
    gen p(GEN2gen((GEN) G[2],vars)),val(longlong(valp(G)));
    return pow(p,val,context0)*(GEN2gen((GEN) G[4],vars)+symbolic(at_of,makesequence(O,symb_pow(p,longlong(precp(G)))))); // removed symb_quote for p-adic
  }

  // WARNING: If g is a matrix this print the transpose of the matrix
  static string GEN2string(const GEN & g){
    // cerr << typ(g) << " " << t_MAT << endl;
    char * ch;
    string s;
    if ((typ(g)==t_MAT) || (typ(g)==t_COL)){
      int taille=lg(g);
      s +="[";
      for (int i=1;i<taille;++i){
	s += GEN2string((long *)g[i]);
	if (i==taille-1)
	  s+="]";
	else
	  s+=",";
      }
      return s;
    }
    ch=GENtostr(g);
    s=ch;
#ifndef HAVE_LIBGC
    free(ch);
#endif
    return s;
  }

  static gen default2gen(const GEN &G){
    string s=GEN2string(G);
    gen g;
    try {
      g=gen(s,context0);
    } catch(...){
      return string2gen(s,false);
    }
    return g;
  }

  gen GEN2gen(const GEN & G,const vecteur & vars){
#ifdef EMCC
    return default2gen(G);
#endif
    switch (typ(G)){
    case t_INT:
      return t_INT2gen(G);
    case t_INTMOD: 
      return t_MOD2gen(G,vars);
    case t_POLMOD:
      return t_POLMOD2gen(G,vars);
    case t_FRAC: case t_RFRAC:
      return t_FRAC2gen(G,vars);
    case t_COMPLEX:
      return t_COMPLEX2gen(G,vars);
    case t_REAL:
      return t_REAL2gen(G);
    case t_POL:
      return t_POL2gen(G,vars);
    case t_VEC: case t_COL:
      return t_VEC2gen(G,vars,1,lg(G));
    case t_VECSMALL: 
      return t_VECSMALL2gen(G);
    case t_MAT:
      return _tran(t_VEC2gen(G,vars,1,lg(G)),context0);
    case t_LIST:
#ifdef PARI23
      return t_VEC2gen(G,vars,2,lgeflist(G));
#else
      return t_VEC2gen(G,vars,2,list_nmax(G));
#endif
    case t_STR:
      return string2gen(GSTR(G),false);
    case t_QUAD:
      return t_QUAD2gen(G,vars);
    case t_PADIC:
      return t_PADIC2gen(G,vars);
    default:
      return default2gen(G);
    }
  }

  static std::string pariprint(const gen & e,int varnum,GIAC_CONTEXT);

  static string pariprint_VECT(const vecteur & v,int varnum,int subtype,GIAC_CONTEXT){
    string s;
    const_iterateur it=v.begin(),itend=v.end();
    if (subtype==_POLY1__VECT){
      if (v.empty())
	return "0";
      string tmp=")*x"+print_INT_(varnum)+"+";
      ++varnum;
      s=pariprint(*it,varnum,contextptr);
      for (++it;it!=itend;++it){
	s="("+s+tmp+pariprint(*it,varnum,contextptr);
      }
      --varnum;
      return s;
    }
    if (subtype!=_SEQ__VECT)
      s="[";
    for (;it!=itend;++it){
      s += pariprint(*it,varnum,contextptr);
      if (it+1!=itend)
	s += ",";      
    }
    if (subtype!=_SEQ__VECT)
      s+="]";
    return s;
  }

  static string pariprintmatrice(const gen & e,int varnum,GIAC_CONTEXT){
    string res = "[";
    const_iterateur it=e._VECTptr->begin(),itend=e._VECTptr->end();
    for (;it!=itend;++it){
      res += pariprint_VECT(*it->_VECTptr,varnum,_SEQ__VECT,contextptr);
      if (it+1!=itend)
	res += ";";
      else
	res +="]";
    }
    return res;
  }

  static string pariprint(const gen & e,int varnum,GIAC_CONTEXT){
    int save_maple_mode=xcas_mode(contextptr);
    xcas_mode(contextptr)=1;
    string res;
    switch (e.type){
    case _INT_:
      res=print_INT_(e.val);
      break;
    case _ZINT:
      res=e.print(contextptr);
      break;
    case _CPLX:
      res=e._CPLXptr->print(contextptr)+"+I*"+(e._CPLXptr+1)->print(contextptr);
      break;
    case _VECT:
      if (ckmatrix(e))
	res=pariprintmatrice(e,varnum,contextptr);
      else
	res=pariprint_VECT(*e._VECTptr,varnum,e.subtype,contextptr);
      break;
    case _SYMB:
      if (e._SYMBptr->sommet==at_rootof && e._SYMBptr->feuille.type==_VECT && e._SYMBptr->feuille._VECTptr->size()==2)
	res="Mod("+pariprint(e._SYMBptr->feuille._VECTptr->front(),varnum,contextptr)+","+pariprint(e._SYMBptr->feuille._VECTptr->back(),varnum,contextptr)+")";
      else {
	if (e._SYMBptr->sommet.ptr()->printsommet) // FIXME
	  res=e.print(contextptr);
	else
	  res=string(e._SYMBptr->sommet.ptr()->s)+"("+pariprint(e._SYMBptr->feuille,varnum,contextptr)+")";
      }
      break;
    case _MOD:
      res= "Mod("+pariprint(*e._MODptr,varnum,contextptr)+","+pariprint(*(e._MODptr+1),varnum,contextptr)+")";
      break;
    case _FRAC:
      res= pariprint(*e._FRACptr,varnum,contextptr)+"/("+pariprint(*(e._FRACptr+1),varnum,contextptr)+")";
      break;
    default: // _SYMB with printsommetasoperator, _DOUBLE_, _REAL, _IDNT
      res=e.print(contextptr);
    }
    xcas_mode(contextptr)=save_maple_mode;
    return res;
  }
  
  static GEN zint2GEN(const gen & g){
    mpz_t * zz=g._ZINTptr;
    int sgn=mpz_sgn(*zz);
    if (!sgn)
      return utoi(0);
    int count=mpz_sizeinbase(*zz,2);
    if (count % (8*sizeof(GEN))) 
      count=count/(8*sizeof(GEN))+3;
    else
      count=count/(8*sizeof(GEN))+2;
    GEN G=cgetg(count,t_INT);
    //size_t countp;
    // mpz_export(&G[2],&countp,-1,sizeof(GEN),0,0,*zz);
    setlgefint(G,count);
    setsigne(G,sgn);
    // return G;
    gen q(abs(g)),tmp,r;
    GEN Gstep(int2n(16));
    vector<int> v;
    for (;!is_zero(q);){
      r=irem(q,gen(1<<16),tmp);
      v.push_back(r.val);
      q=tmp;
    }
    int s=v.size();
    GEN res(utoi(0));
    for (int i=s-1;i>=0;--i){
      res=gmul(res,Gstep);
      res=gadd(res,utoi(v[i]));
    }
    setsigne(res,sgn);
    return res;
  }

  static GEN ingen2GEN(const gen & e,const vecteur & vars,GIAC_CONTEXT);

  static GEN vect2GEN(const gen & g,const vecteur & vars,GIAC_CONTEXT){
    vecteur v (*g._VECTptr);
    int n=v.size(),decal=1;
    GEN res;
    if (g.subtype==_POLY1__VECT){
      decal=2;
      res=cgetg(n+decal,t_POL);
      reverse(v.begin(),v.end());
    }
    else {
      res=cgetg(n+decal,g.subtype==99?t_COL:(g.subtype==98?t_VECSMALL:t_VEC));
    }
    for (int i=0;i<n;++i)
      gel(res,i+decal)=ingen2GEN(v[i],vars,contextptr);
    if (decal==2){
      setsigne(res,1);
      setvarn(res,0);
    }
    return res;
  }

  static GEN mat2GEN(const gen & g,const vecteur & vars,GIAC_CONTEXT){
    matrice M = mtran(*g._VECTptr);
    int n=M.size(),m=M[0]._VECTptr->size();
    GEN res=cgetg(n+1,t_MAT);
    for (int i=1;i<=n;++i){
      GEN resi=gel(res,i)=cgetg(m+1,t_COL);
      vecteur & v = *M[i-1]._VECTptr;
      for (int j=1;j<=m;++j){
	gel(resi,j)=ingen2GEN(v[j-1],vars,contextptr);
      }
    }
    return res;
  }

  static GEN cplx2GEN(const gen & g,const vecteur & vars,GIAC_CONTEXT){
    GEN res=cgetg(3,t_COMPLEX);
    gel(res,1)=ingen2GEN(*g._CPLXptr,vars,contextptr);
    gel(res,2)=ingen2GEN(*(g._CPLXptr+1),vars,contextptr);
    return res;
  }

  static GEN frac2GEN(const gen & g,const vecteur & vars,GIAC_CONTEXT){
    GEN res=cgetg(3,t_FRAC);
    gel(res,1)=ingen2GEN(g._FRACptr->num,vars,contextptr);
    gel(res,2)=ingen2GEN(g._FRACptr->den,vars,contextptr);
    return res;
  }

  static GEN real2GEN(const gen & e_,GIAC_CONTEXT){
    gen e(e_);
    bool neg=false;
    if (is_strictly_positive(-e,contextptr)){
      e=-e;
      neg=true;
    }
    int prec=53;
#ifdef HAVE_LIBMPFR
    if (e.type==_REAL)
      prec=mpfr_get_prec(e._REALptr->inf);
#endif
    const char * s=e.print(contextptr).c_str();
    GEN g=strtor((char *)s,prec);
    if (neg)
      g=gneg(g);
    if (debug_infolevel)
      CERR << "real converted to pari " << GEN2gen(g,vecteur(0)) << endl;
    else e=GEN2gen(g,vecteur(0));
    // for some strange reason, converting g to a gen fixes a bug in conversion
    return g;
  }

  static GEN ingen2GEN(const gen & e,const vecteur & vars,GIAC_CONTEXT){
    switch (e.type){
    case _INT_:
      return stoi(e.val);
    case _ZINT:
      return zint2GEN(e);
    case _CPLX:
      return cplx2GEN(e,vars,contextptr);
    case _FRAC:
      return frac2GEN(e,vars,contextptr);
    case _DOUBLE_: case _REAL:
      return real2GEN(e,contextptr);
    case _VECT:
      if (ckmatrix(e))
	return mat2GEN(e,vars,contextptr);
      else
	return vect2GEN(e,vars,contextptr);
    }
    // add vars to e
    string s=pariprint(e,0,contextptr);
    vecteur vars_(vars);
    // if (vars_.empty()) vars_.push_back(vx_var);
    if (!vars_.empty())
      s="["+(vars_.size()==1?vars_.front().print():print_VECT(vars_,_SEQ__VECT,contextptr))+","+s+"]";
    GEN res= gp_read_str((char *) s.c_str());
    // gp_read_str seems to have problems with large strings (s.size()>2^16?)
    return vars_.empty()?res:gel(res,1+vars_.size());
  }
  GEN gen2GEN(const gen & e,const vecteur & vars,GIAC_CONTEXT){
#ifdef PARI23
    if (setjmp(GP_DATA->env)){
      setsizeerr(gettext("Error in PARI subsystem"));
    }
#else
    cb_pari_err_recover=gp_err_recover;
    if (setjmp(env)){ 
      setsizeerr(gettext("Error in PARI subsystem"));
    }
#endif
    return ingen2GEN(e,vars,contextptr);
  }

  static void pari_cleanup(void * arg) {
#ifdef HAVE_LIBPTHREAD
    if (arg)
      pthread_mutex_unlock((pthread_mutex_t *)arg);
#endif
  }

  int check_pari_mutex(){
#ifdef HAVE_LIBPTHREAD
    if (!pari_mutex_ptr){
      pthread_mutex_t tmp=PTHREAD_MUTEX_INITIALIZER;
      pari_mutex_ptr=new pthread_mutex_t(tmp);
    }
    return pthread_mutex_trylock(pari_mutex_ptr);
#else
    return 0;
#endif
  }

  void abort_if_locked(){
    if (check_pari_mutex())
      setsizeerr(gettext("PARI locked by another thread. Try again later.\nIf PARI is locked by a cancelled thread, you can unlock it with pari_unlock()"));
  }

  gen pari_isprime(const gen & e,int certif){
    gen tmp;
    abort_if_locked();
#ifdef HAVE_LIBPTHREAD
    pthread_cleanup_push(pari_cleanup, (void *) pari_mutex_ptr);
#endif
    long av=get_pari_avma();
    tmp=GEN2gen(gisprime(gen2GEN(e,vecteur(0),0),certif),vecteur(0));
    avma=av;
#ifdef HAVE_LIBPTHREAD
    if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
    pthread_cleanup_pop(0);
#endif
    return tmp;
  }

  string pari_ifactor(const gen & e){
    abort_if_locked();
    string s;
#ifdef HAVE_LIBPTHREAD
    pthread_cleanup_push(pari_cleanup, (void *) pari_mutex_ptr);
#endif
    long av=get_pari_avma();
    GEN g=gen2GEN(e,vecteur(0),0);
    GEN gf=factorint(g,0);
    s=GEN2string(gf);
    avma=av;
#ifdef HAVE_LIBPTHREAD
    if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
    pthread_cleanup_pop(0);
#endif
    return s;
  }

  gen pari_gamma(const gen & e){
    abort_if_locked();
    gen res;
#ifdef HAVE_LIBPTHREAD
    pthread_cleanup_push(pari_cleanup, (void *) pari_mutex_ptr);
#endif
    long av=get_pari_avma();
    GEN g=gen2GEN(e,vecteur(0),0);
    GEN gf=ggamma(g,precision(g));
    res=GEN2gen(gf,vecteur(0));
    avma=av;
#ifdef HAVE_LIBPTHREAD
    if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
    pthread_cleanup_pop(0);
#endif
    return res;
  }

  gen pari_zeta(const gen & e){
    abort_if_locked();
    gen res;
#ifdef HAVE_LIBPTHREAD
    pthread_cleanup_push(pari_cleanup, (void *) pari_mutex_ptr);
#endif
    long av=get_pari_avma();
    GEN g=gen2GEN(e,vecteur(0),0);
    GEN gf=gzeta(g,precision(g));
    res=GEN2gen(gf,vecteur(0));
    avma=av;
#ifdef HAVE_LIBPTHREAD
    if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
    pthread_cleanup_pop(0);
#endif
    return res;
  }

  gen pari_psi(const gen & e){
    abort_if_locked();
    gen res;
#ifdef HAVE_LIBPTHREAD
    pthread_cleanup_push(pari_cleanup, (void *) pari_mutex_ptr);
#endif
    long av=get_pari_avma();
    GEN g=gen2GEN(e,vecteur(0),0);
    GEN gf=gpsi(g,precision(g));
    res=GEN2gen(gf,vecteur(0));
    avma=av;
#ifdef HAVE_LIBPTHREAD
    if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
    pthread_cleanup_pop(0);
#endif
    return res;
  }

  gen pari_ffinit(const gen & p,int n){
    gen tmp;
#ifdef HAVE_LIBPTHREAD
    abort_if_locked();
    pthread_cleanup_push(pari_cleanup, (void *) pari_mutex_ptr);
#endif
    long av=get_pari_avma();
    tmp=GEN2gen(ffinit(gen2GEN(p,vecteur(0),0),n,0),vecteur(0));
    avma=av;
#ifdef HAVE_LIBPTHREAD
    if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
    pthread_cleanup_pop(0);
#endif
    return tmp;
  }

  // for factorization over Z when many modular factors arise
  // This is a call to PARI combine_factors
  // WARNING: You must remove static from the declaration of combine_factors
  // in pari/src/basemath/polarith2.c
  // GEN combine_factors(GEN a, GEN famod, GEN p, long klim, long hint);
  bool pari_lift_combine(const vecteur & a,const vector<vecteur> & factmod,gen & modulo,vector<vecteur> & res){
#ifdef PARI23
    long av=get_pari_avma();
    GEN pari_a=gen2GEN(r2e(a,x__IDNT_e,context0),vecteur(0),0);
    string s("[");
    vector<vecteur>::const_iterator it=factmod.begin(),itend=factmod.end();
    for (;it!=itend;){
      s += r2e(*it,x__IDNT_e,context0).print();
      ++it;
      if (it==itend)
	break;
      s+=",";
    }
    s+="]";
    // cerr << s << endl;
    GEN pari_factmod=gp_read_str((char *) s.c_str());
    GEN pari_modulo=gen2GEN(modulo,vecteur(0),0);
    GEN pari_res=combine_factors(pari_a,pari_factmod,pari_modulo,0,1);
    // back conversion
    string res_s=GENtostr(pari_res);
    gen res_v(res_s.substr(0,res_s.size()-1),context0);
    if (res_v.type!=_VECT)
      setsizeerr();
    const_iterateur jt=res_v._VECTptr->begin(),jtend=res_v._VECTptr->end();
    for (;jt!=jtend;++jt){
      res.push_back(*e2r(*jt,x__IDNT_e,context0)._VECTptr);
    }
    avma=av;
    return true;
#else
    return false;
#endif
  }

  static gen pari_exec(const string & s,GIAC_CONTEXT){
    long av=get_pari_avma();
    void * save_pari_stack_limit = PARI_stack_limit;
    PARI_stack_limit=0; // required since the stack changed
#ifdef PARI23
    if (setjmp(GP_DATA->env))
#else
    cb_pari_err_recover=gp_err_recover;
    if (setjmp(env))
#endif
      { 
#ifdef HAVE_LIBPTHREAD
	if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
#endif
	avma = av;
	*logptr(contextptr) << gettext("Error in PARI subsystem") << endl;
	PARI_stack_limit = save_pari_stack_limit ;
	// setsizeerr();
	return undef;
      } 
    GEN gres= gp_read_str((char *) s.c_str());
    gen res=GEN2gen(gres,vecteur(0));
    avma=av;
    PARI_stack_limit = save_pari_stack_limit ;
    return res;
  }

#include "input_parser.h"
#define _ARGS_ argvec[0], argvec[1], argvec[2], argvec[3],\
               argvec[4], argvec[5], argvec[6], argvec[7], argvec[8]
  // args=pari_function_name, arg1, ...
  // or pari_function_name quoted to define a function

  enum {
    RET_GEN=0,
    RET_VOID=1,
    RET_INT=2,
    RET_LONG=3
  };
  typedef GEN (*PFGEN)(ANYARG);

  extern const unary_function_ptr * const  at_pari;
  static gen in_pari(const gen & args,GIAC_CONTEXT){
    if (args.type<_IDNT)
      return args;
    vecteur v(gen2vecteur(args));
    int vs=v.size();
    if (!vs){ // export all pari functions
      entree * ptr=functions_basic;
      gen tmp; int lextype;
      string redef;
      for (;ptr->name;++ptr){
	pari_function_table[ptr->name]=ptr;
	lextype=find_or_make_symbol(ptr->name,tmp,0,false,contextptr);
	if (lextype==T_SYMBOL)
	  sto(symbolic(at_pari,string2gen(ptr->name,false)),tmp,contextptr);
	else
	  redef += string(ptr->name) + " ";
	find_or_make_symbol(string("pari_")+ptr->name,tmp,0,false,contextptr);
	sto(symbolic(at_pari,string2gen(ptr->name,false)),tmp,contextptr);
      }
      return string2gen("All PARI functions are now defined with the pari_ prefix.\nPARI functions are also defined without prefix except:\n"+redef+"\nWhen working with p-adic numbers use them in a pari() call\nType ?pari for short help\nInside xcas, try Help->Manuals->PARI for HTML help",false);
    }
    if (v[0].is_symb_of_sommet(at_quote)){
      if (vs==1)
	return symbolic(at_pari,args);
      v[0]=v[0]._SYMBptr->feuille;
    }
    for (int i=1;i<vs;i++)
      v[i]=v[i].eval(eval_level(contextptr),contextptr);
    vecteur vars(1,identificateur("O"));
    lidnt(v,vars,false);
    vars.erase(vars.begin());
    bool parse_all=false;
    long av=get_pari_avma();
#ifdef PARI23
    if (setjmp(GP_DATA->env)){ 
      avma = av;
      parse_all=true;
      if (setjmp(GP_DATA->env)){ // if (setjmp(GP_DATA->env)){
#ifdef HAVE_LIBPTHREAD
	if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
#endif
	avma = av;
	*logptr(contextptr) << gettext("Error in PARI subsystem") << endl;
	// setsizeerr();
	return undef;
      } 
    }
#else
    cb_pari_err_recover=gp_err_recover;
    if (setjmp(env)){ 
      avma = av;
      parse_all=true;
      if (setjmp(env)){ // if (setjmp(GP_DATA->env)){
#ifdef HAVE_LIBPTHREAD
	if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
#endif
	avma = av;
	*logptr(contextptr) << gettext("Error in PARI subsystem") << endl;
	// setsizeerr();
	return undef;
      } 
    }
#endif
#ifdef EMCC
    parse_all=true; 
    // otherwise fibonacci or lngamma do not work, probably because of the 'L'
    // conversion (long argument expected) or precision (same)
#endif
    if (!parse_all && v[0].type==_STRNG) {
      string vstr=*v[0]._STRNGptr;
      if (vstr!="") {
	void * save_pari_stack_limit = PARI_stack_limit;
	PARI_stack_limit=0; // required since the stack changed
	if (vs==1)
	  return symbolic(at_pari,args);
	map<string,entree *>::const_iterator i = pari_function_table.find(vstr);
	// look at function prototype for return value
	// and call code, from anal.c line around 1990
	unsigned int ret;
	if (i!=pari_function_table.end()){
	  const char * s =i->second->code;
	  if      (*s <  'a')   ret = RET_GEN;
	  else if (*s == 'v') { ret = RET_VOID; s++; }
	  else if (*s == 'i') { ret = RET_INT;  s++; }
	  else if (*s == 'l') { ret = RET_LONG; s++; }
	  else                  ret = RET_GEN;
	  void * call= i->second->value;
	  // translate gen to GEN's 
	  GEN argvec[9]={0,0,0,0,0,0,0,0,0},res; long m;
	  int k=0;
	  for (int j=1;k<9 && *s && *s!='\n';++s){
	    switch(*s){
	    case 'L': // long
	      if (j==vs) {
		gen res=string2gen("PARI: Bad argument count",false);
		res.subtype=-1;
		return res;
	      }
	      argvec[k]= (GEN) v[j].val;
	      ++j; ++k;
	      break;
	    case 'P': // default precision
	      argvec[k] = (GEN) precdl; k++; break;
	    case 'D': //default param
	      {
		++s;
		switch(*s){
		case 'G': case '&': case 'I': case 'V': 
		  if (j<vs)
		    argvec[k]=ingen2GEN(v[j],vars,contextptr);
#ifdef PARI23
		  else
		    argvec[k]=utoi(0);
#endif
		  ++j; ++k; 
		  break;
		case 'n':
		  if (j<vs){
		    int pos=equalposcomp(vars,v[j]);
		    if (pos)
		      argvec[k]=(long int*)(pos -1);
		  }
		  else
		    argvec[k]=0;
		  ++j; ++k; 
		  break;
		default:
		  if (j<vs)
		    argvec[k]=(long int*) v[j].val;
		  else
		    argvec[k]=0;
		  ++k; ++j;
		  while (*s!= ',') s++;
		  s++;
		  while (*s!= ',') s++;
		}
		break;
	      }
	    case 'p':
	      argvec[k]=(GEN) precreal; 
	      ++k; 
	      break;
	    default:
	      if (j==vs) {
		gen res=string2gen("PARI: Bad argument count",false);
		res.subtype=-1;
		return res;
	      }
	      argvec[k]=ingen2GEN(v[j],vars,contextptr);
	      ++j; ++k;
	      break;
	    } 
	  }
	  switch (ret)
	    {
	    case RET_GEN:
	      res = ((PFGEN)call)(_ARGS_);
	      break;
	      
	    case RET_INT:
	      m = (long)((int (*)(ANYARG))call)(_ARGS_);
	      res = stoi(m); break;
	      
	    case RET_LONG:
	      m = ((long (*)(ANYARG))call)(_ARGS_);
	      res = stoi(m); break;
	      
	    case RET_VOID:
	      ((void (*)(ANYARG))call)(_ARGS_);
	      res = gnil; break;
	    }	  
	  // cerr << GEN2string(res) << endl;
	  gen resg(GEN2gen(res,vars));
	  PARI_stack_limit = save_pari_stack_limit ;
	  avma=av;
	  return resg;
	} // end if (i!=pari_function_table.end())
      } // end if vstr!=""
      if (vstr=="" && vs==2){
	long av=get_pari_avma();
	gen res= GEN2gen(gen2GEN(v[1],vars,contextptr),vars);
	avma=av;
	return res;
      }
    } // end if (!parse_all ...)
    string s;
    if (v[0].type==_FUNC)
      s=v[0]._FUNCptr->ptr()->s;
    else
      s=gen2string(v[0]);
    if (vs>1){
      s+="(";
      for (int i=1;i<vs;){
	s += pariprint(v[i],0,contextptr);
	++i;
	if (i==vs)
	  break;
	s += ",";
      }
      s +=")";
    }
    return pari_exec(s,contextptr);
  }
  gen _pari(const gen & args,GIAC_CONTEXT){
    if ( args.type==_STRNG && args.subtype==-1) return  args;
    abort_if_locked();
    gen res;
#ifdef HAVE_LIBPTHREAD
    pthread_cleanup_push(pari_cleanup, (void *) pari_mutex_ptr);
#endif
    res=in_pari(args,contextptr);
#ifdef HAVE_LIBPTHREAD
    if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
    pthread_cleanup_pop(0);
#endif
    return res;
  }
  static const char _pari_s []="pari";
  static define_unary_function_eval (__pari,&giac::_pari,_pari_s);
  define_unary_function_ptr5( at_pari ,alias_at_pari,&__pari,_QUOTE_ARGUMENTS,true);

  gen _pari_unlock(const gen & args,GIAC_CONTEXT){
    if ( args.type==_STRNG && args.subtype==-1) return  args;
    int locked=check_pari_mutex();
    if (!locked)
      return 0;
#ifdef HAVE_LIBPTHREAD
    delete pari_mutex_ptr;
    pari_mutex_ptr = 0;
#endif
    return 1;
  }
  static const char _pari_unlock_s []="pari_unlock";
  static define_unary_function_eval (__pari_unlock,&giac::_pari_unlock,_pari_unlock_s);
  define_unary_function_ptr5( at_pari_unlock ,alias_at_pari_unlock,&__pari_unlock,_QUOTE_ARGUMENTS,true);

  static std::string cutstring(const std::string & s,int ncol){
    string res;
    if (ncol<20)
      ncol=20;
    int left=s.size(),pos=0,j;
    for (;left>ncol;pos+=j,left-=j){
      for (j=ncol;j>ncol/2;--j){
	if (s[pos+j]==' ')
	  break;
      }
      res=res+s.substr(pos,j)+'\n';
    }
    return res+s.substr(pos,pos+left);
  }

  // Help for g calling PARI online help
  std::string pari_help(const gen & g){
    if (is_zero(g))
      return "Run pari() to export PARI functions.\n?pari(1) to ?pari(11) lists PARI functions by section\n?pari_functionname shows a short help on a function\nInside Xcas, Help->Manual->PARI-GP shows HTML help";
    string res;
    if (g.type==_INT_){
      int section=g.val;
      entree * ptr=functions_basic;
      for (;ptr->name;++ptr){
	if (ptr->menu==section){
	  res += ptr->name;
	  res += " ";
	}
      }
      return cutstring(res,70);
    }
    string gs;
    if (g.type==_FUNC)
      gs=g._FUNCptr->ptr()->s;
    else
      gs=gen2string(g);
    if (gs.size()>5 && gs.substr(0,5)=="pari_")
      gs=gs.substr(5,gs.size()-5);
    entree * ptr=functions_basic;
    for (;ptr->name;++ptr){
      if (ptr->name==gs){
	res = ptr->help;
	return cutstring(res,70);
      }
    }
    return "PARI function not found\nHelp syntax: ?pari(1),...,?pari(12) or ?pari_functionname";
  }

  bool pari_polroots(const vecteur & p,vecteur & res,long prec,GIAC_CONTEXT){
    if (check_pari_mutex())
      return false;
    gen tmp;
#ifdef HAVE_LIBPTHREAD
    pthread_cleanup_push(pari_cleanup, (void *) pari_mutex_ptr);
#endif
    long av=get_pari_avma();
    GEN G=gen2GEN(change_subtype(p,_POLY1__VECT),vecteur(0),contextptr);
    if (debug_infolevel)
      CERR << "pari_polroots " << GEN2gen(G,vecteur(1,vx_var)) << endl;
    G=roots(G,prec);
    tmp=GEN2gen(G,vecteur(0));
    avma=av;
#ifdef HAVE_LIBPTHREAD
    if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
    pthread_cleanup_pop(0);
#endif
    if (tmp.type!=_VECT)
      return false;
    res=*tmp._VECTptr;
    return true;
  }

  bool pari_polresultant(const gen & p,const gen & q,const vecteur & lv,gen & res,GIAC_CONTEXT){
    if (check_pari_mutex())
      return false;
    gen tmp;
#ifdef HAVE_LIBPTHREAD
    pthread_cleanup_push(pari_cleanup, (void *) pari_mutex_ptr);
#endif
    long av=get_pari_avma();
    void * save_pari_stack_limit = PARI_stack_limit;
    PARI_stack_limit=0; 
    GEN P=gen2GEN(p,lv,contextptr);
    GEN Q=gen2GEN(q,lv,contextptr);
    GEN PQ=polresultant0(P,Q,-1,2);
    tmp=GEN2gen(PQ,lv);
    avma=av;
    PARI_stack_limit=save_pari_stack_limit;
#ifdef HAVE_LIBPTHREAD
    if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
    pthread_cleanup_pop(0);
#endif
    res=tmp;
    return true;
  }

  bool pari_nffactor(const gen & p,const gen & pmin,const vecteur & lv,gen & res,GIAC_CONTEXT){
    if (check_pari_mutex())
      return false;
    gen tmp;
#ifdef HAVE_LIBPTHREAD
    pthread_cleanup_push(pari_cleanup, (void *) pari_mutex_ptr);
#endif
    long av=get_pari_avma();
    void * save_pari_stack_limit = PARI_stack_limit;
    PARI_stack_limit=0; 
    GEN P=gen2GEN(p,lv,contextptr);
    GEN Pmin=gen2GEN(pmin,lv,contextptr);
    int prec=decimal_digits(contextptr);
    if (prec<30)
      prec=30;
#if 0
    GEN nf=nfinit0(Pmin,0,prec);
    tmp=GEN2gen(nffactor(nf,P),lv);
#else
    tmp=GEN2gen(nffactor(Pmin,P),lv);
#endif
    avma=av;
    PARI_stack_limit=save_pari_stack_limit;
#ifdef HAVE_LIBPTHREAD
    if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
    pthread_cleanup_pop(0);
#endif
    res=tmp;
    return true;
  }

  bool pari_galoisconj(const gen & g,vecteur & w,GIAC_CONTEXT){
    if (check_pari_mutex())
      return false;
    gen res;
#ifdef HAVE_LIBPTHREAD
    pthread_cleanup_push(pari_cleanup, (void *) pari_mutex_ptr);
#endif
    res=in_pari(makesequence(string2gen("nfgaloisconj",false),g),contextptr);
#ifdef HAVE_LIBPTHREAD
    if (pari_mutex_ptr) pthread_mutex_unlock(pari_mutex_ptr);    
    pthread_cleanup_pop(0);
#endif
    if (res.type!=_VECT)
      return false;
    w=*res._VECTptr;
    gen gp=_symb2poly(makesequence(g,vx_var),contextptr);
    for (int i=0;i<int(w.size());++i){
      gen tmp=w[i];
      tmp=_symb2poly(makesequence(tmp,vx_var),contextptr);
      gen d=1;
      if (tmp.type==_VECT)
	lcmdeno(*tmp._VECTptr,d,contextptr);
      w[i]=symb_rootof(tmp,gp,contextptr)/d;
    }
    return true;
  }
  
#ifndef NO_NAMESPACE_GIAC
}
#endif // ndef NO_NAMESPACE_GIAC

#else
// ! HAVE_LIBPARI
#ifndef NO_NAMESPACE_GIAC
namespace giac {
#endif // ndef NO_NAMESPACE_GIAC

  static gen pari_error(){
    return undeferr(gettext("Not implemented, please recompile giac with PARI"));
  }

  gen pari_isprime(const gen & e,int certif){
    return string2gen("please recompile giac with PARI",false);
  }
  std::string pari_ifactor(const gen & e){
    return "please recompile giac with PARI";
  }
  
  gen pari_ffinit(const gen & p,int n){
    return string2gen("please recompile giac with PARI",false);
  }

  gen pari_gamma(const gen & e){
    return pari_error();
  }
  gen pari_zeta(const gen & e){
    return pari_error();
  }
  gen pari_psi(const gen & e){
    return pari_error();
  }
  
  bool pari_lift_combine(const vecteur& a, const std::vector<vecteur>& factmod,
			 gen& modulo, std::vector<vecteur>& res){
    vecteur tmp(1,pari_error());
    res=std::vector<vecteur>(1,tmp);
    return false;
  }

  bool pari_galoisconj(const gen & g,vecteur & w,GIAC_CONTEXT){
    return false;
  }
  
  gen _pari(const gen & args,GIAC_CONTEXT){
    if ( args.type==_STRNG && args.subtype==-1) return  args;
    return pari_error();
  }
  
  std::string pari_help(const gen & g){
    return "please recompile giac with PARI";
  }

  gen GEN2gen(const GEN & G,const vecteur & vars){
    return gensizeerr("please recompile giac with PARI");
  }
  GEN gen2GEN(const gen & e,const vecteur & vars,GIAC_CONTEXT){
    return 0;
  }

  bool pari_polroots(const vecteur & p,vecteur & res,long prec,GIAC_CONTEXT){
    return false;
  }
  bool pari_polresultant(const gen & p,const gen & q,const vecteur & lv,gen & res,GIAC_CONTEXT){
    return false;
  }
  bool pari_nffactor(const gen & p,const gen & pmin,const vecteur & lv,gen & res,GIAC_CONTEXT){
    return false;
  }
  static const char _pari_s []="pari";
  static define_unary_function_eval (__pari,&giac::_pari,_pari_s);
  define_unary_function_ptr5( at_pari ,alias_at_pari,&__pari,_QUOTE_ARGUMENTS,true);

  static const char _pari_unlock_s []="pari_unlock";
  static define_unary_function_eval (__pari_unlock,&giac::_pari,_pari_unlock_s);
  define_unary_function_ptr5( at_pari_unlock ,alias_at_pari_unlock,&__pari_unlock,_QUOTE_ARGUMENTS,true);


#ifndef NO_NAMESPACE_GIAC
}
#endif // ndef NO_NAMESPACE_GIAC

#endif // HAVE_LIBPARI