size()==2){
gen P=g._VECTptr->front();
gen x=g._VECTptr->back();
gen Px=_e2r(makevecteur(P,x),contextptr);
if (Px.type==_FRAC)
Px=inv(Px._FRACptr->den,contextptr)*Px._FRACptr->num;
if (Px.type!=_VECT)
return gensizeerr();
v=*Px._VECTptr;
}
else
v=*g._VECTptr;
return companion(v);
}
static const char _companion_s []="companion";
static define_unary_function_eval (__companion,&_companion,_companion_s);
define_unary_function_ptr5( at_companion ,alias_at_companion,&__companion,0,true);
gen _border(const gen & g,GIAC_CONTEXT){
if ( g.type==_STRNG && g.subtype==-1) return g;
if (g.type!=_VECT || g._VECTptr->size()!=2 || !ckmatrix(g._VECTptr->front()) || g._VECTptr->back().type!=_VECT)
return gensizeerr();
matrice m=*g._VECTptr->front()._VECTptr,v=*g._VECTptr->back()._VECTptr;
if (m.size()!=v.size())
return gendimerr();
m=mtran(m);
if (ckmatrix(v))
m=mergevecteur(m,mtran(v));
else
m.push_back(v);
return mtran(m);
}
static const char _border_s []="border";
static define_unary_function_eval (__border,&_border,_border_s);
define_unary_function_ptr5( at_border ,alias_at_border,&__border,0,true);
// pade(f(x),x,n,p) or pade(f(x),x,N(x),p)
// Find P/Q such that P/Q=f(x) mod N(x) [1st case N(x)=x^n], deg(P)size()!=4)
return gensizeerr();
vecteur v =*g._VECTptr;
if (v[1].type!=_IDNT || v[3].type!=_INT_)
return gensizeerr();
int n,p=v[3].val,ns;
gen f=v[0],x=v[1],N=v[2];
if (v[2].type==_INT_){
n=v[2].val+1;
N=pow(v[1],n);
ns=n;
}
else {
n=p;
ns=_degree(makevecteur(N,x),contextptr).val;
}
if (p<=0 || n<=0)
return gensizeerr();
f=_series(makevecteur(f,x,zero,ns-1),contextptr);
vecteur l1(lop(f,at_order_size));
vecteur l2(l1.size(),zero);
f=subst(f,l1,l2,false,contextptr);
// Convert f and N to poly1
vecteur l(1,x);
lvar(f,l);
lvar(N,l);
int ls=int(l.size());
gen ff(sym2r(f,l,contextptr));
gen fn,fd;
fxnd(ff,fn,fd);
gen Nf(sym2r(N,l,contextptr));
gen Nn,Nd;
fxnd(Nf,Nn,Nd); // Note nd has no meaning
vecteur fp;
if (fn.type==_POLY)
fp=polynome2poly1(*fn._POLYptr,1);
else
return gensizeerr();
vecteur Np;
if (Nn.type==_POLY)
Np=polynome2poly1(*Nn._POLYptr,1);
else
return gensizeerr();
n=int(Np.size())-1;
// Check that fp is a poly of degree less than n
if (n<1 || p>n || signed(fp.size())>n)
return gendimerr();
vecteur a,b;
if (!egcd_pade(Np,fp,p,a,b,0))
*logptr(contextptr) << gettext("Solution may be wrong since a and b are not prime together: ")+gen(a).print(contextptr)+","+gen(b).print(contextptr) << endl;
gen res=poly12polynome(a,1,ls);
res=res/(fd*gen(poly12polynome(b,1,ls)));
res=r2sym(res,l,contextptr);
return res;
}
static const char _pade_s []="pade";
static define_unary_function_eval (__pade,&_pade,_pade_s);
define_unary_function_ptr5( at_pade ,alias_at_pade,&__pade,0,true);
static const char _interp_s []="interp";
static define_unary_function_eval (__interp,&_lagrange,_interp_s);
define_unary_function_ptr5( at_interp ,alias_at_interp,&__interp,0,true);
static gen lhsrhs(const gen & g,int i){
if (!is_equal(g) && !g.is_symb_of_sommet(at_interval))
return gensizeerr();
gen & f=g._SYMBptr->feuille;
if (f.type!=_VECT || f._VECTptr->size()!=2)
return gensizeerr();
return (*f._VECTptr)[i];
}
gen _lhs(const gen & g,GIAC_CONTEXT){
if ( g.type==_STRNG && g.subtype==-1) return g;
return lhsrhs(g,0);
}
static const char _lhs_s []="lhs";
static define_unary_function_eval (__lhs,&_lhs,_lhs_s);
define_unary_function_ptr5( at_lhs ,alias_at_lhs,&__lhs,0,true);
gen _rhs(const gen & g,GIAC_CONTEXT){
if ( g.type==_STRNG && g.subtype==-1) return g;
return lhsrhs(g,1);
}
static const char _rhs_s []="rhs";
static define_unary_function_eval (__rhs,&_rhs,_rhs_s);
define_unary_function_ptr5( at_rhs ,alias_at_rhs,&__rhs,0,true);
// Given [v_0 ... v_(2n-1)] (begin of the recurrence sequence)
// return [b_n...b_0] such that b_n*v_{n+k}+...+b_0*v_k=0
// Example [1,-1,3,3] -> [1,-3,-6]
// -> the recurrence relation is v_{n+2}=3v_{n+1}+6v_n
gen _reverse_rsolve(const gen & g,GIAC_CONTEXT){
if ( g.type==_STRNG && g.subtype==-1) return g;
if (g.type!=_VECT)
return gensizeerr();
vecteur v=reverse_rsolve(*g._VECTptr);
return v/v.front();
}
static const char _reverse_rsolve_s []="reverse_rsolve";
static define_unary_function_eval (__reverse_rsolve,&_reverse_rsolve,_reverse_rsolve_s);
define_unary_function_ptr5( at_reverse_rsolve ,alias_at_reverse_rsolve,&__reverse_rsolve,0,true);
// Approx fft or exact if args=poly1,omega,n
gen fft(const gen & g_orig,int direct,GIAC_CONTEXT){
if (g_orig.type==_VECT && g_orig.subtype==_SEQ__VECT && g_orig._VECTptr->size()==3 && g_orig._VECTptr->front().type==_VECT){
vecteur & v =*g_orig._VECTptr->front()._VECTptr;
int n=int(v.size());
if (n<2)
return gendimerr();
gen omega=(*g_orig._VECTptr)[1];
gen modulo=(*g_orig._VECTptr)[2];
if (direct==1)
omega=invmod(omega,modulo);
if (omega.type==_INT_ && modulo.type==_INT_ && n==(1<<(sizeinbase2(n)-1))){
if (debug_infolevel)
CERR << CLOCK()*1e-6 << " fft start" << endl;
vector A; A.reserve(n);
for (int i=0;imoduloon=true;
env->modulo=modulo;
fft(v,w,res,env);
delete env;
if (direct==1)
return res;
else
return smod(invmod(n,modulo)*res,modulo);
}
gen g=g_orig;
if (g.type!=_VECT)
return gensizeerr();
int n=int(g._VECTptr->size());
if (n<2)
return gendimerr();
vector< complex > vd;
if (convert(*g._VECTptr,vd,true)){
if (debug_infolevel)
CERR << CLOCK()*1e-6 << " fft start" << endl;
double theta=2.0*M_PI/n;
if (direct) theta=-theta;
bool done=false;
if (n==(1<<(sizeinbase2(n)-1))){
fft2(&vd.front(),n,theta);
done=true;
}
#if 1 //ndef HAVE_LIBGSL
if (!done){
complex w(std::cos(theta),std::sin(theta));
vector< complex > W(n),res(n);
for (int i=0;i(std::cos(i*theta),std::sin(i*theta));
else
W[i]=W[i-1]*w;
}
fft(&vd[0],n,&W[0],n,&res[0]);
done=true;
}
#endif
if (done){
gen r=vecteur(0);
vecteur & v=*r._VECTptr;
v.reserve(n);
if (direct){
for (int i=0;isize()==2){
return sto(g._VECTptr->back(),g._VECTptr->front(),contextptr);
}
if (is_equal(g)){
gen & f=g._SYMBptr->feuille;
if (f.type==_VECT && f._VECTptr->size()==2)
return sto(f._VECTptr->back(),f._VECTptr->front(),contextptr);
}
if (g.type!=_VECT)
return gensizeerr();
return apply(g,_assign,contextptr);
}
static const char _assign_s []="assign";
static define_unary_function_eval (__assign,&_assign,_assign_s);
define_unary_function_ptr5( at_assign ,alias_at_assign,&__assign,0,true);
gen _implicitplot3d(const gen & g,GIAC_CONTEXT){
if ( g.type==_STRNG && g.subtype==-1) return g;
return _plotimplicit(g,contextptr);
}
static const char _implicitplot3d_s []="implicitplot3d";
static define_unary_function_eval (__implicitplot3d,&_implicitplot3d,_implicitplot3d_s);
define_unary_function_ptr5( at_implicitplot3d ,alias_at_implicitplot3d,&__implicitplot3d,0,true);
void convert_double_int(vecteur & v){
for (unsigned i=0;i1 && v[1].type!=_VECT)
v=makevecteur(1,v);
gen g=v[0];
if (ckmatrix(g)){
if (v.size()==2 && v[1].type==_INT_){
sample_rate=v[1].val;
v=*g._VECTptr;
v.insert(v.begin(),makevecteur(1,16,sample_rate));
g=v[0];
}
else {
if (v.size()==3 && v[1].type==_INT_ && v[2].type==_INT_){
sample_rate=v[1].val;
bits_per_sample=v[2].val;
v=*g._VECTptr;
v.insert(v.begin(),makevecteur(1,bits_per_sample,sample_rate));
g=v[0];
}
}
}
// g=channels,bits_per_sample,sample_rate,data_size,
if (g.type==_INT_)
g=makevecteur(g,16);
if (g.type!=_VECT || g._VECTptr->empty())
return false;
vecteur w=*g._VECTptr;
if (w.size()==1)
w.push_back(16);
int ws=int(w.size());
if (w[0].type!=_INT_ || w[1].type!=_INT_)
return false;
channels=w[0].val;
if (channels<=0 || channels>int(v.size())){
channels=int(v.size());
if (v.front().type!=_VECT || !is_integer_vecteur(*v.front()._VECTptr))
return false;
data_size=unsigned(v.front()._VECTptr->size());
for (int i=1;iv[i]._VECTptr->size())
data_size=unsigned(v[i]._VECTptr->size());
}
w=makevecteur(channels,16,44100); ws=3;
g=w;
v.insert(v.begin(),g);
}
if (channels<=0 || channels>4 || int(v.size())<=channels)
return false;
bits_per_sample=(w[1].val/8)*8;
if (bits_per_sample<=0)
return false;
if (ws>=3)
sample_rate=w[2].val;
else
sample_rate=44100;
if (sample_rate<1)
return false;
if (ws>=4 && w[3].type==_INT_)
data_size=w[3].val;
else
data_size=RAND_MAX;
for (int i=1;i<=channels;++i){
if (v[i].type!=_VECT || !is_integer_vecteur(*v[i]._VECTptr))
return false;
if (data_size>v[i]._VECTptr->size())
data_size=unsigned(v[i]._VECTptr->size());
}
return true;
}
#ifdef HAVE_LIBAO
#include
typedef unsigned short aou16;
typedef unsigned int aou32;
gen _playsnd(const gen & args,GIAC_CONTEXT){
if (args.type==_STRNG){
if (args.subtype==-1) return args;
return _playsnd(_readwav(args,contextptr),contextptr);
}
ao_device *device=0;
ao_sample_format format;
int default_driver;
ao_initialize();
default_driver = ao_default_driver_id();
memset(&format, 0, sizeof(format));
format.bits = 16;
format.channels = 2;
format.rate = 44100;
format.byte_format = AO_FMT_LITTLE;
unsigned int data_size=0;
vecteur v;
if (args.type==_VECT && !args._VECTptr->empty()){
// set format
v=*args._VECTptr;
if (!read_audio(v,format.channels,format.rate,format.bits,data_size))
return gensizeerr(gettext("Invalid sound data"));
}
if (data_size){
*logptr(contextptr) << gettext("Using sound parameters: channels, rate, bits, records ") << format.channels << "," << format.rate << "," << format.bits << "," << data_size << endl;
device = ao_open_live(default_driver, &format, NULL /* no options */);
if (device == NULL)
return gensizeerr(gettext("Error opening audio device."));
unsigned n=data_size*format.channels*format.bits/8;
char * buffer=(char *)malloc(n*sizeof(char));
// aou16 * bufshort=(aou16 *) buffer;
aou32 * bufint=(aou32 *) buffer;
if (buffer){
// copy data from v into buffer and play it
unsigned c=format.channels,b=format.bits/8;
for (unsigned i=0;i>8) & 0xff;
}
if (b==4)
bufint[i*c+j]=u;
}
}
ao_play(device, buffer, n);
}
}
ao_close(device);
ao_shutdown();
return 1;
}
static const char _playsnd_s []="playsnd";
static define_unary_function_eval (__playsnd,&_playsnd,_playsnd_s);
define_unary_function_ptr5( at_playsnd ,alias_at_playsnd,&__playsnd,0,true);
#else
#if !defined GIAC_GGB && defined EMCC // must have EM_ASM code javascript inlined (emscripten 1.30.4 at least?)
#include
gen _playsnd(const gen & args,GIAC_CONTEXT){
if (args.type==_STRNG){
if (args.subtype==-1) return args;
return _playsnd(_readwav(args,contextptr),contextptr);
}
int nbits = 16;
int nchannels = 2;
int nrate = 44100;
unsigned int data_size=0;
vecteur v;
if (args.type==_VECT && !args._VECTptr->empty()){
// set format
v=*args._VECTptr;
if (!read_audio(v,nchannels,nrate,nbits,data_size))
return gensizeerr(gettext("Invalid sound data"));
}
if (data_size){
*logptr(contextptr) << gettext("Using sound parameters: channels, rate, bits, records ") << nchannels << "," << nrate << "," << data_size << endl;
unsigned nDataBytes=data_size*nchannels*sizeof(float);
// copy data from v into buffer and play it
unsigned b=nbits/8;
float * ptr = (float *) malloc(nDataBytes);
for (unsigned j=0;jsize()==2){
n=args._VECTptr->front();
rate=args._VECTptr->back();
}
else
n=args;
n=evalf_double(n,1,contextptr);
if (n.type!=_DOUBLE_ || n._DOUBLE_val<=0 || rate.type!=_INT_ || rate.val < 1 )
return gensizeerr(gettext("Invalid sound parameters"));
double r=evalf_double(rate,1,contextptr)._DOUBLE_val;
double nr=r*n._DOUBLE_val;
if (nr>LIST_SIZE_LIMIT)
return gensizeerr("Too many records");
vecteur v;
v.reserve(int(nr));
for (int i=0;ireserve(n/(channels*bits_per_sample));
}
while (n>0 && !feof(f)){
for (int i=1;i<=channels;++i){
u=0;
if (fread(&u,bits_per_sample,1,f)!=1)
return false;
// for (int j=0;jpush_back(int(u));
if (n<=0)
break;
}
}
return true;
}
static bool readwav(const string & s,gen & g){
FILE * f = fopen(s.c_str(),"r");
if (!f)
return false;
bool res=in_readwav(f,g);
fclose(f);
return res;
}
gen _readwav(const gen & g,GIAC_CONTEXT){
if ( g.type==_STRNG && g.subtype==-1) return g;
if (g.type!=_STRNG)
return gensizeerr();
gen res;
bool ok= readwav(*g._STRNGptr,res);
if (!ok)
return gensizeerr(gettext("File not found or unrecognized wav file format"));
return res;
}
static const char _readwav_s []="readwav";
static define_unary_function_eval (__readwav,&_readwav,_readwav_s);
define_unary_function_ptr5( at_readwav ,alias_at_readwav,&__readwav,0,true);
static bool in_writewav(FILE * f,const vecteur & v_){
if (v_.empty())
return false;
vecteur v(v_);
int channels,sample_rate=44100,bits_per_sample=0;
unsigned int u,byte_rate,block_align=0,data_size=1U<<31;
if (!read_audio(v,channels,sample_rate,bits_per_sample,data_size))
return false;
u=0x46464952;
if (fwrite(&u,4,1,f)!=1)
return false;
u=36+data_size*bits_per_sample/8*channels;
if (fwrite(&u,4,1,f)!=1)
return false;
u=0x45564157;
if (fwrite(&u,4,1,f)!=1)
return false;
u=0x20746d66;
if (fwrite(&u,4,1,f)!=1)
return false;
u=0x10;
if (fwrite(&u,4,1,f)!=1)
return false;
fputc(1,f);
fputc(0,f);
fputc(channels,f);
fputc(0,f);
if (fwrite(&sample_rate,4,1,f)!=1)
return false;
byte_rate = sample_rate * channels * bits_per_sample/8;
if (fwrite(&byte_rate,4,1,f)!=1)
return false;
block_align=channels * bits_per_sample/8;
if (fwrite(&block_align,2,1,f)!=1)
return false;
if (fwrite(&bits_per_sample,2,1,f)!=1)
return false;
u=0x61746164; // "data"
if (fwrite(&u,4,1,f)!=1)
return false;
u= bits_per_sample/8* data_size*channels; // should be data_size
if (fwrite(&u,4,1,f)!=1)
return false;
// write data
u /= channels;
bits_per_sample /= 8;
unsigned n = u/ bits_per_sample;
for (unsigned i=0;isize()!=2 || g._VECTptr->front().type!=_STRNG || g._VECTptr->back().type!=_VECT)
return gensizeerr();
bool ok= writewav(*g._VECTptr->front()._STRNGptr,*g._VECTptr->back()._VECTptr);
if (!ok)
return gensizeerr(gettext("Unable to open file or unable to code data"));
return 1;
}
static const char _writewav_s []="writewav";
static define_unary_function_eval (__writewav,&_writewav,_writewav_s);
define_unary_function_ptr5( at_writewav ,alias_at_writewav,&__writewav,0,true);
#endif // RTOS_THREADX
static gen animate2d3d(const gen & g,bool dim3,GIAC_CONTEXT){
int s=0,frames=10;
if (g.type!=_VECT || (s=int(g._VECTptr->size()))<3)
return gensizeerr();
vecteur v = *g._VECTptr;
gen t;
double tmin,tmax;
if (!readrange(v[2],gnuplot_tmin,gnuplot_tmax,t,tmin,tmax,contextptr))
return gensizeerr();
// find a frame argument, otherwise use 10
for (int i=3;ifeuille;
if (f.type==_VECT && f._VECTptr->size()==2 && f._VECTptr->front().type==_INT_ && f._VECTptr->front().val==_FRAMES && f._VECTptr->back().type==_INT_)
frames=f._VECTptr->back().val;
}
}
if (frames==0)
return gensizeerr();
v.erase(v.begin()+2);
gen h=symbolic(dim3?at_plotfunc:at_plot,gen(v,g.subtype));
v=makevecteur(h,t,tmin,tmax,(tmax-tmin)/frames);
gen tmpseq=seqprod(v,0,contextptr);
if (is_undef(tmpseq)) return tmpseq;
return symbolic(at_animation,tmpseq);
}
gen _animate(const gen & g,GIAC_CONTEXT){
if ( g.type==_STRNG && g.subtype==-1) return g;
return animate2d3d(g,false,contextptr);
}
static const char _animate_s []="animate";
static define_unary_function_eval (__animate,&_animate,_animate_s);
define_unary_function_ptr5( at_animate ,alias_at_animate,&__animate,0,true);
gen _animate3d(const gen & g,GIAC_CONTEXT){
if ( g.type==_STRNG && g.subtype==-1) return g;
return animate2d3d(g,true,contextptr);
}
static const char _animate3d_s []="animate3d";
static define_unary_function_eval (__animate3d,&_animate3d,_animate3d_s);
define_unary_function_ptr5( at_animate3d ,alias_at_animate3d,&__animate3d,0,true);
gen _even(const gen & g_,GIAC_CONTEXT){
gen g(g_);
if ( g.type==_STRNG && g.subtype==-1) return g;
if (!is_integral(g)) return gentypeerr(contextptr);
return is_zero(smod(g,2));
}
static const char _even_s []="even";
static define_unary_function_eval (__even,&_even,_even_s);
define_unary_function_ptr5( at_even ,alias_at_even,&__even,0,true);
gen _odd(const gen & g_,GIAC_CONTEXT){
gen g(g_);
if ( g.type==_STRNG && g.subtype==-1) return g;
if (!is_integral(g)) return gentypeerr(contextptr);
return !is_zero(smod(g,2));
}
static const char _odd_s []="odd";
static define_unary_function_eval (__odd,&_odd,_odd_s);
define_unary_function_ptr5( at_odd ,alias_at_odd,&__odd,0,true);
#ifdef HAVE_LIBPNG
int write_png(const char *file_name, void *rows_, int w, int h, int colortype, int bitdepth){
png_bytep * rows=(png_bytep *) rows_;
png_structp png_ptr;
png_infop info_ptr;
FILE *fp = fopen(file_name, "wb");
const char *doing = "open for writing";
if (!(fp = fopen(file_name, "wb"))) goto fail;
doing = "create png write struct";
if (!(png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL,
NULL, NULL))) goto fail;
doing = "create png info struct";
if (!(info_ptr = png_create_info_struct(png_ptr))) goto fail;
if (setjmp(png_jmpbuf(png_ptr))) goto fail;
doing = "init IO";
png_init_io(png_ptr, fp);
doing = "write header";
png_set_IHDR(png_ptr, info_ptr, w, h, bitdepth, colortype,
PNG_INTERLACE_NONE,
PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
doing = "write info";
png_write_info(png_ptr, info_ptr);
doing = "write image";
png_write_image(png_ptr, rows);
doing = "write end";
png_write_end(png_ptr, NULL);
fclose(fp);
return 0;
fail:
printf("Write_png: could not %s\n", doing);
return -1;
}
static bool writergb(const string & filename,const vecteur & v){
if (v.size()==2 && ckmatrix(v[1])){
int w,h;
mdims(*v[1]._VECTptr,h,w);
unsigned i=0,taille=w*h;
unsigned char * screenbuf=new unsigned char[w*h];
unsigned rowbytes = w;
// fill screenbuf with data from v[] order 1,2,4,3
for (;i2 && v[1].type==_INT_) ?v[1].val:0),
( (s>2 && v[2].type==_INT_) ?v[2].val:0),
res);
if (!ok)
return gensizeerr(gettext("File not found or unrecognized image file format"));
return res;
}
static const char _readrgb_s []="readrgb";
static define_unary_function_eval (__readrgb,&_readrgb,_readrgb_s);
define_unary_function_ptr5( at_readrgb ,alias_at_readrgb,&__readrgb,0,true);
gen linear_apply(const gen & e,const gen & x,const gen & l,gen & remains, GIAC_CONTEXT, gen (* f)(const gen &,const gen &,const gen &,gen &,const context *)){
if (is_constant_wrt(e,x,contextptr) || (e==x) )
return f(e,x,l,remains,contextptr);
// e must be of type _SYMB
if (e.type!=_SYMB) return gensizeerr(gettext("in linear_apply"));
unary_function_ptr u(e._SYMBptr->sommet);
gen arg(e._SYMBptr->feuille);
gen res;
if (u==at_neg){
res=-linear_apply(arg,x,l,remains,contextptr,f);
remains=-remains;
return res;
} // end at_neg
if (u==at_plus){
if (arg.type!=_VECT)
return linear_apply(arg,x,l,remains,contextptr,f);
const_iterateur it=arg._VECTptr->begin(),itend=arg._VECTptr->end();
for (gen tmp;it!=itend;++it){
res = res + linear_apply(*it,x,l,tmp,contextptr,f);
remains =remains + tmp;
}
return res;
} // end at_plus
if (u==at_prod){
if (arg.type!=_VECT)
return linear_apply(arg,x,l,remains,contextptr,f);
// find all constant terms in the product
vecteur non_constant;
gen prod_constant;
decompose_prod(*arg._VECTptr,x,non_constant,prod_constant,true,contextptr);
if (non_constant.empty()) return gensizeerr(gettext("in linear_apply 2")); // otherwise the product would be constant
if (non_constant.size()==1)
res = linear_apply(non_constant.front(),x,l,remains,contextptr,f);
else
res = f(symbolic(at_prod,non_constant),x,l,remains,contextptr);
remains = prod_constant * remains;
return prod_constant * res;
} // end at_prod
return f(e,x,l,remains,contextptr);
}
// discrete product primitive of a polynomial P
gen product(const polynome & P,const vecteur & v,const gen & n,gen & remains,GIAC_CONTEXT){
// factor P, for a linear factor a*n+b, multiply res by a^n*gamma(n+b/a)
// for other factors multiply remains by the factor
polynome Pcont;
factorization f;
gen divan=1,res,extra_div=1;
if (!factor(P,Pcont,f,/* is_primitive*/false,/* with_sqrt*/true,/* complex */true,divan,extra_div) || extra_div!=1){
remains=r2e(P,v,contextptr);
return 1;
}
res = pow(divan,-n,contextptr);
factorization::const_iterator it=f.begin(),itend=f.end();
for (;it!=itend;++it){
gen tmp=r2e(it->fact,v,contextptr);
if (it->fact.lexsorted_degree()!=1){
remains = remains * pow(tmp,it->mult);
}
else {
gen a=derive(tmp,n,contextptr);
if (is_undef(a))
return a;
gen b=normal(tmp-a*n,contextptr);
res = res * pow(a,it->mult*n,contextptr) * pow(symbolic(at_factorial,n+b/a-1),it->mult,contextptr);
}
}
return res*pow(r2e(Pcont,v,contextptr),n,contextptr);
}
// product(P,n=a..b) where the first variable in v is n
gen product(const polynome & P,const vecteur & v,const gen & n,const gen & a,const gen & b,GIAC_CONTEXT){
gen remains(1),res=product(P,v,n,remains,contextptr);
res=subst(res,n,b+1,false,contextptr)/subst(res,n,a,false,contextptr);
if (is_one(remains))
return res;
else
return res*symbolic(at_product,gen(makevecteur(remains,n,a,b),_SEQ__VECT));
}
// solve vnext*sol[n+1]+v*sol[n]+vcst=0
// example: u_{n+1}=u_n/n+1 -> n*u_{n+1}-u_n-n=0
// vnext=[1,0], v=[1], vcst=[-1,0] -> n solution
static bool rsolve_particular(const modpoly & vnext,const modpoly & v,const modpoly & vcst,modpoly & sol,GIAC_CONTEXT){
// find majoration of degree for solution
int vnextdeg=int(vnext.size())-1,vdeg=int(v.size())-1,vcstdeg=int(vcst.size())-1,soldeg;
soldeg=vcstdeg-giacmax(vnextdeg,vdeg);
if (vnextdeg==vdeg && is_zero(vnext.front()+v.front()))
soldeg++;
if (soldeg<0)
return false;
modpoly vars(soldeg+1);
for (int i=0;i<=soldeg;++i){
vars[i]=identificateur("rsolve_x"+print_INT_(i));
}
modpoly equations=operator_plus(operator_times(vnext,taylor(vars,1,0),0),
operator_times(v,vars,0),0);
equations=operator_plus(equations,vcst,0);
sol=linsolve(equations,vars,contextptr);
return !sol.empty() && !is_undef(sol);
}
// solve a recurrence relation x_{n+1}=l0*x_n+e
static gen rsolve(const gen & e,const gen & n,const gen & l0,gen & remains,GIAC_CONTEXT){
if (is_zero(e)){
remains=0;
return e;
}
vecteur v;
polynome P,Q,R;
if (!is_hypergeometric(e,*n._IDNTptr,v,P,Q,R,contextptr)){
*logptr(contextptr) << gettext("Cst part must be hypergeometric") << endl;
remains=e;
return 0;
}
if (Q.lexsorted_degree() || R.lexsorted_degree()){
*logptr(contextptr) << gettext("Cst part must be of type a^n*P(n)") << endl;
remains=e;
return 0;
}
gen a=r2e(Q,v,contextptr)/r2e(R,v,contextptr);
gen P0=r2e(P,v,contextptr);
if (is_zero(P0))
return 0;
for (int i=0;;i++){
gen tmp=normal(subst(e/P0,n,i,false,contextptr),contextptr);
if (!is_undef(tmp)){
P0=tmp*P0/pow(a,i);
break;
}
}
v.clear();
v.push_back(n);
lvar(P0,v);
lvar(l0,v);
lvar(a,v);
vecteur v1(v.begin()+1,v.end());
gen l=e2r(l0,v1,contextptr);
P0=e2r(P0,v,contextptr);
gen P0n,P0d;
fxnd(P0,P0n,P0d);
if (P0n.type!=_POLY)
P=polynome(P0n,int(v.size()));
else
P=*P0n._POLYptr;
a=e2r(a,v1,contextptr);
remains=0;
// solve u(n+1)=l*u(n)+a^n*P(n)
if (a==l){ // let v(n)=u(n)/l^n, then v(n+1)=v(n)+P(n)/l
return pow(l0,n,contextptr)*sum(r2e(P,v,contextptr),n,remains,contextptr)/r2e(l,v1,contextptr)/r2e(P0d,v,contextptr);
}
// search u(n)=a^n*Q(n) with Q a polynomial of same degree than P
// we have u(n+1)=a^(n+1)*Q(n+1)=l*a^n*Q(n)+a^n*P(n)
// hence a*Q(n+1)-l*Q(n)=P(n)
vecteur p=polynome2poly1(P,1);
int pdeg=int(p.size())-1;
vecteur q(pdeg+1);
reverse(p.begin(),p.end());
for (int j=pdeg;j>=0;j--){
gen tmp;
for (int k=j+1;k<=pdeg;++k){
tmp = tmp + comb(k,j)*q[k];
}
q[j]=(p[j]-a*tmp)/(a-l);
}
reverse(q.begin(),q.end());
gen res=r2e(q,v1,contextptr);
if (res.type!=_VECT)
return gentypeerr();
res=symb_horner(*res._VECTptr,n);
res=res*pow(r2e(a,v1,contextptr),n,contextptr)/r2e(P0d,v,contextptr);
return res;
}
/*
void gen2fracpoly1(const gen & g,vecteur & num,vecteur & den){
if (g.type==_FRAC){
num=gen2vecteur(g._FRACptr->num);
den=gen2vecteur(g._FRACptr->den);
}
else {
num=gen2vecteur(g);
den=vecteur(1,1);
}
}
*/
// solveseq(expression,[var,value])
// var is a vector of dim the number of terms in the recurrence
gen _seqsolve(const gen & args,GIAC_CONTEXT){
if ( args.type==_STRNG && args.subtype==-1) return args;
gen f,x=vx_var,uzero=zero,n;
int dim=1,udim=1;
if (args.type==_VECT){
vecteur & v=*args._VECTptr;
int s=int(v.size());
if (!s)
return gentoofewargs("");
f=v[0];
if (s>1)
x=v[1];
if (s>2)
uzero=eval(v[2],eval_level(contextptr),contextptr);
if (x.type==_VECT){
dim=int(x._VECTptr->size());
if (uzero.type==_VECT)
udim=int(uzero._VECTptr->size());
}
}
else
f=args;
x=gen2vecteur(x);
vecteur fv=gen2vecteur(f);
fv=quote_eval(fv,*x._VECTptr,contextptr);
bool fv1=false;
if (fv.size()==1 && fv.front().type==_VECT){
fv1=true;
vecteur tmp=*fv.front()._VECTptr;
fv=tmp;
}
if (dim==udim+1){
n=x._VECTptr->back();
if (n.type!=_IDNT){
identificateur nn(" rsolve_n");
x._VECTptr->back()=nn;
f=quotesubst(f,n,nn,contextptr);
gen res=_seqsolve(gen(makevecteur(f,x,uzero),_SEQ__VECT),contextptr);
return quotesubst(res,nn,n,contextptr);
}
x=vecteur(x._VECTptr->begin(),x._VECTptr->end()-1);
--dim;
}
else {
if (dim!=udim)
return gendimerr();
n=gen(identificateur("rsolve_n"));
}
int fs=int(fv.size());
int add=dim-fs;
if (f.type!=_VECT && !fv1){
f=vecteur(x._VECTptr->end()-add,x._VECTptr->end());
f._VECTptr->push_back(fv[0]);
}
else {
for (int i=0;isize()==1 && x._VECTptr->size()==1 && is_zero(derive(f,n,contextptr))){
// homographic?
gen tmp=ratnormal(f._VECTptr->front(),contextptr),var(x._VECTptr->front());
gen tmpn=_getNum(tmp,contextptr),tmpd=_getDenom(tmp,contextptr),a,b,c,d;
if (is_linear_wrt(tmpn,var,a,b,contextptr) && is_linear_wrt(tmpd,var,c,d,contextptr)&& !is_zero(c)){
// u_{n+1}=(a*u_n+b)/(c*u_n+d)
// find fixed points
gen fixed=solve(a*n+b-(c*n+d)*n,n,1,contextptr);
if (fixed.type==_VECT && fixed._VECTptr->size()==2){
gen r1=fixed._VECTptr->front(),r2=fixed._VECTptr->back();
// (u_n-r1)/(u_n-r2) is geometric, compute ratio
gen un1=(a*n+b)/(c*n+d);
gen r=normal((un1-r1)*(n-r2)/((un1-r2)*(n-r1)),contextptr);
if (is_zero(derive(r,n,contextptr))){
un1=pow(r,n,contextptr)*(var-r1)/(var-r2);
return (r1-un1*r2)/(1-un1);
}
}
if (fixed.type==_VECT && fixed._VECTptr->size()==1){
gen r1=fixed._VECTptr->front();
// 1/(u_n-r1) is arithmetic
gen un1=(a*n+b)/(c*n+d);
gen r=normal(inv(un1-r1,contextptr)-inv(n-r1,contextptr),contextptr);
if (is_zero(derive(r,n,contextptr))){
un1=r*n+inv(var-r1,contextptr);
return inv(un1,contextptr)+r1;
}
}
}
}
return symbolic(at_seqsolve,args);
}
vecteur cst=*normal(subvecteur(*f._VECTptr,multmatvecteur(m,*x._VECTptr)),contextptr)._VECTptr;
if (m.size()==1 && cst.size()==1){
gen l(m[0][0]),c(cst[0]),remains;
gen u0=gen2vecteur(uzero)[0];
if (n.type!=_IDNT)
return gentypeerr();
if (!is_zero(derive(l,n,contextptr))){
vecteur v;
polynome P,Q,R;
if (!is_hypergeometric(l,*n._IDNTptr,v,P,Q,R,contextptr)){
*logptr(contextptr) << gettext("Cst part must be hypergeometric") << endl;
return symbolic(at_seqsolve,args);
}
// l(n+1)/l(n) = P(n+1)/P(n)*Q(n)/R(n+1)
std::vector< monomial >::const_iterator it=Q.coord.begin();
polynome q0=Tnextcoeff(it,Q.coord.end()).untrunc1();
it=R.coord.begin();
polynome r0=Tnextcoeff(it,R.coord.end()).untrunc1();
if (!is_zero(r0*Q-q0*R)){
*logptr(contextptr) << gettext("Unable to handle coeff of homogeneous part") << endl;
return symbolic(at_seqsolve,args);
}
// -> l(n)=l(0)*P(n)/P(0) -> product(l(n))
gen pn=r2e(P,v,contextptr)/r2e(Q,v,contextptr);
gen q0r0=r2e(q0,v,contextptr)/r2e(r0,v,contextptr);
gen res=pow(subst(normal(l/pn,contextptr),n,0,false,contextptr),n,contextptr)*simplify(product(P,v,n,0,n-1,contextptr)/product(Q,v,n,0,n-1,contextptr),contextptr)*pow(q0r0,n*(n-1)/2,contextptr);
// then we might search for a polynomial particular solution to e
if (is_zero(c))
return u0/subst(res,n,0,false,contextptr)*res;
if (!is_one(q0r0))
return gensizeerr("Unable to find particular solution, general solution is "+res.print(contextptr));
// u_{n+1}=l*u_{n}+c
gen tmp=l*x[0]+c,tmpnum,tmpden;
vecteur tmpv(1,n);
lvar(tmp,tmpv);
tmp=e2r(tmp,tmpv,contextptr);
fxnd(tmp,tmpnum,tmpden);
tmpnum=r2e(tmpnum,tmpv,contextptr);
tmpden=r2e(tmpden,tmpv,contextptr);
tmpnum=_e2r(gen(makevecteur(tmpnum,n),_SEQ__VECT),contextptr);
tmpden=_e2r(gen(makevecteur(tmpden,n),_SEQ__VECT),contextptr);
if (is_zero(derive(tmpnum,n,contextptr)) &&
is_zero(derive(tmpnum,n,contextptr)) &&
tmpnum.type==_VECT && tmpden.type==_VECT){
vecteur tmpn,tmpd,ln,cn;
tmpn=*tmpnum._VECTptr;
tmpd=*tmpden._VECTptr;
gen difftmpn=derive(tmpn,x[0],contextptr);
if (is_undef(difftmpn) || difftmpn.type!=_VECT)
return gensizeerr(contextptr);
ln=trim(*difftmpn._VECTptr,0);
cn=trim(subst(tmpn,x[0],0,false,contextptr),0);
if (rsolve_particular(-tmpd,ln,cn,tmpn,contextptr)){
gen sol=_r2e(gen(makevecteur(tmpn,n),_SEQ__VECT),contextptr);
// general solution is sol+C*res, sol(0)+C*res(0)=u0
gen C=subst((u0-sol)/res,n,0,false,contextptr);
return sol+C*res;
}
}
*logptr(contextptr) << gettext("Unable to find a particular solution for inhomogeneous part") << endl;
return symbolic(at_seqsolve,args);
}
gen d=linear_apply(c,n,l,remains,contextptr,rsolve);
if (!is_zero(remains))
*logptr(contextptr) << gettext("Unable to solve recurrence") << endl;
// d is a particular solution of u(n+1)=l*u(n)+c(n)
// add a general solution d(n)+C*l^n
// such that at n=0 we get u0 -> C+d(0)=u0
gen C=normal(u0-quotesubst(d,n,0,contextptr),contextptr);
return d+C*pow(l,n,contextptr);
}
// u(n+1)=M*u(n)+cst
// Let M=PDP^-1, v=P^-1*u, v satisfies v(n+1)=Dv(n)+P^-1*cst
// solve for v then for u
if (!is_zero(derive(m,n,contextptr)))
return symbolic(at_seqsolve,args);
if (has_num_coeff(m))
m=*evalf(m,1,contextptr)._VECTptr;
matrice P,Pinv,D;
bool b=complex_mode(contextptr);
complex_mode(true,contextptr);
egv(m,P,D,contextptr,true,false,false);
complex_mode(b,contextptr);
Pinv=minv(P,contextptr);
if (is_undef(Pinv))
return Pinv;
// if (fs==1) uzero=_revlist(uzero);
vecteur vzero=multmatvecteur(Pinv,*uzero._VECTptr);
vecteur vcst=multmatvecteur(Pinv,cst);
int taille=int(m.size());
vecteur res(taille);
for (int i=taille-1;i>=0;--i){
// find cst coefficient
gen c=vcst[i],l=D[i][i];
for (int j=i+1;jsize()!=u._VECTptr->size())
return gendimerr();
}
else {
if (u.type==_VECT)
return gendimerr();
}
vecteur uv(gen2vecteur(u));
int uvs=int(uv.size());
vecteur initcond;
aplatir(*apply(initcond0,equal2diff)._VECTptr,initcond);
gen f=apply(f0,equal2diff);
if (n.type!=_IDNT){
identificateur N(" rsolve_N");
gen F=quotesubst(f,n,N,contextptr);
gen tmp=crsolve(F,u,N,initcond0,contextptr);
return quotesubst(tmp,N,n,contextptr);
}
vecteur vof0(lop(f,at_of));
// keep only those with u
vecteur vofa,vofb,vopos,vofun;
bool all_one=true,all_bzero=true;
for (const_iterateur it=vof0.begin();it!=vof0.end();++it){
gen & itf=it->_SYMBptr->feuille;
int pos;
if (itf.type==_VECT && itf._VECTptr->size()==2 && (pos=equalposcomp(uv,itf._VECTptr->front())) ){
gen tmp=it->_SYMBptr->feuille._VECTptr->back(),a,b;
if (is_linear_wrt(tmp,n,a,b,contextptr) && !is_zero(a)){
if (!is_one(a))
all_one=false;
if (!is_zero(b))
all_bzero=false;
vofa.push_back(a);
vofb.push_back(b);
vopos.push_back(pos-1);
vofun.push_back(*it);
}
else
return gensizeerr(gettext("Unable to handle this kind of recurrence"));
}
}
if (vofa.empty())
return undef;
if (all_one){
gen bmin=vofb.front(),bmax=vofb.front();
int nmax=0;
for (const_iterateur it=vofb.begin();it!=vofb.end();++it){
const gen & tmp = *it;
if (is_strictly_greater(tmp,bmax,contextptr)){
nmax=int(it-vofb.begin());
bmax=tmp;
}
if (is_strictly_greater(bmin,tmp,contextptr))
bmin=tmp;
}
if (bmin.type!=_INT_ || bmax.type!=_INT_)
return gensizeerr(gettext("Bad indexes in recurrence"));
int m=bmin.val,M=bmax.val;
if (uvs>1){
if (M!=m+1)
return gendimerr(gettext("Multi-recurrences implemented only for 1 step"));
// generate identifiers x0,...,x_{uvs-1} and y0,...,y_{uvs-1}
vecteur idx,idX;
for (int i=0;itype!=_INT_)
return gensizeerr();
bool isx=it->val==m;
int pos=vopos[it-vofb.begin()].val;
vofid.push_back(isx?idx[pos]:idX[pos]);
}
gen F=quotesubst(f,vofun,vofid,contextptr);
if (F.type!=_VECT)
return gensizeerr();
// solves for idX in terms of idx
vecteur Fv=gsolve(*F._VECTptr,idX,1/* complex */,0/* approx mode=no*/,contextptr),res;
if (is_undef(Fv))
return Fv;
for (const_iterateur it=Fv.begin();it!=Fv.end();++it){
vecteur idxn(idx);
idxn.push_back(n);
gen vn=_seqsolve(makevecteur(*it,idxn,idx),contextptr);
gen un=quotesubst(vn,n,n+1-M,contextptr);
if (un.type!=_VECT)
return gensizeerr("Unable to solve this recurrence");
if (int(un._VECTptr->size())!=uvs)
return gendimerr(contextptr);
if (initcond.empty())
return vecteur(1,un);
// solve initial conditions
vecteur vout;
for (int i=0;itype!=_INT_)
return gensizeerr();
vofid.push_back(ids[it->val-m]);
}
gen F=quotesubst(f,vofun,vofid,contextptr),a,b;
if (!is_linear_wrt(F,ids.back(),a,b,contextptr))
return gensizeerr();
F=-b/a; // xM in terms of xm to xM-1
ids.back()=n;
vecteur uinit(ids);
uinit.pop_back();
// let v_n=u_{n+m} -> u_n=v_{n-m}
// old code
// gen vn=_seqsolve(makevecteur(F,ids,uinit),contextptr);
// gen un=quotesubst(vn,n,n-m,contextptr);
// new code that solves rsolve(u(n)=(n)*u(n-1),u(n),u(1)=3);
gen Fm=quotesubst(F,n,n-m,contextptr);
gen un=_seqsolve(makevecteur(Fm,ids,uinit),contextptr);
// end new code
return rsolve_initcond(initcond,un,u,n,uinit,contextptr);
}
// divide and conquier?
if (uvs==1 && vofa.size()==2 && all_bzero){
gen idy=identificateur("rsolve_y"),x(identificateur("rsolve_x"));
vecteur ids(makevecteur(x,idy));
// replace vofun
vecteur vofid;
gen b;
if (is_one(vofa[0])){
vofid=ids;
b=vofa[1];
}
if (is_one(vofa[1])){
vofid=makevecteur(ids[1],ids[0]);
b=vofa[0];
}
if (vofid.empty() || ck_is_greater(1,b,contextptr))
return gensizeerr();
vofun.push_back(n);
vofid.push_back(pow(b,n,contextptr));
gen F=quotesubst(f,vofun,vofid,contextptr);
// F(x,y,n)=0 where y=t(b*n) x=t(n)
// auxiliary sequence u(n)=t(b^n), verifies F(u(n+1),u(n),b^n)=0
gen A,B;
if (is_linear_wrt(F,idy,A,B,contextptr)){
F=-B/A; // u(n+1) in terms of u(n)=x
gen vn=_seqsolve(makevecteur(F,makevecteur(x,n),x),contextptr);
gen un=quotesubst(vn,n,ln(n,contextptr)/ln(b,contextptr),contextptr);
// solve initial cond
return rsolve_initcond(initcond,un,u,n,makevecteur(x),contextptr);
}
}
return gensizeerr(gettext("Not yet implemented"));
}
static gen rsolve(const gen & f0,const gen &u,const gen & n,vecteur & initcond0,int st,GIAC_CONTEXT){
bool b=complex_mode(contextptr);
complex_mode(true,contextptr);
gen res=crsolve(f0,u,n,initcond0,contextptr);
if (!b){
complex_mode(b,contextptr);
// FIXME take real part for res
}
return res;
}
// example args=u(n+1)=2*u(n)+n,u(n),u(0)=1
gen _rsolve(const gen & args,GIAC_CONTEXT){
if ( args.type==_STRNG && args.subtype==-1) return args;
vecteur varg=gen2vecteur(args);
if (debug_infolevel>20)
varg.dbgprint();
int s=int(varg.size());
if (!s)
return gendimerr();
gen f,u,n;
vecteur initcond;
if (s>1){
gen un=varg[1];
if (un.is_symb_of_sommet(at_of)){
gen & unf=un._SYMBptr->feuille;
if (unf.type==_VECT && unf._VECTptr->size()==2){
u=unf._VECTptr->front();
n=unf._VECTptr->back();
}
}
if (un.type==_VECT){
vecteur & unv=*un._VECTptr;
vecteur uv;
for (const_iterateur it=unv.begin();it!=unv.end();++it){
const gen & itg=*it;
if (itg.is_symb_of_sommet(at_of)){
gen & unf=itg._SYMBptr->feuille;
if (unf.type==_VECT && unf._VECTptr->size()==2){
uv.push_back(unf._VECTptr->front());
if (is_zero(n))
n=unf._VECTptr->back();
else
if (n!=unf._VECTptr->back())
return gentypeerr();
}
}
}
u=uv;
}
if (is_zero(n))
return gentypeerr();
}
else {
u=gen(identificateur("rsolve_u"));
n=gen(identificateur("rsolve_n"));
}
vecteur quoted=gen2vecteur(u);
quoted.push_back(n);
varg=quote_eval(varg,quoted,contextptr);
f=varg[0];
if (s>2){
initcond=vecteur(varg.begin()+2,varg.end());
if (initcond.size()==1 && initcond.front().type==_VECT){
// use a temporary vector,
// otherwise the source is destroyed when copying the first element of the vector
// another way to do that could be to have a gen with a copy of initcond.front
vecteur tmp=*initcond.front()._VECTptr;
initcond=tmp;
}
}
else {
if (f.type==_VECT && !f._VECTptr->empty()){
if (u.type!=_VECT || f._VECTptr->size()!=u._VECTptr->size()){
initcond=*f._VECTptr;
f=initcond.front();
initcond.erase(initcond.begin());
}
}
}
int st=step_infolevel(contextptr);
step_infolevel(0,contextptr);
gen res=rsolve(f,u,n,initcond,st,contextptr);
step_infolevel(st,contextptr);
return res;
}
static const char _rsolve_s []="rsolve";
static define_unary_function_eval_quoted (__rsolve,&_rsolve,_rsolve_s);
define_unary_function_ptr5( at_rsolve ,alias_at_rsolve,&__rsolve,_QUOTE_ARGUMENTS,true);
static bool islessthanf(const gen & a,const gen & b){
if (a.type!=_VECT || b.type!=_VECT)
return is_strictly_greater(b,a,context0);
vecteur & av =*a._VECTptr;
vecteur & bv =*b._VECTptr;
int avs=int(av.size()),bvs=int(bv.size());
if (avs!=bvs)
return avs indexbegin,indexsize;
int nindexes=1;
gen initv(vecteur(0));
vecteur & init = *initv._VECTptr;
int shift = xcas_mode(contextptr)!=0 || abs_calc_mode(contextptr)==38;
for (int i=0;ifeuille;
if (f.type!=_VECT || f._VECTptr->size()!=2)
return gendimerr(args[i].print(contextptr));
vecteur & fv= *f._VECTptr;
if (fv[0].type!=_INT_ || fv[1].type!=_INT_ || fv[1].val>24)
return gendimerr(gettext("Array too large")+print_INT_(nindexes));
int is=int(indexsize.size());
for (int i=0;isize()>posj)
initval=(*initval._VECTptr)[posj];
else
initialize=false;
curidx[j] = posj + indexbegin[j];
pos /= indexsize[j];
}
// initialize m[curidx]
if (initialize)
m[curidx]=initval;
else
m[curidx]=0;
}
gen res=m;
res.subtype=1;
return res;
}
static const char _array_s []="array";
static define_unary_function_eval (__array,&_array,_array_s);
define_unary_function_ptr5( at_array ,alias_at_array,&__array,0,true);
static const char _whattype_s []="whattype";
static define_unary_function_eval (__whattype,&_type,_whattype_s);
define_unary_function_ptr5( at_whattype ,alias_at_whattype,&__whattype,0,true);
static const char _modp_s []="modp";
static define_unary_function_eval (__modp,&_irem,_modp_s);
define_unary_function_ptr5( at_modp ,alias_at_modp,&__modp,0,true);
gen _makemod(const gen & args,GIAC_CONTEXT){
if ( args.type==_STRNG && args.subtype==-1) return args;
if (args.type!=_VECT || args._VECTptr->size()!=2)
return gentypeerr();
gen a=args._VECTptr->front(),b=args._VECTptr->back();
if (is_zero(b))
return unmod(a);
if (!is_integer(a) || !is_integer(b))
return gentypeerr();
return makemod(a,b);
}
static const char _makemod_s []="makemod";
static define_unary_function_eval (__makemod,&_makemod,_makemod_s);
define_unary_function_ptr5( at_makemod ,alias_at_makemod,&__makemod,0,true);
gen cpp_convert_0(const gen &g,GIAC_CONTEXT){
return g;
}
longlong cpp_convert_2(const gen & g,GIAC_CONTEXT){
gen h=g;
if (!is_integral(h)){
gensizeerr(contextptr);
return 0;
}
if (h.type==_INT_)
return h.val;
#ifdef USE_GMP_REPLACEMENTS
gensizeerr(contextptr);
return 0;
#else
if (mpz_sizeinbase(*h._ZINTptr,2)>62){
gensizeerr(contextptr);
return 0;
}
if (is_greater(0,g,context0))
return -cpp_convert_2(-g,contextptr);
unsigned int lo, hi;
mpz_t tmp;
mpz_init( tmp );
mpz_mod_2exp( tmp, *h._ZINTptr, 64 ); /* tmp = (lower 64 bits of n) */
lo = mpz_get_ui( tmp ); /* lo = tmp & 0xffffffff */
mpz_div_2exp( tmp, tmp, 32 ); /* tmp >>= 32 */
hi = mpz_get_ui( tmp ); /* hi = tmp & 0xffffffff */
mpz_clear( tmp );
longlong i=(((unsigned long long)hi) << 32) + lo;
return i;
#endif
}
double cpp_convert_1(const gen & g,GIAC_CONTEXT){
gen h=evalf_double(g,1,context0);
if (h.type!=_DOUBLE_){
gensizeerr(contextptr);
return 0;
}
return h._DOUBLE_val;
}
complex cpp_convert_4(const gen & g,GIAC_CONTEXT){
gen h=evalf_double(g,1,context0);
if (h.type==_DOUBLE_)
return h._DOUBLE_val;
if (h.type!=_CPLX || h.subtype!=3){
gensizeerr(contextptr);
return 0;
}
return complex(h._CPLXptr->_DOUBLE_val,(h._CPLXptr+1)->_DOUBLE_val);
}
vecteur cpp_convert_7(const gen & g,GIAC_CONTEXT){
if (g.type!=_VECT){
gensizeerr(contextptr);
return 0;
}
return *g._VECTptr;
}
std::string cpp_convert_12(const gen & g,GIAC_CONTEXT){
if (g.type!=_STRNG){
gensizeerr(contextptr);
return "";
}
return *g._STRNGptr;
}
int is_int_or_double(const gen & g){
if (g.type==_INT_)
return 2;
if (g.type==_DOUBLE_)
return _DOUBLE_;
if (g.type==_IDNT){
const char * ch =g._IDNTptr->id_name;
int s=int(strlen(ch));
if (s>=3 && ch[s-2]=='_') {
switch (ch[s-1]){
case 'i': case 'l':
return 2;
case 'd':
return _DOUBLE_;
}
}
return 0;
}
if (g.type==_SYMB){
if (g._SYMBptr->sommet==at_size)
return 2;
vecteur v=lvar(g._SYMBptr->feuille);
return is_int_or_double(v); // this assumes that the function has the same type as the arguments
}
if (g.type!=_VECT)
return 0;
const vecteur & v=*g._VECTptr;
for (int i=0;isommet==at_size)
return 2;
if (g._SYMBptr->sommet==at_abs && cpp_vartype(g._SYMBptr->feuille)==3)
return 2;
vecteur v=lvar(g._SYMBptr->feuille);
const_iterateur it=v.begin(),itend=v.end();
int cur=-1;
for (;it!=itend;++it){
int i=cpp_vartype(*it);
if (i==0)
return 0;
if (cur==2 && i==_DOUBLE_)
cur=_DOUBLE_;
if (cur==i || (cur==_DOUBLE_ && i==2))
continue;
if (cur==-1)
cur=i;
else
return 0;
}
return cur;
}
if (g.type!=_IDNT)
return 0;
const char * ch=g._IDNTptr->id_name;
int cl=int(strlen(ch));
if (cl>=3 && ch[cl-2]=='_'){
switch (ch[cl-1]){
case 'i':case 'l':
return 2;
case 'd':
return _DOUBLE_;
case 'c':
return _CPLX;
case 'v':
return _VECT;
case 's':
return _STRNG;
}
}
return 0;
}
std::string cprintvars(const gen & f0,bool cst,const string & sep,int * typeptr,GIAC_CONTEXT){
vecteur fv=gen2vecteur(f0);
string res;
for (int i=0;;){
gen fvi=fv[i];
bool eq=fvi.is_symb_of_sommet(at_equal);
bool eqsto=fvi.is_symb_of_sommet(at_sto);
if (eq)
fvi=fvi[1];
if (eqsto)
fvi=fvi[2];
if (fvi.type!=_IDNT)
return "Invalid parameter "+fv[i].print(contextptr);
const char * ch=fvi._IDNTptr->id_name;
int cl=int(strlen(ch));
string vtype=cst?"const giac::gen & ":"giac::gen ";
if (typeptr) *typeptr=0;
if (cl>=3 && ch[cl-2]=='_'){
switch (ch[cl-1]){
case 'i': case 'l':
vtype="long ";
if (typeptr) *typeptr=2;
break;
case 'd':
vtype="double ";
if (typeptr) *typeptr=_DOUBLE_;
break;
case 'c':
vtype="complex ";
if (typeptr) *typeptr=_CPLX;
break;
case 'v':
vtype="giac::vecteur ";
if (typeptr) *typeptr=_VECT;
break;
case 's':
vtype="std::string ";
if (typeptr) *typeptr=_STRNG;
break;
}
}
if (typeptr) ++typeptr;
vtype = vtype+ch;
if (eq || eqsto){
gen g=(fv[i])[1];
if (eq)
g=(fv[i])[2];
#if 1
int i=cpp_vartype(fvi);
if (!i || i==cpp_vartype(g))
vtype = vtype + " = " +cprint(g,0,contextptr);
else
vtype = vtype + " = cpp_convert_"+print_INT_(i)+"("+cprint(g,0,contextptr)+",contextptr)";
#else
vtype = vtype+'=';
int i=is_int_or_double(fvi);
if (!i || is_int_or_double(g))
vtype = vtype +g.print(contextptr);
else {
vtype = vtype + "("+g.print(contextptr);
vtype += (i==_DOUBLE_)?")._DOUBLE_val":").val";
}
#endif
}
++i;
if (i==int(fv.size()))
return res+vtype;
res = res + vtype + sep;
}
}
std::string cpp_stoprint(const gen & g,GIAC_CONTEXT){
if (g.is_symb_of_sommet(at_at) && g._SYMBptr->feuille.type==_VECT && g._SYMBptr->feuille._VECTptr->size()==2){
gen f=g._SYMBptr->feuille._VECTptr->front();
gen i=g._SYMBptr->feuille._VECTptr->back();
int t=cpp_vartype(i);
if (t!=2)
return cprint(f,0,contextptr)+"[cpp_convert_2("+cprint(i,0,contextptr)+",contextptr)]";
}
return cprint(g,0,contextptr);
}
// name is the program name if args==program(..) or 1 for vectors of instructions or 0 otherwise or -1 for vector to put in a makesequence
std::string cprint(const gen & args,const gen & name,GIAC_CONTEXT){
if (args.type==_VECT){
string sep=name==1?";\n":",";
string s;
if (name==0 && args.subtype!=_SEQ__VECT)
s="giac::makevecteur(";
if (name==-1)
s="giac::makesequence(";
vecteur::const_iterator it=args._VECTptr->begin(), itend=args._VECTptr->end();
if (it==itend)
return "gen(vecteur(0),"+print_INT_(args.subtype)+");";
for(int i=0;;++i){
s += cprint(*it,(name.type==_VECT && isize())?name[i]:zero,contextptr);
++it;
if (it==itend){
break;
}
s += sep;
}
if (name==0 && args.subtype!=_SEQ__VECT) s += ")";
if (name==-1) s += ")";
if (name==1) s += ";";
return s;
}
if (args.type==_CPLX)
return "gen("+args._CPLXptr->print(contextptr)+","+(args._CPLXptr+1)->print(contextptr)+")";
if (args.type==_FUNC){
if (args==at_break)
return "break";
if (args==at_continue)
return "continue";
string name=args.print(contextptr);
if (name.size()>2 && name[0]=='\'' && name[name.size()-1]=='\'')
name=name.substr(1,name.size()-2);
return "at_"+name;
}
if (args.type!=_SYMB)
return args.print(contextptr);
unary_function_ptr & u=args._SYMBptr->sommet;
if (u==at_break)
return "break";
if (u==at_continue)
return "continue";
gen f=args._SYMBptr->feuille;
if ( (u==at_sto || u==at_array_sto) && f.type==_VECT && f._VECTptr->size()==2){
#if 1
int i=cpp_vartype(f[1]);
if (!i || i==cpp_vartype(f[0]))
return cpp_stoprint(f[1],contextptr)+'='+cprint(f[0],0,contextptr);;
string res= cpp_stoprint(f[1],contextptr)+" = cpp_convert_"+print_INT_(i)+"("+cprint(f[0],0,contextptr)+",contextptr)";
return res;
#else
int i=is_int_or_double(f[1]);
if (!i || is_int_or_double(f[0]))
return f[1].print(contextptr)+'='+cprint(f[0],0,contextptr);;
string res= f[1].print(contextptr)+"=("+cprint(f[0],0,contextptr);
res += (i==_DOUBLE_)?")._DOUBLE_val":").val";
return res;
#endif
}
if (u==at_equal && f.type==_VECT && f._VECTptr->size()==2){
#if 1
int i=cpp_vartype(f[0]);
if (!i || i==cpp_vartype(f[1]))
return cpp_stoprint(f[0],contextptr)+'='+cprint(f[1],0,contextptr);
string res= cpp_stoprint(f[0],contextptr)+" = cpp_convert_"+print_INT_(i)+"("+cprint(f[1],0,contextptr)+",contextptr)";
return res;
#else
int i=is_int_or_double(f[0]);
if (!i || is_int_or_double(f[1]))
return f[0].print(contextptr)+'='+cprint(f[1],0,contextptr);;
string res= f[0].print(contextptr)+"=("+cprint(f[1],0,contextptr);
res += (i==_DOUBLE_)?")._DOUBLE_val":").val";
return res;
#endif
}
if (u==at_size){
int i=cpp_vartype(f);
if (i==_VECT)
return cprint(f,0,contextptr)+".size()";
return "cpp_convert_7("+cprint(f,0,contextptr)+",contextptr).size()";
}
if (u==at_at && f.type==_VECT && f._VECTptr->size()==2){
int i=cpp_vartype(f._VECTptr->front());
int j=cpp_vartype(f._VECTptr->back());
if (i==0 || j==2)
return f._VECTptr->front().print(contextptr)+"["+cprint(f._VECTptr->back(),0,contextptr)+"]";
return f._VECTptr->front().print(contextptr)+"[cpp_convert_2("+cprint(f._VECTptr->back(),0,contextptr)+",contextptr)]";
}
if (u==at_program && f.type==_VECT && f._VECTptr->size()==3){
// header
bool head=name.type==_IDNT;
string heads;
int convert=0;
if (head){
// IMPROVE if name=="main" make it argv,argc->int
string rtype="giac::gen ";
convert=cpp_vartype(name);
if (convert){
switch (convert){
case _DOUBLE_:
rtype="double ";
break;
case 2:
rtype="long ";
break;
case _VECT:
rtype="giac::vecteur ";
break;
case _STRNG:
rtype="std::string ";
break;
}
}
heads=rtype+name._IDNTptr->id_name+"(";
}
else
heads="giac::gen f(";
vector argtype(gen2vecteur(f[0]).size());
heads += cprintvars(f[0],false,",",&argtype.front(),contextptr)+",const giac::context * contextptr=0)";
// core
string core=cprint(f[2],0,contextptr);
string res= heads+"{\n"+core+"\n}\n";
heads="f";
if (head)
heads=name._IDNTptr->id_name;
if (f[0].type!=_VECT || f[0]._VECTptr->size()<=1){
if (f[0].type==_VECT && f[0]._VECTptr->empty())
core ="return "+heads+"(contextptr);";
else {
#if 1
core = "return "+heads+"(cpp_convert_"+print_INT_(cpp_vartype(f[0]._VECTptr->front()))+"(g,contextptr),contextptr);";
#else
switch (convert){
case 0:
core="return "+heads+"(g,contextptr);";
break;
case 2:
core="if (g.type!=_INT_) return gensizeerr(contextptr); else return "+heads+"(g.val,contextptr);";
break;
case _DOUBLE_:
core="if (g.type!=_DOUBLE_) return gensizeerr(contextptr); else return "+heads+"(g._DOUBLE_val,contextptr);";
break;
}
#endif
}
}
else {
int nargs=int(f[0]._VECTptr->size());
core = "if (g.type!=_VECT || g.subtype!=_SEQ__VECT || g._VECTptr->size()!="+print_INT_(nargs)+") return gendimerr(contextptr);\n";
core += "vecteur v = *g._VECTptr;";
for (int i=0;isize()==4){
string res="for(";
res = res +cprint(f[0],0,contextptr);
res = res +';';
res = res +cprint(f[1],0,contextptr);
res = res +';';
res = res +cprint(f[2],0,contextptr);
res = res +"){\n";
res = res +cprint(f[3],0,contextptr);
res = res +";\n}\n";
return res;
}
if (u==at_bloc)
return cprint(f,1,contextptr);
if ( (u==at_ifte || u==at_si) && f.type==_VECT && f._VECTptr->size()==3){
string res="if(";
res = res +cprint(f[0],0,contextptr);
res = res +"){\n";
res = res +cprint(f[1],0,contextptr);
res = res +";\n}";
if (!is_zero(f[2])){
res = res + " else {\n";
res = res + cprint(f[2],0,contextptr) + ";\n}";
}
return res;
}
if (u==at_try_catch && f.type==_VECT && f._VECTptr->size()>=3){
string res="try {\n";
res = res +cprint(f[0],0,contextptr);
res = res +";\n} catch(";
res = res +cprint(f[1],0,contextptr);
res = res +"){\n";
res = res +cprint(f[2],0,contextptr);
res = res +";\n}";
return res;
}
if (u==at_return){
if (f.type==_VECT)
return "return giac::makevecteur("+cprint(f,0,contextptr)+")";
else
return "return "+cprint(f,0,contextptr);
}
if (u.ptr()->cprint)
return u.ptr()->cprint(args._SYMBptr->feuille,u.ptr()->s,contextptr);
// operators: if all vars are double or int args.print()
// otherwise or not an operator giac::opname(,contextptr)
string res=u.ptr()->print(contextptr);
int idf0=is_int_or_double(f);
bool idf= idf0 || cpp_vartype(f)==_CPLX;
if (u==at_sin || u==at_cos || u==at_tan || u==at_asin || u==at_acos || u==at_atan || u==at_exp || u==at_ln || u==at_abs){
if (idf)
return "std::"+args.print(contextptr);
}
else
res = '_'+res;
if (u.ptr()->printsommet==printsommetasoperator || u==at_division ||
(idf0 && (u==at_inferieur_strict || u==at_inferieur_egal || u==at_superieur_strict || u==at_superieur_egal))){
int rs=int(res.size());
if ( (idf && (rs!=2 || res[1]!='^')) || (rs==2 && (res[1]=='+' || res[1]=='*' || res[1]=='>' || res[1]=='<')) || (rs==3 && (res[1]=='>' || res[1]=='<') && res[2]=='=') ){
res=u.ptr()->print(contextptr);
string sres;
vecteur v=gen2vecteur(f);
for (int i=0;;){
bool b=need_parenthesis(v[i]);
if (b) sres += '(';
sres += cprint(v[i],0,contextptr);
if (b) sres += ')';
++i;
if (i>=v.size())
break;
sres += res;
}
return sres;
}
if (rs==2){
switch (res[1]){
case '+':
res="_plus";
break;
case '*':
res="_prod";
break;
case '/':
res="_division";
break;
case '<':
res="_inferieur_strict";
break;
case '>':
res="_superieur_strict";
break;
case '=':
res="_equal";
break;
case '^':
res="_pow";
break;
}
}
}
if (u==at_increment || u==at_decrement){
if (f.type!=_VECT && is_int_or_double(f)){
return cprint(f,0,contextptr)+(u==at_increment?"++":"--");
}
if (f.type==_VECT && f._VECTptr->size()==2 && is_int_or_double(f._VECTptr->front()) && is_int_or_double(f._VECTptr->back()) ){
return cprint(f._VECTptr->front(),0,contextptr)+(u==at_increment?"+=":"-=")+cprint(f._VECTptr->back(),0,contextptr);
}
}
if (u==at_same || u==at_different || u==at_and || u==at_ou ){
return cprint(args._SYMBptr->feuille[0],-1,contextptr)+u.ptr()->s+cprint(args._SYMBptr->feuille[1],-1,contextptr);
}
if (u==at_neg || u==at_not)
return u.ptr()->s+cprint(args._SYMBptr->feuille,-1,contextptr);
if (u==at_inferieur_egal)
return "giac::is_greater("+cprint(args._SYMBptr->feuille[1],-1,contextptr)+","+cprint(args._SYMBptr->feuille[0],0,contextptr)+",contextptr)";
if (u==at_superieur_egal)
return "giac::is_greater("+cprint(args._SYMBptr->feuille[0],-1,contextptr)+","+cprint(args._SYMBptr->feuille[1],0,contextptr)+",contextptr)";
if (u==at_inferieur_strict)
return "!giac::is_greater("+cprint(args._SYMBptr->feuille[0],-1,contextptr)+","+cprint(args._SYMBptr->feuille[1],0,contextptr)+",contextptr)";
if (u==at_superieur_strict)
return "!giac::is_greater("+cprint(args._SYMBptr->feuille[1],-1,contextptr)+","+cprint(args._SYMBptr->feuille[0],0,contextptr)+",contextptr)";
res = "giac::"+res + '(';
res = res + cprint(args._SYMBptr->feuille,-1,contextptr);
res = res +",contextptr)";
return res;
// ...
}
std::string cprint(const gen & args,GIAC_CONTEXT){
return cprint(args,0,contextptr);
}
gen _cprint(const gen & args,GIAC_CONTEXT){
if ( args.type==_STRNG && args.subtype==-1) return args;
int x=xcas_mode(contextptr),l=language(contextptr);
xcas_mode(0,contextptr);
language(2,contextptr);
gen args_evaled=eval(args,1,contextptr);
string s=cprint(args_evaled,args,contextptr);
xcas_mode(x,contextptr);
language(l,contextptr);
return string2gen(s,false);
}
static const char _cprint_s []="cprint";
static define_unary_function_eval (__cprint,&_cprint,_cprint_s);
define_unary_function_ptr5( at_cprint ,alias_at_cprint,&__cprint,_QUOTE_ARGUMENTS,true);
#ifndef GIAC_HAS_STO_38
// make a cpp translation file
gen _cpp(const gen & args,GIAC_CONTEXT){
if ( args.type==_STRNG && args.subtype==-1) return args;
gen tmp=check_secure();
if (is_undef(tmp)) return tmp;
if (args.type==_IDNT){
int x=xcas_mode(contextptr),l=language(contextptr);
xcas_mode(0,contextptr);
language(2,contextptr);
gen args_evaled=eval(args,1,contextptr);
string s=cprint(args_evaled,args,contextptr);
language(l,contextptr);
xcas_mode(x,contextptr);
string funcname=args.print(contextptr);
string filename="giac_"+funcname+".cpp";
ofstream of(filename.c_str());
#ifdef __APPLE__
of << "// -*- mode:c++; compile-command:\" c++ -I/Applications/usr/include -I.. -I. -fPIC -DPIC -g -O2 -dynamiclib giac_" << funcname << ".cpp -o libgiac_" << funcname << ".dylib -L/Applications/usr/lib -L/Applications/usr/64/local/lib -lgiac \" -*-" << endl;
#else
of << "// -*- mode:c++; compile-command:\" c++ -I.. -I. -fPIC -DPIC -g -O2 -c giac_" << funcname << ".cpp -o giac_" << funcname << ".lo && cc -shared giac_"<" << endl;
of << "#include " << endl;
of << "using namespace std;" << endl;
of << "namespace giac {" << endl;
of << "inline bool operator <= (const gen & a,const gen &b){ return is_greater(b,a,giac::context0); }\ninline bool operator >= (const gen & a,const gen &b){ return is_greater(a,b,giac::context0); }" <size()!=2 || args._VECTptr->back().type!=_STRNG)
return gensizeerr(contextptr);
int x=xcas_mode(contextptr);
xcas_mode(0,contextptr);
gen arg=args._VECTptr->front();
gen args_evaled=eval(arg,1,contextptr);
string s=cprint(args_evaled,arg,contextptr);
xcas_mode(x,contextptr);
string filename=*args._VECTptr->back()._STRNGptr;
ofstream of(filename.c_str());
of << "#include " << endl;
of << s << endl;
of.close();
string cmd="indent -br -brf "+filename;
system_no_deprecation(cmd.c_str());
return 1;
}
static const char _cpp_s []="cpp";
static define_unary_function_eval (__cpp,&_cpp,_cpp_s);
define_unary_function_ptr5( at_cpp ,alias_at_cpp,&__cpp,_QUOTE_ARGUMENTS,true);
#endif
gen _hexprint(const gen & g,GIAC_CONTEXT){
if ( g.type==_STRNG && g.subtype==-1) return g;
if (g.type == _INT_)
return string2gen(hexa_print_INT_(g.val), false);
if (g.type == _ZINT)
return string2gen(hexa_print_ZINT(*g._ZINTptr), false);
return gentypeerr();
}
static const char _hexprint_s []="hexprint";
static define_unary_function_eval (__hexprint,&_hexprint,_hexprint_s);
define_unary_function_ptr5( at_hexprint ,alias_at_hexprint,&__hexprint,0,true);
gen _octprint(const gen & g,GIAC_CONTEXT){
if ( g.type==_STRNG && g.subtype==-1) return g;
if (g.type == _INT_)
return string2gen(octal_print_INT_(g.val), false);
if (g.type == _ZINT)
return string2gen(octal_print_ZINT(*g._ZINTptr), false);
return gentypeerr();
}
static const char _octprint_s []="octprint";
static define_unary_function_eval (__octprint,&_octprint,_octprint_s);
define_unary_function_ptr5( at_octprint ,alias_at_octprint,&__octprint,0,true);
gen _binprint(const gen & g,GIAC_CONTEXT){
if ( g.type==_STRNG && g.subtype==-1) return g;
if (g.type == _INT_)
return string2gen(binary_print_INT_(g.val), false);
if (g.type == _ZINT)
return string2gen(binary_print_ZINT(*g._ZINTptr), false);
return gentypeerr();
}
static const char _binprint_s[]="binprint";
static define_unary_function_eval(__binprint,&_binprint,_binprint_s);
define_unary_function_ptr5(at_binprint,alias_at_binprint,&__binprint,0,true);
// const unary_function_ptr at_binprint (&__binprint,0,true);
#ifndef NO_NAMESPACE_GIAC
} // namespace giac
#endif // ndef NO_NAMESPACE_GIAC