1 // -*- mode:C++ ; compile-command: "g++ -DHAVE_CONFIG_H -m32 -I. -I.. -DIN_GIAC -DGIAC_GENERIC_CONSTANTS -g -c -fno-strict-aliasing moyal.cc" -*- 2 #include "giacPCH.h" 3 /* 4 * Copyright (C) 2000,2014 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 5 * 6 * This program is free software; you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation; either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. 18 */ 19 using namespace std; 20 #include <stdexcept> 21 #include <cmath> 22 #include <math.h> 23 #include <cstdlib> 24 #include "sym2poly.h" 25 #include "usual.h" 26 #include "moyal.h" 27 #include "solve.h" 28 #include "intg.h" 29 #include "permu.h" 30 #include "giacintl.h" 31 #ifdef HAVE_LIBGSL 32 #include <gsl/gsl_sf_airy.h> 33 #include <gsl/gsl_sf_erf.h> 34 #endif 35 36 #ifndef NO_NAMESPACE_GIAC 37 namespace giac { 38 #endif // ndef NO_NAMESPACE_GIAC 39 40 #ifdef HAVE_SSCL moyal(const gen & a,const gen & b,const gen & vars,const gen & order)41 gen moyal(const gen & a,const gen & b, const gen &vars,const gen & order){ 42 return symb_moyal(a,b,vars,order); 43 } 44 45 #else // HAVE_SSCL 46 gen moyal(const gen & a,const gen & b, const gen &vars,const gen & order){ 47 return symb_moyal(a,b,vars,order); 48 } 49 50 #endif // HAVE_SSCL 51 52 // "unary" version _moyal(const gen & args,GIAC_CONTEXT)53 gen _moyal(const gen & args,GIAC_CONTEXT){ 54 if ( args.type==_STRNG && args.subtype==-1) return args; 55 int s=int(args._VECTptr->size()); 56 if (s!=4) return gensizeerr(gettext("moyal.cc/_moyal")); 57 return moyal( (*(args._VECTptr))[0],(*(args._VECTptr))[1],(*(args._VECTptr))[2],(*(args._VECTptr))[3]); 58 } 59 static const char _moyal_s []="moyal"; texprintasmoyal(const gen & g,const char * s,GIAC_CONTEXT)60 static string texprintasmoyal(const gen & g,const char * s,GIAC_CONTEXT){ 61 return texprintsommetasoperator(g,"#",contextptr); 62 } 63 static define_unary_function_eval4 (__moyal,&_moyal,_moyal_s,0,&texprintasmoyal); 64 define_unary_function_ptr5( at_moyal ,alias_at_moyal,&__moyal,0,true); 65 incomplete_beta(double a,double b,double p,bool regularize)66 gen incomplete_beta(double a,double b,double p,bool regularize){ // regularize=true by default 67 // I_p(a,b)=1/B(a,b)*int(t^(a-1)*(1-t)^(b-1),t=0..p) 68 // =p^a*(1-p)^(b-1)/B(a,b)*continued fraction expansion 69 // 1/(1+e2/(1+e3/(1+...))) 70 // e_2m=-(a+m-1)*(b-m)/(a+2*m-2)/(a+2*m-1)*(x/(1-x)) 71 // e_2m+1= m*(a+b-1+m)/(a+2*m-1)/(a+2*m)*(x/(1-x)) 72 // assumes p in [0,1] 73 if (p<=0) 74 return 0; 75 if (a<=0 || b<=0) 76 return 1; 77 // add here test for returning 1 if b>>a 78 if (p>a/double(a+b)){ 79 gen tmp=incomplete_beta(b,a,1-p,true); 80 if (regularize) 81 return 1-tmp; 82 return Beta(a,b,context0)*(1-tmp); 83 } 84 // Continued fraction expansion: a1/(b1+a2/(b2+...))) 85 // P0=1, P1=a1, Q0=1, Q1=b1 86 // j>=2: Pj=bj*Pj-1+aj*Pj-2, Qj=bj*Qj-1+aj*Qj-2 87 // Here bm=1, am=em, etc. 88 long_double Pm2=0,Pm1=1,Pm,Qm2=1,Qm1=1,Qm,am,x=p/(1-p); 89 long_double deux=9007199254740992.,invdeux=1/deux; 90 for (long_double m=1;m<100;++m){ 91 // odd term 92 am=-(a+m-1)*(b-m)/(a+2*m-2)/(a+2*m-1)*x; 93 Pm=Pm1+am*Pm2; 94 Qm=Qm1+am*Qm2; 95 Pm2=Pm1; Pm1=Pm; 96 Qm2=Qm1; Qm1=Qm; 97 // even term 98 am=m*(a+b-1+m)/(a+2*m-1)/(a+2*m)*x; 99 Pm=Pm1+am*Pm2; 100 Qm=Qm1+am*Qm2; 101 // cerr << Pm/Qm << " " << Pm2/Qm2 << "\n"; 102 if (absdouble(Pm/Qm-Pm2/Qm2)<1e-16*absdouble(Pm/Qm)){ 103 double res=Pm/Qm; 104 #if 0 // def VISUALC // no lgamma available 105 gen r=res/a*std::pow(p,a)*std::pow(1-p,b-1); 106 if (regularize) 107 r=r*Gamma(a+b,context0)/Gamma(a,context0)/Gamma(b,context0); 108 return r; 109 #else 110 if (regularize) 111 return res/a*std::exp(a*std::log(p)+(b-1)*std::log(1-p)+lngamma(a+b)-lngamma(a)-lngamma(b)); 112 return res/a*std::exp(a*std::log(p)+(b-1)*std::log(1-p)); 113 #endif 114 } 115 Pm2=Pm1; Pm1=Pm; 116 Qm2=Qm1; Qm1=Qm; 117 if (absdouble(Pm)>deux){ 118 Pm2 *= invdeux; Qm2 *= invdeux; Pm1 *= invdeux; Qm1 *= invdeux; 119 } 120 if (absdouble(Pm)<invdeux){ 121 Pm2 *= deux; Qm2 *= deux; Pm1 *= deux; Qm1 *= deux; 122 } 123 } 124 return undef; //error 125 } 126 beta_mult(gen & res,gen & a,GIAC_CONTEXT)127 static void beta_mult(gen &res,gen & a,GIAC_CONTEXT){ 128 for (;;){ 129 gen a1=a-1; 130 if (!is_positive(a1,contextptr)) 131 return; 132 res=a1*res; 133 a=a1; 134 } 135 } Beta(const gen & a,const gen & b,GIAC_CONTEXT)136 gen Beta(const gen & a,const gen& b,GIAC_CONTEXT){ 137 if (a.type==_DOUBLE_ || b.type==_DOUBLE_ || 138 a.type==_FLOAT_ || b.type==_FLOAT_ || 139 a.type==_CPLX || b.type==_CPLX ){ 140 gen A=evalf_double(a,1,contextptr); 141 gen B=evalf_double(b,1,contextptr); 142 gen C=lngamma(A+B,contextptr); 143 A=lngamma(A,contextptr); 144 B=lngamma(B,contextptr); 145 C=A+B-C; 146 C=exp(C,contextptr); 147 return C; 148 } 149 gen n; 150 if (a.type==_FRAC && b.type==_FRAC && is_positive(a,contextptr) && is_positive(b,contextptr) && is_integer( (n=a+b) )){ 151 gen res=1,a_(a),b_(b); 152 beta_mult(res,a_,contextptr); 153 beta_mult(res,b_,contextptr); 154 if (a_+b_==1){ 155 return ratnormal(res*cst_pi/sin(cst_pi*a_,contextptr)/Gamma(n,contextptr),contextptr); 156 } 157 } 158 return Gamma(a,contextptr)*Gamma(b,contextptr)/Gamma(a+b,contextptr); 159 } _Beta(const gen & args,GIAC_CONTEXT)160 gen _Beta(const gen & args,GIAC_CONTEXT){ 161 if ( args.type==_STRNG && args.subtype==-1) return args; 162 if (args.type!=_VECT) 163 return symbolic(at_Beta,args); 164 vecteur v=*args._VECTptr; 165 int s=int(v.size()); 166 if (s>2 && (v[0].type==_DOUBLE_ || v[1].type==_DOUBLE_ || v[2].type==_DOUBLE_ || v[0].type==_REAL || v[1].type==_REAL || v[2].type==_REAL)){ 167 gen tmp=evalf_double(v,1,contextptr); 168 if (tmp.type==_VECT) 169 v=*tmp._VECTptr; 170 s=int(v.size()); 171 } 172 if ( (s==3 || s==4) && v[0].type==_DOUBLE_ && v[1].type==_DOUBLE_ && v[2].type==_DOUBLE_ ){ 173 return incomplete_beta(v[0]._DOUBLE_val,v[1]._DOUBLE_val,v[2]._DOUBLE_val, s==4 && !is_zero(v[3]) ); 174 } 175 if (s<2 || s>4) 176 return gendimerr(contextptr); 177 if (s==4){ 178 if (is_zero(v[3])) 179 return symbolic(at_Beta,makesequence(v[0],v[1],v[2])); 180 return symbolic(at_Beta,makesequence(v[0],v[1],v[2]))/Beta(v[0],v[1],contextptr); 181 } 182 if (s!=2) 183 return symbolic(at_Beta,args); 184 return Beta(v[0],v[1],contextptr); 185 } 186 static const char _Beta_s []="Beta"; 187 static define_unary_function_eval (__Beta,&_Beta,_Beta_s); 188 define_unary_function_ptr5( at_Beta ,alias_at_Beta,&__Beta,0,true); 189 _upper_incomplete_gamma(const gen & args,GIAC_CONTEXT)190 gen _upper_incomplete_gamma(const gen & args,GIAC_CONTEXT){ 191 if ( args.type==_STRNG && args.subtype==-1) return args; 192 if (args.type!=_VECT) 193 return symbolic(at_upper_incomplete_gamma,args); 194 vecteur v=*args._VECTptr; 195 int s=int(v.size()); 196 if (s>=2 && (v[0].type==_DOUBLE_ || v[1].type==_DOUBLE_)){ 197 v[0]=evalf_double(v[0],1,contextptr); 198 v[1]=evalf_double(v[1],1,contextptr); 199 } 200 if ( (s==2 || s==3) && v[0].type==_DOUBLE_ && v[1].type==_DOUBLE_ ){ 201 double res=upper_incomplete_gammad(v[0]._DOUBLE_val,v[1]._DOUBLE_val,s==3?!is_zero(v[2]):false); 202 if (res==-1){ 203 if (s==3 && !is_zero(v[2])) 204 return 1-lower_incomplete_gamma(v[0]._DOUBLE_val,v[1]._DOUBLE_val,true,contextptr); 205 return Gamma(v[0]._DOUBLE_val,contextptr)-lower_incomplete_gamma(v[0]._DOUBLE_val,v[1]._DOUBLE_val,false,contextptr); 206 // return gensizeerr(contextptr); 207 } 208 return res; 209 } 210 if (s<2 || s>3) 211 return gendimerr(contextptr); 212 if (abs_calc_mode(contextptr)!=38) // check may be removed if ugamma declared 213 return symbolic(at_upper_incomplete_gamma,args); 214 if (s==3){ 215 if (is_zero(v[2])) 216 return Gamma(v[0],contextptr)-symbolic(at_lower_incomplete_gamma,makesequence(v[0],v[1])); 217 return symbolic(at_Gamma,makesequence(v[0],v[1],1)); 218 } 219 return symbolic(at_upper_incomplete_gamma,args); 220 } 221 static const char _upper_incomplete_gamma_s []="ugamma"; 222 static define_unary_function_eval (__upper_incomplete_gamma,&_upper_incomplete_gamma,_upper_incomplete_gamma_s); 223 define_unary_function_ptr5( at_upper_incomplete_gamma ,alias_at_upper_incomplete_gamma,&__upper_incomplete_gamma,0,true); 224 _polygamma(const gen & args,GIAC_CONTEXT)225 gen _polygamma(const gen & args,GIAC_CONTEXT){ 226 if ( args.type==_STRNG && args.subtype==-1) return args; 227 if (args.type!=_VECT) 228 return symbolic(at_polygamma,args); 229 vecteur v=*args._VECTptr; 230 int s=int(v.size()); 231 if (args.subtype==_SEQ__VECT && s==2) 232 return _Psi(makesequence(v[1],v[0]),contextptr); 233 return symbolic(at_polygamma,args); 234 } 235 static const char _polygamma_s []="polygamma"; 236 static define_unary_function_eval (__polygamma,&_polygamma,_polygamma_s); 237 define_unary_function_ptr5( at_polygamma ,alias_at_polygamma,&__polygamma,0,true); 238 Airy_Ai(const gen & x,GIAC_CONTEXT)239 gen Airy_Ai(const gen & x,GIAC_CONTEXT){ 240 gen e=x.evalf(1,contextptr); 241 #ifdef HAVE_LIBGSL 242 if (e.type==_DOUBLE_) 243 return gsl_sf_airy_Ai(e._DOUBLE_val,GSL_PREC_DOUBLE); 244 #endif 245 return symbolic(at_Airy_Ai,x); 246 } _Airy_Ai(const gen & args,GIAC_CONTEXT)247 gen _Airy_Ai(const gen & args,GIAC_CONTEXT){ 248 if ( args.type==_STRNG && args.subtype==-1) return args; 249 return apply(args,Airy_Ai,contextptr); 250 } 251 static const char _Airy_Ai_s []="Airy_Ai"; 252 static define_unary_function_eval (__Airy_Ai,&_Airy_Ai,_Airy_Ai_s); 253 define_unary_function_ptr5( at_Airy_Ai ,alias_at_Airy_Ai,&__Airy_Ai,0,true); 254 Airy_Bi(const gen & x,GIAC_CONTEXT)255 gen Airy_Bi(const gen & x,GIAC_CONTEXT){ 256 gen e=x.evalf(1,contextptr); 257 #ifdef HAVE_LIBGSL 258 if (e.type==_DOUBLE_) 259 return gsl_sf_airy_Bi(e._DOUBLE_val,GSL_PREC_DOUBLE); 260 #endif 261 return symbolic(at_Airy_Bi,x); 262 } _Airy_Bi(const gen & args,GIAC_CONTEXT)263 gen _Airy_Bi(const gen & args,GIAC_CONTEXT){ 264 if ( args.type==_STRNG && args.subtype==-1) return args; 265 return apply(args,Airy_Bi,contextptr); 266 } 267 static const char _Airy_Bi_s []="Airy_Bi"; 268 static define_unary_function_eval (__Airy_Bi,&_Airy_Bi,_Airy_Bi_s); 269 define_unary_function_ptr5( at_Airy_Bi ,alias_at_Airy_Bi,&__Airy_Bi,0,true); 270 _UTPN(const gen & args,GIAC_CONTEXT)271 gen _UTPN(const gen & args,GIAC_CONTEXT){ 272 if ( args.type==_STRNG && args.subtype==-1) return args; 273 if (args.type!=_VECT) 274 return erfc(args/plus_sqrt2,contextptr)/2; 275 vecteur & v=*args._VECTptr; 276 int s=int(v.size()); 277 if (s!=3 || is_zero(v[1])) 278 return gensizeerr(contextptr); 279 return erfc((v[2]-v[0])/sqrt(2*v[1],contextptr),contextptr)/2; 280 } 281 static const char _UTPN_s []="UTPN"; 282 static define_unary_function_eval (__UTPN,&_UTPN,_UTPN_s); 283 define_unary_function_ptr5( at_UTPN ,alias_at_UTPN,&__UTPN,0,true); 284 285 #ifndef USE_GMP_REPLACEMENTS randdiscrete(const vecteur & m,GIAC_CONTEXT)286 gen randdiscrete(const vecteur &m,GIAC_CONTEXT) { 287 int n; 288 if (m.empty() || !m.front().is_integer() || (n=m.front().val)<1) 289 return gensizeerr(contextptr); 290 double ran1=giac_rand(contextptr)/(rand_max2+1.0); 291 double ran2=giac_rand(contextptr)/(rand_max2+1.0); 292 int i=std::floor(n*ran1); 293 int index=is_strictly_greater(m[i+1],ran2,contextptr)?i:m[n+i+1].val; 294 if (int(m.size()-1)==3*n) 295 return m[2*n+index+1]; 296 return index+array_start(contextptr); 297 } 298 _discreted(const gen & g,GIAC_CONTEXT)299 gen _discreted(const gen &g,GIAC_CONTEXT) { 300 if (g.type==_STRNG && g.subtype==-1) return g; 301 return symbolic(at_discreted,g); 302 } 303 static const char _discreted_s []="discreted"; 304 static define_unary_function_eval (__discreted,&_discreted,_discreted_s); 305 define_unary_function_ptr5( at_discreted,alias_at_discreted,&__discreted,0,true); 306 #endif 307 randNorm(GIAC_CONTEXT)308 double randNorm(GIAC_CONTEXT){ 309 /* 310 double d=rand()/(rand_max2+1.0); 311 d=2*d-1; 312 identificateur x(" x"); 313 return newton(erf(x)-d,x,d); 314 */ 315 double u=giac_rand(contextptr)/(rand_max2+1.0); 316 double d=giac_rand(contextptr)/(rand_max2+1.0); 317 return std::sqrt(-2*std::log(u))*std::cos(2*M_PI*d); 318 } randnorm2(double & r1,double & r2,GIAC_CONTEXT)319 void randnorm2(double & r1,double & r2,GIAC_CONTEXT){ 320 /* 321 double d=rand()/(rand_max2+1.0); 322 d=2*d-1; 323 identificateur x(" x"); 324 return newton(erf(x)-d,x,d); 325 */ 326 for (;;){ 327 double u=giac_rand(contextptr)/(rand_max2+1.0); 328 double v=giac_rand(contextptr)/(rand_max2+1.0); 329 double w=u*u+v*v; 330 if (w>0 && w<=1){ 331 w=std::sqrt(-2*std::log(w)/w); 332 r1=u*w; 333 r2=v*w; 334 return; 335 } 336 } 337 } 338 _randNorm(const gen & args,GIAC_CONTEXT)339 gen _randNorm(const gen & args,GIAC_CONTEXT){ 340 if ( args.type==_STRNG && args.subtype==-1) return args; 341 if (args.type==_VECT && args._VECTptr->empty()) 342 return randNorm(contextptr); 343 if (args.type!=_VECT || args._VECTptr->size()!=2) 344 return gensizeerr(contextptr); 345 vecteur & v=*args._VECTptr; 346 if (v[1].type==_VECT){ 347 if (!is_squarematrix(v[1])) 348 return gendimerr(contextptr); 349 int n=int(v[1]._VECTptr->size()); 350 vecteur w(n); 351 for (int i=0;i<n;++i) 352 w[i]=randNorm(contextptr); 353 return evalf(v[0]+v[1]*w,1,contextptr); 354 } 355 return evalf(v[0]+v[1]*gen(randNorm(contextptr)),1,contextptr); 356 } 357 static const char _randNorm_s []="randNorm"; 358 static define_unary_function_eval (__randNorm,&_randNorm,_randNorm_s); 359 define_unary_function_ptr5( at_randNorm ,alias_at_randNorm,&__randNorm,0,true); 360 361 static const char _randnormald_s []="randnormald"; 362 static define_unary_function_eval (__randnormald,&_randNorm,_randnormald_s); 363 define_unary_function_ptr5( at_randnormald ,alias_at_randnormald,&__randnormald,0,true); 364 365 static const char _normalvariate_s []="normalvariate"; 366 static define_unary_function_eval (__normalvariate,&_randNorm,_normalvariate_s); 367 define_unary_function_ptr5( at_normalvariate ,alias_at_normalvariate,&__normalvariate,0,true); 368 randchisquare(int k,GIAC_CONTEXT)369 double randchisquare(int k,GIAC_CONTEXT){ 370 double res=0.0; 371 for (int i=0;i<k;++i){ 372 double u=giac_rand(contextptr)/(rand_max2+1.0); 373 double d=giac_rand(contextptr)/(rand_max2+1.0); 374 u=-2*std::log(u); 375 d=std::cos(2*M_PI*d); 376 d=d*d*u; 377 res += d; 378 } 379 return res; 380 } _randchisquare(const gen & args,GIAC_CONTEXT)381 gen _randchisquare(const gen & args,GIAC_CONTEXT){ 382 if ( args.type==_STRNG && args.subtype==-1) return args; 383 gen g(args); 384 if (!is_integral(g) || g.type!=_INT_ || g.val<=0 || g.val>1000) 385 return gensizeerr(contextptr); 386 return randchisquare(g.val,contextptr); 387 } 388 static const char _randchisquare_s []="randchisquare"; 389 static define_unary_function_eval (__randchisquare,&_randchisquare,_randchisquare_s); 390 define_unary_function_ptr5( at_randchisquare ,alias_at_randchisquare,&__randchisquare,0,true); 391 392 static const char _randchisquared_s []="randchisquared"; 393 static define_unary_function_eval (__randchisquared,&_randchisquare,_randchisquared_s); 394 define_unary_function_ptr5( at_randchisquared ,alias_at_randchisquared,&__randchisquared,0,true); 395 randstudent(int k,GIAC_CONTEXT)396 double randstudent(int k,GIAC_CONTEXT){ 397 return randNorm(contextptr)/std::sqrt(randchisquare(k,contextptr)/k); 398 } _randstudent(const gen & args,GIAC_CONTEXT)399 gen _randstudent(const gen & args,GIAC_CONTEXT){ 400 if ( args.type==_STRNG && args.subtype==-1) return args; 401 gen g(args); 402 if (!is_integral(g) || g.type!=_INT_ || g.val<=0 || g.val>1000) 403 return gensizeerr(contextptr); 404 return randstudent(g.val,contextptr); 405 } 406 static const char _randstudent_s []="randstudent"; 407 static define_unary_function_eval (__randstudent,&_randstudent,_randstudent_s); 408 define_unary_function_ptr5( at_randstudent ,alias_at_randstudent,&__randstudent,0,true); 409 static const char _randstudentd_s []="randstudentd"; 410 static define_unary_function_eval (__randstudentd,&_randstudent,_randstudentd_s); 411 define_unary_function_ptr5( at_randstudentd ,alias_at_randstudentd,&__randstudentd,0,true); 412 randfisher(int k1,int k2,GIAC_CONTEXT)413 double randfisher(int k1,int k2,GIAC_CONTEXT){ 414 return randchisquare(k1,contextptr)/k1/(randchisquare(k2,contextptr)/k2); 415 } _randfisher(const gen & args,GIAC_CONTEXT)416 gen _randfisher(const gen & args,GIAC_CONTEXT){ 417 if ( args.type==_STRNG && args.subtype==-1) return args; 418 if (args.type!=_VECT || args._VECTptr->size()!=2) 419 return gensizeerr(contextptr); 420 gen g1(args._VECTptr->front()),g2(args._VECTptr->back()); 421 if (!is_integral(g1) || g1.type!=_INT_ || g1.val<=0 || g1.val>1000 || 422 !is_integral(g2) || g2.type!=_INT_ || g2.val<=0 || g2.val>1000 423 ) 424 return gensizeerr(contextptr); 425 return randfisher(g1.val,g2.val,contextptr); 426 } 427 static const char _randfisher_s []="randfisher"; 428 static define_unary_function_eval (__randfisher,&_randfisher,_randfisher_s); 429 define_unary_function_ptr5( at_randfisher ,alias_at_randfisher,&__randfisher,0,true); 430 static const char _randfisherd_s []="randfisherd"; 431 static define_unary_function_eval (__randfisherd,&_randfisher,_randfisherd_s); 432 define_unary_function_ptr5( at_randfisherd ,alias_at_randfisherd,&__randfisherd,0,true); 433 normald(const gen & m,const gen & s,const gen & x,GIAC_CONTEXT)434 static gen normald(const gen & m,const gen & s,const gen & x,GIAC_CONTEXT){ 435 gen v(s*s); 436 return inv(sqrt(2*cst_pi*v,contextptr),contextptr)*exp(-pow(x-m,2)/(2*v),contextptr); 437 } _normald(const gen & g,GIAC_CONTEXT)438 gen _normald(const gen & g,GIAC_CONTEXT){ 439 if ( g.type==_STRNG && g.subtype==-1) return g; 440 if (g.type!=_VECT) 441 return normald(0,1,g,contextptr); 442 vecteur & v=*g._VECTptr; 443 int s=int(v.size()); 444 if (s==2) 445 return symbolic(at_normald,g); 446 if (s==3) 447 return normald(v[0],v[1],v[2],contextptr); 448 return gensizeerr(contextptr); 449 } 450 static const char _normald_s []="normald"; 451 static define_unary_function_eval (__normald,&_normald,_normald_s); 452 define_unary_function_ptr5( at_normald ,alias_at_normald,&__normald,0,true); 453 454 static const char _NORMALD_s []="NORMALD"; 455 static define_unary_function_eval (__NORMALD,&_normald,_NORMALD_s); 456 define_unary_function_ptr5( at_NORMALD ,alias_at_NORMALD,&__NORMALD,0,true); 457 unif_rand(GIAC_CONTEXT)458 double unif_rand(GIAC_CONTEXT){ 459 return giac_rand(contextptr)/(rand_max2+1.0); 460 } exp_rand(GIAC_CONTEXT)461 double exp_rand(GIAC_CONTEXT){ 462 double u=giac_rand(contextptr)/(rand_max2+1.0); 463 return -std::log(1-u); 464 } _randexp(const gen & args,GIAC_CONTEXT)465 gen _randexp(const gen & args,GIAC_CONTEXT){ 466 if ( args.type==_STRNG && args.subtype==-1) return args; 467 return gen(exp_rand(contextptr))/args; 468 } 469 static const char _randexp_s []="randexp"; 470 static define_unary_function_eval (__randexp,&_randexp,_randexp_s); 471 define_unary_function_ptr5( at_randexp ,alias_at_randexp,&__randexp,0,true); 472 473 static const char _expovariate_s []="expovariate"; 474 static define_unary_function_eval (__expovariate,&_randexp,_expovariate_s); 475 define_unary_function_ptr5( at_expovariate ,alias_at_expovariate,&__expovariate,0,true); 476 477 // Normal cumulative distribution function 478 // proba that X<x for X following a normal distrib of mean mean and dev dev 479 // arg = vector [mean,dev,x] or x alone (mean=0, dev=1) normal_cdf(const gen & g,GIAC_CONTEXT)480 static gen normal_cdf(const gen & g,GIAC_CONTEXT){ 481 if (g.type==_DOUBLE_){ 482 double x=-g._DOUBLE_val*std::sqrt(2.0)/2; 483 #if 1 484 if (x>0){ 485 // erf is odd, (erf(sqrt(2)/2*g)+1)/2=(1-erf(x))/2=erfc(x)/2 486 return .5*erfc(x,contextptr); 487 } 488 #else 489 if (x>5){ 490 double p=.3275911,a1=.254829592,a2=-.284496736,a3=1.421413741,a4=-1.453152027,a5=1.061405429; 491 double t=1.0/(1+p*x); 492 return .5*std::exp(-x*x)*t*(a1+t*(a2+t*(a3+t*(a4+t*a5)))); 493 //double p=0.47047,a1=.3480242,a2=-0.0958798,a3=.7478556; 494 //double t=1.0/(1+p*x); 495 //return .5*std::exp(-x*x)*t*(a1+t*(a2+t*a3)); 496 } 497 #endif 498 } 499 return rdiv(erf(ratnormal(plus_sqrt2_2*g,contextptr),contextptr)+plus_one,2,contextptr); 500 } _normal_cdf(const gen & g,GIAC_CONTEXT)501 gen _normal_cdf(const gen & g,GIAC_CONTEXT){ 502 if ( g.type==_STRNG && g.subtype==-1) return g; 503 if (g.type!=_VECT) 504 return normal_cdf(g,contextptr); 505 vecteur v=*g._VECTptr; 506 int s=int(v.size()); 507 if (s==2){ 508 v.insert(v.begin(),1); 509 v.insert(v.begin(),0); 510 s +=2; 511 // return normal_cdf(v[1],contextptr)-normal_cdf(v[0],contextptr); 512 } 513 if (s==3) 514 return normal_cdf((v[2]-v[0])/v[1],contextptr); 515 if (s==4){ 516 // precision: 517 if (is_strictly_greater(v[2],v[3],contextptr)) 518 return gensizeerr(contextptr); 519 if (is_strictly_greater(v[3]-v[0],v[0]-v[2],contextptr)){ 520 v[0]=-v[0]; 521 v[2]=-v[2]; 522 v[3]=-v[3]; 523 swapgen(v[2],v[3]); 524 } 525 return normal_cdf((v[3]-v[0])/v[1],contextptr)-normal_cdf((v[2]-v[0])/v[1],contextptr); 526 } 527 return gensizeerr(contextptr); 528 } 529 static const char _normal_cdf_s []="normal_cdf"; 530 static define_unary_function_eval (__normal_cdf,&_normal_cdf,_normal_cdf_s); 531 define_unary_function_ptr5( at_normal_cdf ,alias_at_normal_cdf,&__normal_cdf,0,true); 532 533 static const char _normald_cdf_s []="normald_cdf"; 534 static define_unary_function_eval (__normald_cdf,&_normal_cdf,_normald_cdf_s); 535 define_unary_function_ptr5( at_normald_cdf ,alias_at_normald_cdf,&__normald_cdf,0,true); 536 537 // returns x s.t. UTPN(x) ~ y 538 // ref Abramowitz & Stegun equation 26.2.22 utpn_initial_guess(double y)539 static double utpn_initial_guess(double y){ 540 double t=std::sqrt(-2*std::log(y)); 541 t=t-(2.30753+.27061*t)/(1+0.99229*t+0.04481*t*t); 542 return t; 543 } utpn_inverse(double y,GIAC_CONTEXT)544 static gen utpn_inverse(double y,GIAC_CONTEXT){ 545 identificateur x(" x"); 546 return newton(erf(x/std::sqrt(2.0),contextptr)+2*y-1,x,utpn_initial_guess(y),NEWTON_DEFAULT_ITERATION,1e-5,1e-12,true,1,0,1,0,.5,contextptr); 547 } normal_icdf(const gen & g_orig,GIAC_CONTEXT)548 static gen normal_icdf(const gen & g_orig,GIAC_CONTEXT){ 549 gen g=evalf_double(g_orig,1,contextptr); 550 if (g.type!=_DOUBLE_ ) 551 return symbolic(at_normal_icdf,g); 552 if (g._DOUBLE_val<0 || g._DOUBLE_val>1) 553 return gensizeerr(contextptr); 554 if (g._DOUBLE_val==0) 555 return minus_inf; 556 if (g._DOUBLE_val==1) 557 return plus_inf; 558 return utpn_inverse(1-g._DOUBLE_val,contextptr); 559 // identificateur x(" x"); 560 // plus_sqrt2*newton(erf(x)+1-2*g,x,0); 561 } _normal_icdf(const gen & g,GIAC_CONTEXT)562 gen _normal_icdf(const gen & g,GIAC_CONTEXT){ 563 if ( g.type==_STRNG && g.subtype==-1) return g; 564 if (g.type!=_VECT) 565 return normal_icdf(g,contextptr); 566 vecteur v=*g._VECTptr; 567 int s=v.size(); 568 if (s==2 && (v.back()==at_left || v.back()==at_right || v.back()==at_centre|| v.back()==at_tail)){ 569 v=makevecteur(0,1,v[0],v[1]); 570 s=4; 571 } 572 if (s<3) 573 return gensizeerr(contextptr); 574 if (s==4 && v.back()==at_centre) 575 v[2]=(1-v[2])/2; 576 if (s==4 && v.back()==at_tail) 577 v[2]=v[2]/2; 578 gen g2(normal_icdf(v[2],contextptr)); 579 gen g1=v[0]-v[1]*g2; 580 g2=v[0]+v[1]*g2; 581 if (s==4 && (v.back()==at_centre || v.back()==at_tail)) 582 return makevecteur(g2,g1); 583 if (s==4 && v.back()==at_right) 584 return g1; 585 return g2; 586 } 587 static const char _normal_icdf_s []="normal_icdf"; 588 static define_unary_function_eval (__normal_icdf,&_normal_icdf,_normal_icdf_s); 589 define_unary_function_ptr5( at_normal_icdf ,alias_at_normal_icdf,&__normal_icdf,0,true); 590 591 static const char _normald_icdf_s []="normald_icdf"; 592 static define_unary_function_eval (__normald_icdf,&_normal_icdf,_normald_icdf_s); 593 define_unary_function_ptr5( at_normald_icdf ,alias_at_normald_icdf,&__normald_icdf,0,true); 594 apply3rd(const gen & e1,const gen & e2,const gen & e3,const context * contextptr,gen (* f)(const gen &,const gen &,const gen &,const context *))595 gen apply3rd(const gen & e1, const gen & e2,const gen & e3,const context * contextptr, gen (* f) (const gen &, const gen &,const gen &,const context *) ){ 596 if (e3.type!=_VECT) 597 return f(e1,e2,e3,contextptr); 598 const_iterateur it=e3._VECTptr->begin(),itend=e3._VECTptr->end(); 599 vecteur v; 600 v.reserve(itend-it); 601 for (;it!=itend;++it){ 602 gen tmp=f(e1,e2,*it,contextptr); 603 if (is_undef(tmp)) 604 return gen2vecteur(tmp); 605 v.push_back(tmp); 606 } 607 return gen(v,e3.subtype); 608 } 609 binomial(const gen & N,const gen & K,const gen & P,GIAC_CONTEXT)610 gen binomial(const gen & N,const gen & K,const gen & P,GIAC_CONTEXT){ 611 gen n(N),k(K),p(P); 612 is_integral(n); is_integral(k); is_integral(p); 613 if (p.type==_VECT) 614 return apply3rd(n,k,p,contextptr,binomial); 615 if ( (is_zero(p) && is_zero(k)) || (is_one(p) && n==k)) 616 return 1; 617 if (is_strictly_positive(-n,contextptr)) 618 return _negbinomial(makesequence(-n,k,p),contextptr); 619 if (is_strictly_positive(k,contextptr) && is_strictly_greater(1,k,contextptr)){ 620 if (is_strictly_positive(p,contextptr) && is_strictly_greater(1,p,contextptr)) 621 return gensizeerr(contextptr); 622 return binomial(n,p,k,contextptr); 623 } 624 if (k.type==_DOUBLE_ || k.type==_FLOAT_ || k.type==_FRAC){ 625 if (p.type==_DOUBLE_ || p.type==_FLOAT_ || p.type==_FRAC) 626 return gensizeerr(contextptr); 627 return binomial(n,p,k,contextptr); 628 } 629 if (p.type==_DOUBLE_ || p.type==_FLOAT_){ 630 #if 0 //def VISUALC 631 gen nd=evalf2bcd(n,1,contextptr); 632 gen kd=evalf2bcd(k,1,contextptr); 633 gen pd=evalf2bcd(p,1,contextptr); 634 if (nd.type==_FLOAT_ && kd.type==_FLOAT_ && pd.type==_FLOAT_){ 635 return Gamma(nd+1,contextptr)/Gamma(kd+1,contextptr)/Gamma(nd-kd+1,contextptr)*pow(pd,kd,contextptr)*pow(1-pd,nd-kd,contextptr); 636 } 637 #else 638 gen nd=evalf_double(n,1,contextptr); 639 gen kd=evalf_double(k,1,contextptr); 640 gen pd=evalf_double(p,1,contextptr); 641 if (nd.type==_DOUBLE_ && kd.type==_DOUBLE_ && pd.type==_DOUBLE_){ 642 double ndd=nd._DOUBLE_val,kdd=kd._DOUBLE_val,pdd=pd._DOUBLE_val,nk=ndd-kdd; 643 return std::exp(lngamma(ndd+1)-lngamma(kdd+1)-lngamma(nk+1)+kdd*std::log(pdd)+nk*std::log(1-pdd)); 644 } 645 #endif 646 } 647 if (!is_positive(p,contextptr) || !is_greater(1,p,contextptr)){ 648 if (abs_calc_mode(contextptr)==38) 649 return gensizeerr(contextptr); 650 if (calc_mode(contextptr)!=1) 651 *logptr(contextptr) << "Assuming probability=" << p << "\n"; 652 } 653 return comb(n,k,contextptr)*pow(p,k,contextptr)*pow(1-p,n-k,contextptr); 654 } _binomial(const gen & g,GIAC_CONTEXT)655 gen _binomial(const gen & g,GIAC_CONTEXT){ 656 if ( g.type==_STRNG && g.subtype==-1) return g; 657 if (g.type!=_VECT) 658 return gensizeerr(contextptr); 659 vecteur & v=*g._VECTptr; 660 int s=int(v.size()); 661 if (s==2){ 662 if (is_strictly_positive(v[1],contextptr) && is_strictly_greater(1,v[1],contextptr)) 663 return symbolic(at_binomial,g); // inert form for binomial pseudo-random generation 664 gen v0=evalf_double(v[0],1,contextptr),v1=evalf_double(v[1],1,contextptr); 665 if (v0.type!=_DOUBLE_ || v1.type!=_DOUBLE_) 666 return _factorial(v[0],contextptr)/(_factorial(v[1],contextptr)*_factorial(v[0]-v[1],contextptr)); 667 return comb(v[0],v[1],contextptr); 668 } 669 if (s==3){ 670 if (0 && calc_mode(contextptr)==1) 671 return binomial(v[0],v[2],v[1],contextptr); 672 return binomial(v[0],v[1],v[2],contextptr); 673 } 674 return gensizeerr(contextptr); 675 } 676 static const char _binomial_s []="binomial"; 677 static define_unary_function_eval (__binomial,&_binomial,_binomial_s); 678 define_unary_function_ptr5( at_binomial ,alias_at_binomial,&__binomial,0,true); 679 680 static const char _BINOMIAL_s []="BINOMIAL"; 681 static define_unary_function_eval (__BINOMIAL,&_binomial,_BINOMIAL_s); 682 define_unary_function_ptr5( at_BINOMIAL ,alias_at_BINOMIAL,&__BINOMIAL,0,true); 683 _multinomial(const gen & g,GIAC_CONTEXT)684 gen _multinomial(const gen & g,GIAC_CONTEXT){ 685 if ( g.type==_STRNG && g.subtype==-1) return g; 686 if (g.type!=_VECT) 687 return gensizeerr(contextptr); 688 vecteur & v=*g._VECTptr; 689 int s=int(v.size()); 690 if (s==3){ 691 gen n=v[0],K=v[1],P=v[2]; 692 if (!is_zero(1-_plus(P,contextptr),contextptr)) 693 swapgen(K,P); 694 if (_plus(K,contextptr)!=n || K.type!=_VECT || P.type!=_VECT || K._VECTptr->size()!=P._VECTptr->size()) 695 return gensizeerr(contextptr); 696 vecteur k=*K._VECTptr,p=*P._VECTptr; 697 unsigned s=unsigned(k.size()); 698 gen res=_factorial(n,contextptr); 699 for (unsigned i=0;i<s;++i){ 700 res=res/_factorial(k[i],contextptr); 701 } 702 for (unsigned i=0;i<s;++i){ 703 res=res*pow(p[i],k[i],contextptr); 704 } 705 return res; 706 } 707 return gensizeerr(contextptr); 708 } 709 static const char _multinomial_s []="multinomial"; 710 static define_unary_function_eval (__multinomial,&_multinomial,_multinomial_s); 711 define_unary_function_ptr5( at_multinomial ,alias_at_multinomial,&__multinomial,0,true); 712 _randmultinomial(const gen & args,GIAC_CONTEXT)713 gen _randmultinomial(const gen & args,GIAC_CONTEXT){ 714 if ( args.type==_STRNG && args.subtype==-1) return args; 715 if (args.type!=_VECT || args._VECTptr->empty()) 716 return gensizeerr(contextptr); 717 gen g1(args._VECTptr->front()); 718 if (args._VECTptr->size()==2 && g1.type==_VECT){ 719 gen g2(args._VECTptr->back()); 720 if (g2.type!=_VECT || g2._VECTptr->size()!=g1._VECTptr->size() || !is_zero(1-_sum(g1,contextptr)) ) 721 return gensizeerr(contextptr); 722 double u=giac_rand(contextptr)/(rand_max2+1.0); 723 gen somme=0; 724 for (unsigned i=0;i<g1._VECTptr->size();++i){ 725 somme += g1[i]; 726 if (is_greater(somme,u,contextptr)) 727 return g2[i]; 728 } 729 return undef; 730 } 731 if (!is_zero(1-_sum(args,contextptr))) 732 return gensizeerr(contextptr); 733 double u=giac_rand(contextptr)/(rand_max2+1.0); 734 gen somme=0; 735 for (unsigned i=0;i<args._VECTptr->size();++i){ 736 somme += (*args._VECTptr)[i]; 737 if (is_greater(somme,u,contextptr)) 738 return int(i); 739 } 740 return undef; 741 } 742 static const char _randmultinomial_s []="randmultinomial"; 743 static define_unary_function_eval (__randmultinomial,&_randmultinomial,_randmultinomial_s); 744 define_unary_function_ptr5( at_randmultinomial ,alias_at_randmultinomial,&__randmultinomial,0,true); 745 _negbinomial(const gen & g,GIAC_CONTEXT)746 gen _negbinomial(const gen & g,GIAC_CONTEXT){ 747 if ( g.type==_STRNG && g.subtype==-1) return g; 748 if (g.type!=_VECT) 749 return gensizeerr(contextptr); 750 vecteur & v=*g._VECTptr; 751 int s=int(v.size()); 752 if (s==2) return symbolic(at_negbinomial,g); // for random 753 if (s==3){ 754 gen r=v[0],p=v[1],k=v[2]; 755 gen tmp=evalf_double(k,1,contextptr); 756 if (tmp.type==_DOUBLE_ && tmp._DOUBLE_val<1 && tmp._DOUBLE_val>0) 757 swapgen(k,p); 758 if (is_integral(p)) 759 swapgen(k,p); 760 if (is_zero(k)) 761 return pow(p,r,contextptr); 762 return binomial(r+k-1,r,p,contextptr)*(1-p)*r/k; 763 } 764 return gensizeerr(contextptr); 765 } 766 static const char _negbinomial_s []="negbinomial"; 767 static define_unary_function_eval (__negbinomial,&_negbinomial,_negbinomial_s); 768 define_unary_function_ptr5( at_negbinomial ,alias_at_negbinomial,&__negbinomial,0,true); 769 _negbinomial_cdf(const gen & g,GIAC_CONTEXT)770 gen _negbinomial_cdf(const gen & g,GIAC_CONTEXT){ 771 if ( g.type==_STRNG && g.subtype==-1) return g; 772 if (g.type!=_VECT) 773 return gensizeerr(contextptr); 774 vecteur & v=*g._VECTptr; 775 int s=int(v.size()); 776 if (s==3){ 777 gen n=v[0],p=v[1],k=v[2]; 778 return _Beta(makesequence(n,k+1,p,1),contextptr); 779 } 780 if (s==4){ 781 gen n=v[0],p=v[1],k1=v[2],k2=v[3]; 782 return _Beta(makesequence(n,k2+1,p,1),contextptr)-_Beta(makesequence(n,k1+1,p,1),contextptr); 783 } 784 return gensizeerr(contextptr); 785 } 786 static const char _negbinomial_cdf_s []="negbinomial_cdf"; 787 static define_unary_function_eval (__negbinomial_cdf,&_negbinomial_cdf,_negbinomial_cdf_s); 788 define_unary_function_ptr5( at_negbinomial_cdf ,alias_at_negbinomial_cdf,&__negbinomial_cdf,0,true); 789 _negbinomial_icdf(const gen & g,GIAC_CONTEXT)790 gen _negbinomial_icdf(const gen & g,GIAC_CONTEXT){ 791 if ( g.type==_STRNG && g.subtype==-1) return g; 792 if (g.type!=_VECT) 793 return gensizeerr(contextptr); 794 vecteur & v=*g._VECTptr; 795 int s=int(v.size()); 796 if (s==3){ 797 gen R=v[0],P=evalf_double(v[1],1,contextptr),T=v[2]; 798 if (!is_integral(R) || R.val<=0 || P._DOUBLE_val<=0 || P._DOUBLE_val>=1) 799 return gensizeerr(contextptr); 800 int r=R.val; 801 long_double p=P._DOUBLE_val,t=T._DOUBLE_val; 802 if (t<=0) 803 return 0; 804 if (t>=1) 805 return 1; 806 long_double cumul=std::pow(p,r),current=cumul; 807 if (cumul==0){ 808 *logptr(contextptr) << gettext("Underflow") <<"\n"; 809 return undef; 810 } 811 // negbinomial(r,p,k+1)/negbinomial(r,p,k)))=(1-p)*(k+r)/(k+1) 812 for (int k=0;;){ 813 if (cumul>=t) 814 return k; 815 current=current*(k+r)*(1-p)/(k+1); 816 if (cumul==cumul+current) 817 return k; 818 ++k; 819 cumul += current; 820 } 821 } 822 return gensizeerr(contextptr); 823 } 824 static const char _negbinomial_icdf_s []="negbinomial_icdf"; 825 static define_unary_function_eval (__negbinomial_icdf,&_negbinomial_icdf,_negbinomial_icdf_s); 826 define_unary_function_ptr5( at_negbinomial_icdf ,alias_at_negbinomial_icdf,&__negbinomial_icdf,0,true); 827 binomial_cdf(const gen & n,const gen & p,const gen & x0,const gen & x,GIAC_CONTEXT)828 gen binomial_cdf(const gen & n,const gen &p,const gen & x0,const gen & x,GIAC_CONTEXT){ 829 gen fx=_floor(x,contextptr),fx0=_ceil(x0,contextptr); 830 if (fx.type==_FLOAT_) 831 fx=get_int(fx._FLOAT_val); 832 if (fx0.type==_FLOAT_) 833 fx0=get_int(fx0._FLOAT_val); 834 if (fx.type!=_INT_ || fx.val<0 || fx0.type!=_INT_ || fx0.val<0) 835 return gensizeerr(contextptr); 836 if (fx0.val>fx.val) 837 return 0; 838 gen pd=p; 839 if (pd.type==_FLOAT_) pd=evalf_double(p,1,contextptr); 840 if (n.type==_INT_ && pd.type==_DOUBLE_ && (fx.val-fx0.val)>100){ 841 // improve if n large using Abramowitz-Stegun p. 944-945 842 // [Not very useful, perhaps for much larger values?] 843 // sum(binomial(n,p,k),k=a..n)=Ip(a,n-a+1) 844 // hence binomial_cdf(n,p,fx0,fx)=Ip(fx0,n-fx0+1)-Ip(fx+1,n-(fx+1)+1) 845 // If more than 100 terms to compute use 846 // I_p(a,b)=p^a*(1-p)^b/a/B(a,b)*[1+sum(B(a+1,n+1)/B(a+b,n+1)*p^(n+1))] 847 // B(a,b)=Gamma(a)*Gamma(b)/Gamma(a+b) 848 double p=pd._DOUBLE_val; 849 gen res; 850 if (fx0.val<=0){ 851 if (fx.val>=n.val) 852 return 1; 853 res=1-incomplete_beta(fx.val+1,n.val-fx.val,p); 854 } 855 else { 856 if (fx.val>=n.val) 857 res=incomplete_beta(fx0.val,n.val-fx0.val+1,p); 858 else 859 res=incomplete_beta(fx0.val,n.val-fx0.val+1,p)-incomplete_beta(fx.val+1,n.val-fx.val,p); 860 } 861 if (!is_undef(res)) 862 return res; 863 // Other formula 864 // n=a+b-1 865 // I_{1-p}(b,a)=sum(i=0,a-1,comb(a+b-1,i)*p^i*(1-p)^(a+b-1-i)=binomial_cdf(a+b-1,p,0,a-1) 866 // I_p(a,b)=1-I_{1-p}(b,a) 867 // I_{1-p}(b,a)= Gamma(b,y)/Gamma(b)-1/24/N^2*y^b*exp(-y)/(b-2)!*(b+1+y) + 868 // +1/5760/N^4*y^b*exp(-y)/(b-2)!*[(b-3)*(b-2)*(5*b+7)*(b+1+y)-(5b-7)*(b+3+y)*y^2] 869 // N=a+b/2-1/2, y=-N*ln(p), a>>b 870 // Gamma(b,y) = y^(b-1)*exp(-y)*(1+(b-1)/y+(b-1)*(b-2)/y^2+...) 871 } 872 gen last=binomial(n,fx0.val,pd,contextptr); 873 if (last.type==_FLOAT_) 874 last=evalf_double(last,1,contextptr); 875 gen res=last; 876 if (n.type==_INT_ && pd.type==_DOUBLE_ && last.type==_DOUBLE_ && last._DOUBLE_val!=0){ 877 double tmp=pd._DOUBLE_val/(1-pd._DOUBLE_val); 878 double lastd=last._DOUBLE_val,resd=res._DOUBLE_val; 879 for (int i=fx0.val+1;i<=fx.val;++i){ 880 if (i%25 == 0) 881 lastd=evalf_double(binomial(n,i,pd,contextptr),1,contextptr)._DOUBLE_val; // avoid loss of precision 882 else 883 lastd *= (n.val-i+1)*tmp/i; 884 resd += lastd; 885 } 886 return resd; 887 } 888 if (n.type==_INT_ && pd.type==_FRAC && last.type==_FRAC){ 889 gen tmp=pd/(1-pd); 890 for (int i=fx0.val+1;i<=fx.val;++i){ 891 last = last*gen((n.val-i+1)*tmp/i); 892 res += last; 893 } 894 } 895 else { 896 for (int i=fx0.val+1;i<=fx.val;++i){ 897 res += binomial(n,i,pd,contextptr); 898 } 899 } 900 return res; 901 } _binomial_cdf(const gen & g,GIAC_CONTEXT)902 gen _binomial_cdf(const gen & g,GIAC_CONTEXT){ 903 if ( g.type==_STRNG && g.subtype==-1) return g; 904 if (g.type!=_VECT) 905 return gensizeerr(contextptr); 906 vecteur & v=*g._VECTptr; 907 int s=int(v.size()); 908 if (s==3){ 909 if (v[2].type==_IDNT) 910 return symbolic(at_binomial_cdf,makesequence(v[0],v[1],v[2])); 911 return binomial_cdf(v[0],v[1],0,v[2],contextptr); 912 } 913 if (s==4) 914 return binomial_cdf(v[0],v[1],v[2],v[3],contextptr); 915 return gensizeerr(contextptr); 916 } 917 static const char _binomial_cdf_s []="binomial_cdf"; 918 static define_unary_function_eval (__binomial_cdf,&_binomial_cdf,_binomial_cdf_s); 919 define_unary_function_ptr5( at_binomial_cdf ,alias_at_binomial_cdf,&__binomial_cdf,0,true); 920 binomial_icdf(const gen & n0,const gen & p0,const gen & x_orig,GIAC_CONTEXT)921 gen binomial_icdf(const gen & n0,const gen &p0,const gen & x_orig,GIAC_CONTEXT){ 922 gen x=evalf_double(x_orig,1,contextptr); 923 gen p=evalf_double(p0,1,contextptr); 924 gen n=_floor(n0,contextptr); 925 if (n.type==_FLOAT_) 926 n=get_int(n._FLOAT_val); 927 if (!is_zero(n-n0)) 928 return gensizeerr(contextptr); 929 if (x._DOUBLE_val==0) 930 return zero; 931 if (x._DOUBLE_val==1) 932 return n; 933 if (is_strictly_greater(p,1,contextptr) || is_strictly_greater(0,p,contextptr)) 934 return gensizeerr(contextptr); 935 if (n.type!=_INT_ || p.type!=_DOUBLE_ || x.type!=_DOUBLE_ || x._DOUBLE_val<0 || x._DOUBLE_val>1 ) 936 return symbolic(at_binomial_icdf,makesequence(n,p,x)); 937 int N=n.val; 938 long_double P=p._DOUBLE_val; 939 if (N*P>=30 && N*(1-P)>=30){ 940 // use approximation by normal law as a starting point 941 gen g=_floor(_normal_icdf(makesequence(n*p,sqrt(n*p*(1-p),contextptr),x),contextptr),contextptr); 942 int G=g.val; 943 gen cdf=evalf_double(_binomial_cdf(makesequence(n,p,g),contextptr),1,contextptr); 944 long_double CDF=cdf._DOUBLE_val; 945 long_double T=x._DOUBLE_val; 946 long_double current=std::exp(lngamma(N+1)-lngamma(G+1)-lngamma(N-G+1)+G*std::log(P)+(N-G)*std::log(1-P)); 947 long_double P1=P/(1-P); 948 // increase G until CDF>=T 949 for (;T>CDF;){ 950 ++G; 951 current *= (N-G+1)*P1/G; 952 CDF += current; 953 } 954 if (G!=g.val) 955 return G; 956 // decrease G until CDF<T 957 for (;T<=CDF;){ 958 CDF -= current; 959 current /= (N-G+1)*P1/G; 960 --G; 961 } 962 return G+1; 963 } 964 gen last=pow(1-p0,n,contextptr); 965 gen b=last; 966 int k=1; 967 if (last.type==_FLOAT_) 968 last=evalf_double(last,1,contextptr); 969 if (last.type==_DOUBLE_){ 970 double tmp=(p/(1-p))._DOUBLE_val,lastd=last._DOUBLE_val,bd=lastd; 971 for (;k<=n.val;++k){ 972 if (x._DOUBLE_val<=bd) 973 return k-1; 974 if (k%25 == 0) 975 lastd=evalf_double(binomial(n,k,p,contextptr),1,contextptr)._DOUBLE_val; // avoid loss of precision 976 else 977 lastd *= (n.val-k+1)*tmp/k; 978 bd += lastd; 979 } 980 return n; 981 } 982 #if 1 983 gen tmp=p0/(1-p0); 984 for (;k<=n.val;++k){ 985 if (!ck_is_strictly_greater(x_orig,b,contextptr)) 986 return k-1; 987 last = last * ((n.val-k+1)*tmp)/k; 988 b += last; 989 } 990 #else 991 for (;k<=n.val;++k){ 992 if (!ck_is_strictly_greater(x,b,contextptr)) 993 return k-1; 994 b=b+binomial(n,k,p,contextptr); 995 } 996 #endif 997 return n; 998 } _binomial_icdf(const gen & g,GIAC_CONTEXT)999 gen _binomial_icdf(const gen & g,GIAC_CONTEXT){ 1000 if ( g.type==_STRNG && g.subtype==-1) return g; 1001 if (g.type!=_VECT) 1002 return gensizeerr(contextptr); 1003 vecteur & v=*g._VECTptr; 1004 int s=int(v.size()); 1005 if (s==3) 1006 return binomial_icdf(v[0],v[1],v[2],contextptr); 1007 if (s==4) 1008 return binomial_icdf(v[0],v[1],v[3],contextptr)-binomial_icdf(v[0],v[1],v[2],contextptr); 1009 return gensizeerr(contextptr); 1010 } 1011 static const char _binomial_icdf_s []="binomial_icdf"; 1012 static define_unary_function_eval (__binomial_icdf,&_binomial_icdf,_binomial_icdf_s); 1013 define_unary_function_ptr5( at_binomial_icdf ,alias_at_binomial_icdf,&__binomial_icdf,0,true); 1014 randbinomial(int n,double P,GIAC_CONTEXT)1015 gen randbinomial(int n,double P,GIAC_CONTEXT){ 1016 if (P<=0) 1017 return 0; 1018 if (P>=1) 1019 return n; 1020 if (n>1000) 1021 return binomial_icdf(n,P,double(giac_rand(contextptr))/rand_max2,contextptr); 1022 int ok=0; 1023 P*=rand_max2; 1024 for (int i=0;i<n;++i){ 1025 if (giac_rand(contextptr)<=P) 1026 ok++; 1027 } 1028 return ok; 1029 } 1030 // randbinomial(n,p) returns k in [0..n] with proba binomial(n,k,p) _randbinomial(const gen & g,GIAC_CONTEXT)1031 gen _randbinomial(const gen & g,GIAC_CONTEXT){ 1032 if ( g.type==_STRNG && g.subtype==-1) return g; 1033 if (g.type!=_VECT) 1034 return gensizeerr(contextptr); 1035 vecteur & v=*g._VECTptr; 1036 int s=int(v.size()); 1037 if (s<2 1038 || s>2 // 3 1039 ) 1040 return gensizeerr(contextptr); 1041 gen n=v[0]; 1042 gen p=v[1]; 1043 int k=1; 1044 if (s==3){ 1045 gen K=v[2]; 1046 if (!is_integral(K) || K.type!=_INT_) 1047 return gensizeerr(contextptr); 1048 k=K.val; 1049 } 1050 if (!is_integral(n) || n.type!=_INT_ || n.val<=0 || ck_is_strictly_greater(0,p,contextptr) || ck_is_strictly_greater(p,1,contextptr)) 1051 return gensizeerr(contextptr); 1052 p=evalf_double(p,1,contextptr); 1053 return randbinomial(n.val,p._DOUBLE_val,contextptr); 1054 } 1055 static const char _randbinomial_s []="randbinomial"; 1056 static define_unary_function_eval (__randbinomial,&_randbinomial,_randbinomial_s); 1057 define_unary_function_ptr5( at_randbinomial ,alias_at_randbinomial,&__randbinomial,0,true); 1058 poisson(const gen & m,const gen & k,GIAC_CONTEXT)1059 gen poisson(const gen & m,const gen & k,GIAC_CONTEXT){ 1060 if (k.type==_VECT) 1061 return apply2nd(m,k,contextptr,poisson); 1062 gen M=evalf_double(m,1,contextptr); 1063 if (M.type==_DOUBLE_){ 1064 gen K=evalf_double(k,1,contextptr); 1065 if (K.type==_DOUBLE_) 1066 return std::exp(-M._DOUBLE_val + K._DOUBLE_val*std::log(M._DOUBLE_val)-lngamma(K._DOUBLE_val+1)); 1067 } 1068 return exp(-m,contextptr)*pow(m,k,contextptr)/_factorial(k,contextptr); 1069 } _poisson(const gen & g,GIAC_CONTEXT)1070 gen _poisson(const gen & g,GIAC_CONTEXT){ 1071 if ( g.type==_STRNG && g.subtype==-1) return g; 1072 if (g.type!=_VECT) 1073 return symbolic(at_poisson,g); 1074 vecteur & v=*g._VECTptr; 1075 int s=int(v.size()); 1076 if (s==2) 1077 return poisson(v[0],v[1],contextptr); 1078 return gensizeerr(contextptr); 1079 } 1080 static const char _poisson_s []="poisson"; 1081 static define_unary_function_eval (__poisson,&_poisson,_poisson_s); 1082 define_unary_function_ptr5( at_poisson ,alias_at_poisson,&__poisson,0,true); 1083 1084 static const char _POISSON_s []="POISSON"; 1085 static define_unary_function_eval (__POISSON,&_poisson,_POISSON_s); 1086 define_unary_function_ptr5( at_POISSON ,alias_at_POISSON,&__POISSON,0,true); 1087 1088 // exp(-lambda)*sum(lambda^k/k!,k=0..x) 1089 // or 1-exp(-lambda)*sum(lambda^k/k!,k=x+1..inf) poisson_cdf(double lambda,double x)1090 double poisson_cdf(double lambda,double x){ 1091 long_double N=lambda; 1092 long_double res=0,prod=1; 1093 int fx=int(std::floor(x)); 1094 if (fx>=lambda){ 1095 for (int i=fx+1;prod>1e-17;){ 1096 res += prod; 1097 prod *= N; 1098 ++i; 1099 prod /= long_double(i); 1100 } 1101 res *= std::exp(-N+(fx+1)*std::log(N)-lngamma(fx+2.)); 1102 return 1-res; 1103 } 1104 #if 1 1105 for (int i=fx;i>=0 && prod>1e-17;--i){ 1106 res += prod; 1107 prod /= N; 1108 prod *= long_double(i); 1109 } 1110 res *= std::exp(-N+fx*std::log(N)-lngamma(fx+1.)); 1111 return res; 1112 #else 1113 for (int i=0;i<=fx;){ 1114 res += prod; 1115 prod *= N; 1116 ++i; 1117 prod /= long_double(i); 1118 } 1119 res *= std::exp(-N); 1120 return res; 1121 #endif 1122 } poisson_cdf(const gen & lambda_,const gen & x,GIAC_CONTEXT)1123 gen poisson_cdf(const gen & lambda_,const gen & x,GIAC_CONTEXT){ 1124 gen fx=_floor(x,contextptr); 1125 gen lambda=evalf_double(lambda_,1,contextptr); 1126 if (fx.type==_INT_ && fx.val>=0 && lambda.type==_DOUBLE_) 1127 return poisson_cdf(lambda._DOUBLE_val,fx.val); 1128 if (is_zero(fx-x)) 1129 return _upper_incomplete_gamma(makesequence(x+1,lambda,1),contextptr); 1130 else 1131 return _upper_incomplete_gamma(makesequence(evalf(fx,1,contextptr),lambda,1),contextptr); 1132 #if 0 1133 gen res=0; 1134 for (int i=0;i<=fx.val;++i){ 1135 res +=poisson(lambda,i,contextptr); 1136 } 1137 return res; 1138 #endif 1139 //identificateur k(" k"); 1140 //return sum(poisson(n,k,contextptr),k,0,_floor(x,contextptr),contextptr); 1141 } _poisson_cdf(const gen & g,GIAC_CONTEXT)1142 gen _poisson_cdf(const gen & g,GIAC_CONTEXT){ 1143 if ( g.type==_STRNG && g.subtype==-1) return g; 1144 if (g.type!=_VECT) 1145 return gensizeerr(contextptr); 1146 vecteur & v=*g._VECTptr; 1147 int s=int(v.size()); 1148 if (s==2) 1149 return poisson_cdf(v[0],v[1],contextptr); 1150 if (s==3) 1151 return poisson_cdf(v[0],v[2],contextptr)-poisson_cdf(v[0],v[1]-1,contextptr); 1152 return gensizeerr(contextptr); 1153 } 1154 static const char _poisson_cdf_s []="poisson_cdf"; 1155 static define_unary_function_eval (__poisson_cdf,&_poisson_cdf,_poisson_cdf_s); 1156 define_unary_function_ptr5( at_poisson_cdf ,alias_at_poisson_cdf,&__poisson_cdf,0,true); 1157 1158 // randpoisson(lambda) returns k>0 with proba poisson(lambda,k) randpoisson(double lambda,GIAC_CONTEXT)1159 gen randpoisson(double lambda,GIAC_CONTEXT){ 1160 if (lambda>700) 1161 return poisson_icdf(lambda,double(giac_rand(contextptr))/rand_max2,contextptr); 1162 int k=0; 1163 if (lambda<200){ 1164 double seuil=std::exp(-lambda); 1165 double res=1.0; 1166 for (;;++k){ 1167 res *= giac_rand(contextptr)/(rand_max2+1.0); 1168 if (res<seuil) 1169 return k; 1170 } 1171 } 1172 double res=0.0; 1173 for (;;++k){ 1174 double u = giac_rand(contextptr)/(rand_max2+1.0); 1175 res += -std::log(1-u)/lambda; 1176 if (res>=1.0) 1177 return k; 1178 } 1179 } _randpoisson(const gen & g,GIAC_CONTEXT)1180 gen _randpoisson(const gen & g,GIAC_CONTEXT){ 1181 if ( g.type==_STRNG && g.subtype==-1) return g; 1182 gen G=evalf_double(g,1,contextptr); 1183 if (G.type!=_DOUBLE_) 1184 return gensizeerr(contextptr); 1185 double lambda=G._DOUBLE_val; 1186 if (lambda<=0) 1187 return gensizeerr(contextptr); 1188 return randpoisson(lambda,contextptr); 1189 } 1190 static const char _randpoisson_s []="randpoisson"; 1191 static define_unary_function_eval (__randpoisson,&_randpoisson,_randpoisson_s); 1192 define_unary_function_ptr5( at_randpoisson ,alias_at_randpoisson,&__randpoisson,0,true); 1193 poisson_icdf(double m,double t,GIAC_CONTEXT)1194 gen poisson_icdf(double m,double t,GIAC_CONTEXT){ 1195 if (t==0) 1196 return zero; 1197 if (t==1) 1198 return plus_inf; 1199 #if 1 1200 if (m>90){ 1201 // 170.! =7e306 we must insure that the naive definition does not return >170 1202 // hence the test since poisson_cdf(90.,170.)=1.0 to double precision 1203 // approximation using normal_icdf 1204 gen g=_ceil(_normal_icdf(makesequence(m,sqrt(m,contextptr),t),contextptr),contextptr); 1205 if (is_undef(g)) 1206 return gensizeerr("Underflow"); 1207 int G=g.val; 1208 // check that poisson_cdf(m,g)>=t, if not increase g 1209 gen pg=evalf_double(_poisson_cdf(makesequence(m,g),contextptr),1,contextptr); 1210 long_double CDF=pg._DOUBLE_val; 1211 long_double M=m; 1212 long_double current= std::exp(-M+G*std::log(M)-lngamma(G+1)); 1213 long_double T=t; 1214 for (;T>CDF;){ 1215 ++G; 1216 current *= M/G; 1217 CDF += current; // std::exp(-m._DOUBLE_val+G*std::log(m._DOUBLE_val)-lngamma(G+1)); 1218 } 1219 if (G!=g.val) 1220 return G; 1221 // decrease G until cdf<t 1222 for (;T<=CDF;){ 1223 CDF -= current; // pg -= std::exp(-m._DOUBLE_val+G*std::log(m._DOUBLE_val)-lngamma(G+1)); 1224 current *= G/M; 1225 --G; 1226 } 1227 return G+1; 1228 } 1229 long_double M=m; 1230 int k=0; 1231 long_double T=t*std::exp(M),B=0,prod=1; 1232 for (;;){ 1233 B += prod; 1234 if (B>=T) 1235 return k; 1236 ++k; 1237 prod *= M; 1238 prod /= k; 1239 } 1240 #else 1241 if (m>300) 1242 return gensizeerr(gettext("Overflow")); 1243 gen b; 1244 for (;;++k){ 1245 b=b+poisson(m,k,contextptr); 1246 if (!ck_is_strictly_greater(t,b,contextptr)) 1247 return k; 1248 } 1249 return t; 1250 #endif 1251 } poisson_icdf(const gen & m_orig,const gen & t_orig,GIAC_CONTEXT)1252 gen poisson_icdf(const gen & m_orig,const gen & t_orig,GIAC_CONTEXT){ 1253 gen t=evalf_double(t_orig,1,contextptr); 1254 gen m=evalf_double(m_orig,1,contextptr); 1255 if (t.type!=_DOUBLE_ || t._DOUBLE_val<0 || t._DOUBLE_val>1) 1256 return gensizeerr(contextptr); 1257 if (m.type!=_DOUBLE_ ) 1258 return symbolic(at_poisson_icdf,makesequence(m,t)); 1259 return poisson_icdf(m._DOUBLE_val,t._DOUBLE_val,contextptr); 1260 } 1261 _poisson_icdf(const gen & g,GIAC_CONTEXT)1262 gen _poisson_icdf(const gen & g,GIAC_CONTEXT){ 1263 if ( g.type==_STRNG && g.subtype==-1) return g; 1264 if (g.type!=_VECT) 1265 return gensizeerr(contextptr); 1266 vecteur & v=*g._VECTptr; 1267 int s=int(v.size()); 1268 if (s==2) 1269 return poisson_icdf(v[0],v[1],contextptr); 1270 if (s==3) 1271 return poisson_icdf(v[0],v[2],contextptr)-poisson_icdf(v[0],v[1],contextptr); 1272 return gensizeerr(contextptr); 1273 } 1274 static const char _poisson_icdf_s []="poisson_icdf"; 1275 static define_unary_function_eval (__poisson_icdf,&_poisson_icdf,_poisson_icdf_s); 1276 define_unary_function_ptr5( at_poisson_icdf ,alias_at_poisson_icdf,&__poisson_icdf,0,true); 1277 student(const gen & n0,const gen & x,GIAC_CONTEXT)1278 gen student(const gen & n0,const gen & x,GIAC_CONTEXT){ 1279 if (x.type==_VECT) 1280 return apply2nd(n0,x,contextptr,student); 1281 gen n(n0); 1282 if (!is_integral(n) || n.val<1) 1283 return gensizeerr(contextptr); 1284 return Gamma(rdiv(n+1,2,contextptr),contextptr)/Gamma(rdiv(n,2,contextptr),contextptr)/sqrt(n*cst_pi,contextptr)*pow((1+pow(x,2)/n),-rdiv(n+1,2,contextptr),contextptr); 1285 } _student(const gen & g,GIAC_CONTEXT)1286 gen _student(const gen & g,GIAC_CONTEXT){ 1287 if ( g.type==_STRNG && g.subtype==-1) return g; 1288 if (g.type!=_VECT){ 1289 if (abs_calc_mode(contextptr)==38) 1290 return symbolic(at_student,g); 1291 return symbolic(at_studentd,g); 1292 } 1293 vecteur v=*g._VECTptr; 1294 int s=int(v.size()); 1295 if (s==2){ 1296 if (v[1].type==_DOUBLE_ || v[1].type==_FLOAT_) 1297 return evalf(student(v[0],v[1],contextptr),1,contextptr); 1298 return student(v[0],v[1],contextptr); 1299 } 1300 return gensizeerr(contextptr); 1301 } 1302 static const char _student_s []="student"; 1303 static define_unary_function_eval (__student,&_student,_student_s); 1304 define_unary_function_ptr5( at_student ,alias_at_student,&__student,0,true); 1305 static const char _studentd_s []="studentd"; 1306 static define_unary_function_eval (__studentd,&_student,_studentd_s); 1307 define_unary_function_ptr5( at_studentd ,alias_at_studentd,&__studentd,0,true); 1308 1309 /* statically handled by derive.cc 1310 static gen d2_UTPT(const gen & e,GIAC_CONTEXT ){ 1311 return -_student(e,contextptr); 1312 } 1313 static const partial_derivative_onearg D_at_UTPT(d2_UTPT); 1314 static define_unary_function_eval (d2_UTPT_eval,&d2_UTPT,""); 1315 define_unary_function_ptr( D2_UTPT,alias_D2_UTPT,&d2_UTPT_eval); 1316 static unary_function_ptr d_UTPT(int i){ 1317 if (i==1) 1318 return at_zero; 1319 if (i==2) 1320 return D2_UTPT; 1321 return gensizeerr(contextptr); 1322 return 0; 1323 } 1324 static const partial_derivative_multiargs D_UTPT(&d_UTPT); 1325 */ FTS(int ndf,double cs2,double term,int j,double sum)1326 static double FTS(int ndf,double cs2,double term, int j,double sum){ 1327 for(;;){ 1328 term=term*cs2*(ndf+j)/(j+2); 1329 double oldsum=sum; 1330 sum=sum+term; 1331 if (oldsum==sum) 1332 return sum; 1333 j +=2; 1334 } 1335 } FCS(double fac,double term,int j,double res)1336 static double FCS(double fac,double term,int j,double res){ 1337 for (;;){ 1338 j -= 2; 1339 if (j<=1) 1340 return res; 1341 res = ((j-1)*fac*res)/j + term; 1342 } 1343 } 1344 FSS(double ofs,double term,double fac,int j,double sum)1345 static double FSS(double ofs,double term,double fac,int j,double sum){ 1346 for (;;){ 1347 j -= 2; 1348 if (j<=1) 1349 return sum; 1350 sum=(ofs+j-2)/j*fac*sum+term; 1351 } 1352 } FSS2(double ofs,double term,double fac,int j,double sum)1353 static double FSS2(double ofs,double term,double fac,int j,double sum){ 1354 for (;;){ 1355 j -= 2; 1356 if (j<=0) 1357 return sum; 1358 sum=(ofs+j)/j*fac*sum+term; 1359 } 1360 } TTS(int dof,int j,double term,double res,double cs2)1361 static double TTS(int dof,int j,double term,double res,double cs2){ 1362 for (;;){ 1363 j +=2; 1364 term= (term*cs2*(j-1))/j; 1365 double ores= res; 1366 res= res + term/(dof+j); 1367 if (ores==res) return res; 1368 } 1369 } UTPT(const gen & n_orig,const gen & x0,GIAC_CONTEXT)1370 gen UTPT(const gen & n_orig,const gen & x0,GIAC_CONTEXT){ 1371 gen n=n_orig; 1372 if (!is_integral(n)) 1373 return gensizeerr(contextptr); 1374 if (x0==plus_inf) 1375 return 0; 1376 if (x0==minus_inf) 1377 return 1; 1378 gen x1=evalf_double(x0,1,contextptr); 1379 if (n.type!=_INT_ || x1.type!=_DOUBLE_) 1380 return symbolic(at_UTPT,makesequence(n,x0)); 1381 int dof=n.val; 1382 if (dof<=0) 1383 return gendimerr(contextptr); 1384 double x=x1._DOUBLE_val,x2=x*x,y2= x2/dof; 1385 if (0 && dof>=100){ 1386 double y=std::log(y2)+1, a=dof-0.5, b=48*a*a; 1387 y=a*y; 1388 double res = (((((-.4*y - 3.3)*y -24)*y - 85.5)/(.8*y*y + 100 + b)+ y + 3)/b + 1)*std::sqrt(y); 1389 if (x<0) 1390 res=-res; 1391 return _UTPN(res,contextptr); 1392 } 1393 double y=std::sqrt(y2),b= 1+y2,cs2=1/b; 1394 if (x2<25){ 1395 double res; 1396 if (dof==1) 1397 res=0; 1398 else 1399 res=FCS(cs2,y,dof,y); 1400 if (dof %2) 1401 res=2/M_PI*(std::atan(y)+cs2*res); 1402 else 1403 res=res*std::sqrt(cs2); 1404 if (x>0) 1405 return (1-res)/2; 1406 else 1407 return (1+res)/2; 1408 } 1409 else { 1410 double res= TTS(dof,0,dof,1,cs2); 1411 res= FCS(cs2,0,dof+2,res); 1412 if (dof %2) 1413 res= 2/M_PI*std::sqrt(cs2)*res; 1414 res /=2; 1415 if (x<0) 1416 res=1-res; 1417 return res; 1418 } 1419 } _UTPT(const gen & args,GIAC_CONTEXT)1420 gen _UTPT(const gen & args,GIAC_CONTEXT){ 1421 if ( args.type==_STRNG && args.subtype==-1) return args; 1422 if (args.type!=_VECT) 1423 return gensizeerr(contextptr); 1424 vecteur & v=*args._VECTptr; 1425 int s=int(v.size()); 1426 if (s!=2) 1427 return gensizeerr(contextptr); 1428 return UTPT(v[0],v[1],contextptr); 1429 } 1430 static const char _UTPT_s []="UTPT"; 1431 static define_unary_function_eval (__UTPT,&_UTPT,_UTPT_s); 1432 define_unary_function_ptr5( at_UTPT ,alias_at_UTPT,&__UTPT,0,true); 1433 1434 // 26.7.5 in Abramowitz & Stegun utpt_initial_guess(int n,double y,GIAC_CONTEXT)1435 static double utpt_initial_guess(int n,double y,GIAC_CONTEXT){ 1436 // double xp=utpn_initial_guess(y); 1437 double xp=utpn_inverse(y,contextptr)._DOUBLE_val; 1438 double xp2=xp*xp; 1439 double g1xp=xp*(xp2+1)/4; 1440 double g2xp=((5*xp2+16)*xp2+3)*xp/96; 1441 xp=xp+g1xp/n+g2xp/(n*n); 1442 return xp; 1443 } 1444 1445 // dof=degree of freedom student_cdf(const gen & dof0,const gen & x1,const gen & x2,GIAC_CONTEXT)1446 gen student_cdf(const gen & dof0,const gen & x1,const gen & x2,GIAC_CONTEXT){ 1447 gen X2=evalf_double(x2,1,contextptr); 1448 gen X1=evalf_double(x1,1,contextptr); 1449 gen dof(dof0); 1450 if (!is_integral(dof) || dof.val<1 || X1.type!=_DOUBLE_ || X2.type!=_DOUBLE_){ 1451 if (!is_inf(X1) && !is_inf(X2)) 1452 return symbolic(at_student_cdf,gen(makevecteur(dof0,x1,x2),_SEQ__VECT)); 1453 } 1454 return UTPT(dof,X1,contextptr)-UTPT(dof,X2,contextptr); 1455 } _student_cdf(const gen & g,GIAC_CONTEXT)1456 gen _student_cdf(const gen & g,GIAC_CONTEXT){ 1457 if ( g.type==_STRNG && g.subtype==-1) return g; 1458 if (g.type!=_VECT) 1459 return gensizeerr(contextptr); 1460 vecteur & v=*g._VECTptr; 1461 int s=int(v.size()); 1462 if (s==2) 1463 return student_cdf(v[0],minus_inf,v[1],contextptr); 1464 if (s==3) 1465 return student_cdf(v[0],v[1],v[2],contextptr); 1466 return gensizeerr(contextptr); 1467 } 1468 static const char _student_cdf_s []="student_cdf"; 1469 static define_unary_function_eval (__student_cdf,&_student_cdf,_student_cdf_s); 1470 define_unary_function_ptr5( at_student_cdf ,alias_at_student_cdf,&__student_cdf,0,true); 1471 1472 static const char _studentd_cdf_s []="studentd_cdf"; 1473 static define_unary_function_eval (__studentd_cdf,&_student_cdf,_studentd_cdf_s); 1474 define_unary_function_ptr5( at_studentd_cdf ,alias_at_studentd_cdf,&__studentd_cdf,0,true); 1475 student_icdf(const gen & m0,const gen & t_orig,GIAC_CONTEXT)1476 gen student_icdf(const gen & m0,const gen & t_orig,GIAC_CONTEXT){ 1477 gen t=evalf_double(t_orig,1,contextptr); 1478 gen m(m0); 1479 if (!is_integral(m) || m.val<1 || t.type!=_DOUBLE_ || t._DOUBLE_val<0 || t._DOUBLE_val>1) 1480 return symbolic(at_student_icdf,makesequence(m,t)); 1481 if (t._DOUBLE_val==0) 1482 return zero; 1483 if (t._DOUBLE_val==1) 1484 return plus_inf; 1485 double y=t._DOUBLE_val; 1486 double x0=utpt_initial_guess(m.val,1-y,contextptr); 1487 // return x0; 1488 // FIXME: use an iterative method to improve the initial guess 1489 identificateur x(" x"); 1490 gen res=newton(_student_cdf(makesequence(m,x),contextptr)-y,x,x0,NEWTON_DEFAULT_ITERATION,1e-5,1e-12,true,1,0,1,0,.5,contextptr); 1491 if (!is_undef(res)) 1492 return res; 1493 // for example student_icdf(100,0.95) 1494 *logptr(contextptr) << "Low accuracy" << "\n"; 1495 return x0; 1496 } _student_icdf(const gen & g,GIAC_CONTEXT)1497 gen _student_icdf(const gen & g,GIAC_CONTEXT){ 1498 if ( g.type==_STRNG && g.subtype==-1) return g; 1499 if (g.type!=_VECT) 1500 return gensizeerr(contextptr); 1501 vecteur & v=*g._VECTptr; 1502 int s=int(v.size()); 1503 if (s==2) 1504 return student_icdf(v[0],v[1],contextptr); 1505 if (s==3){ 1506 if (v[2]==at_left) 1507 return student_icdf(v[0],v[1],contextptr); 1508 if (v[2]==at_right) 1509 return -student_icdf(v[0],v[1],contextptr); 1510 if (v[2]==at_centre){ 1511 gen t=student_icdf(v[0],(1-v[1])/2,contextptr); 1512 return makevecteur(-t,t); 1513 } 1514 if (v[2]==at_tail){ 1515 gen t=student_icdf(v[0],v[1]/2,contextptr); 1516 return makevecteur(-t,t); 1517 } 1518 return student_icdf(v[0],v[2],contextptr)-student_icdf(v[0],v[1],contextptr); 1519 } 1520 return gensizeerr(contextptr); 1521 } 1522 static const char _student_icdf_s []="student_icdf"; 1523 static define_unary_function_eval (__student_icdf,&_student_icdf,_student_icdf_s); 1524 define_unary_function_ptr5( at_student_icdf ,alias_at_student_icdf,&__student_icdf,0,true); 1525 1526 static const char _studentd_icdf_s []="studentd_icdf"; 1527 static define_unary_function_eval (__studentd_icdf,&_student_icdf,_studentd_icdf_s); 1528 define_unary_function_ptr5( at_studentd_icdf ,alias_at_studentd_icdf,&__studentd_icdf,0,true); 1529 chisquare(const gen & n,const gen & x,GIAC_CONTEXT)1530 gen chisquare(const gen & n,const gen & x,GIAC_CONTEXT){ 1531 if (x.type==_VECT) 1532 return apply2nd(n,x,contextptr,chisquare); 1533 gen n2=n/2; 1534 return rdiv(pow(x,n2-1,contextptr)*exp(-x/2,contextptr),Gamma(n2,contextptr)*pow(2,n2,contextptr),contextptr); 1535 } _chisquare(const gen & g,GIAC_CONTEXT)1536 gen _chisquare(const gen & g,GIAC_CONTEXT){ 1537 if ( g.type==_STRNG && g.subtype==-1) return g; 1538 if (g.type!=_VECT){ 1539 if (abs_calc_mode(contextptr)==38) 1540 return symbolic(at_chisquare,g); 1541 return symbolic(at_chisquared,g); 1542 } 1543 vecteur & v=*g._VECTptr; 1544 int s=int(v.size()); 1545 if (s==2) 1546 return chisquare(v[0],v[1],contextptr); 1547 return gensizeerr(contextptr); 1548 } 1549 static const char _chisquare_s []="chisquare"; 1550 static define_unary_function_eval (__chisquare,&_chisquare,_chisquare_s); 1551 define_unary_function_ptr5( at_chisquare ,alias_at_chisquare,&__chisquare,0,true); 1552 1553 static const char _chisquared_s []="chisquared"; 1554 static define_unary_function_eval (__chisquared,&_chisquare,_chisquared_s); 1555 define_unary_function_ptr5( at_chisquared ,alias_at_chisquared,&__chisquared,0,true); 1556 1557 /* statically handled by derive.cc 1558 static gen d2_UTPC(const gen & e,GIAC_CONTEXT){ 1559 return -_chisquare(e,contextptr); 1560 } 1561 static const partial_derivative_onearg D_at_UTPC(d2_UTPC); 1562 static define_unary_function_eval (d2_UTPC_eval,&d2_UTPC,""); 1563 define_unary_function_ptr( D2_UTPC,alias_D2_UTPC,&d2_UTPC_eval); 1564 static unary_function_ptr d_UTPC(int i){ 1565 if (i==1) 1566 return at_zero; 1567 if (i==2) 1568 return D2_UTPC; 1569 return gensizeerr(contextptr); 1570 } 1571 static const partial_derivative_multiargs D_UTPC(&d_UTPC); 1572 */ UTPC(const gen & n_orig,const gen & x0,GIAC_CONTEXT)1573 gen UTPC(const gen & n_orig,const gen & x0,GIAC_CONTEXT){ 1574 gen dof=n_orig; 1575 if (x0==plus_inf) 1576 return 0; 1577 if (is_zero(x0)) 1578 return 1; 1579 gen x1=evalf_double(x0,1,contextptr); 1580 if (!is_integral(dof) || x1.type!=_DOUBLE_) 1581 return symbolic(at_UTPC,gen(makevecteur(dof,x0),_SEQ__VECT)); // gensizeerr(contextptr); 1582 int n=dof.val; 1583 double x=x1._DOUBLE_val; 1584 if (x<0) 1585 return 1; 1586 if (x>10000) 1587 return 0.0; 1588 if (n<1) 1589 return gensizeerr(contextptr); 1590 if (n==1) 1591 return 2*_UTPN(sqrt(x,contextptr),contextptr); 1592 if (n>100){ 1593 } 1594 double res=1; 1595 if (x>2){ 1596 int r=n%2+2; 1597 res=std::exp(-x/2); 1598 double term = res; 1599 for (;r<n;r += 2){ 1600 term = term*x/r; 1601 res += term; 1602 } 1603 } 1604 else { 1605 int r=n-2; 1606 for (;r>1;r-=2){ 1607 res = res*x/r+1; 1608 } 1609 res *= std::exp(-x/2); 1610 } 1611 if (n%2) 1612 return std::sqrt(2*x/M_PI)*res+2*_UTPN(sqrt(x,contextptr),contextptr); 1613 else 1614 return res; 1615 } _UTPC(const gen & args,GIAC_CONTEXT)1616 gen _UTPC(const gen & args,GIAC_CONTEXT){ 1617 if ( args.type==_STRNG && args.subtype==-1) return args; 1618 if (args.type!=_VECT) 1619 return gensizeerr(contextptr); 1620 vecteur & v=*args._VECTptr; 1621 int s=int(v.size()); 1622 if (s!=2) 1623 return gensizeerr(contextptr); 1624 return UTPC(v[0],v[1],contextptr); 1625 } 1626 static const char _UTPC_s []="UTPC"; 1627 static define_unary_function_eval (__UTPC,&_UTPC,_UTPC_s); 1628 define_unary_function_ptr5( at_UTPC ,alias_at_UTPC,&__UTPC,0,true); 1629 chisquare_cdf(const gen & dof,const gen & x1,const gen & x2,GIAC_CONTEXT)1630 gen chisquare_cdf(const gen & dof,const gen & x1,const gen & x2,GIAC_CONTEXT){ 1631 return UTPC(dof,x1,contextptr)-UTPC(dof,x2,contextptr); 1632 } _chisquare_cdf(const gen & g,GIAC_CONTEXT)1633 gen _chisquare_cdf(const gen & g,GIAC_CONTEXT){ 1634 if ( g.type==_STRNG && g.subtype==-1) return g; 1635 if (g.type!=_VECT) 1636 return gensizeerr(contextptr); 1637 vecteur & v=*g._VECTptr; 1638 int s=int(v.size()); 1639 if (s==2) 1640 return chisquare_cdf(v[0],0,v[1],contextptr); 1641 if (s==3) 1642 return chisquare_cdf(v[0],v[1],v[2],contextptr); 1643 return gensizeerr(contextptr); 1644 } 1645 static const char _chisquare_cdf_s []="chisquare_cdf"; 1646 static define_unary_function_eval (__chisquare_cdf,&_chisquare_cdf,_chisquare_cdf_s); 1647 define_unary_function_ptr5( at_chisquare_cdf ,alias_at_chisquare_cdf,&__chisquare_cdf,0,true); 1648 1649 static const char _chisquared_cdf_s []="chisquared_cdf"; 1650 static define_unary_function_eval (__chisquared_cdf,&_chisquare_cdf,_chisquared_cdf_s); 1651 define_unary_function_ptr5( at_chisquared_cdf ,alias_at_chisquared_cdf,&__chisquared_cdf,0,true); 1652 1653 // Abramowitz & Stegun 26.4.17 utpc_initial_guess(int n,double y,GIAC_CONTEXT)1654 static double utpc_initial_guess(int n,double y,GIAC_CONTEXT){ 1655 // double xp=utpn_initial_guess(y); 1656 if (n==2) 1657 return -2*std::log(y); 1658 if (n==1) 1659 y=y/2; 1660 double xp=utpn_inverse(y,contextptr)._DOUBLE_val; 1661 if (n==1) 1662 return xp*xp; 1663 double d=2/(9.0*n); 1664 d=1+xp*std::sqrt(d)-d; 1665 return n*d*d*d; 1666 } 1667 chisquare_icdf(const gen & m0,const gen & t_orig,GIAC_CONTEXT)1668 gen chisquare_icdf(const gen & m0,const gen & t_orig,GIAC_CONTEXT){ 1669 gen t=evalf_double(t_orig,1,contextptr); 1670 gen m(m0); 1671 if (!is_integral(m) || t.type!=_DOUBLE_ || t._DOUBLE_val<0 || t._DOUBLE_val>1) 1672 return gensizeerr(contextptr); 1673 if (t._DOUBLE_val==0) 1674 return zero; 1675 if (t._DOUBLE_val==1) 1676 return plus_inf; 1677 // return utpc_initial_guess(m.val,1-t._DOUBLE_val); 1678 double x0=utpc_initial_guess(m.val,1-t._DOUBLE_val,contextptr); 1679 // FIXME 1680 identificateur x(" z"); 1681 return newton(1-UTPC(m,x,contextptr)-t,x,x0,NEWTON_DEFAULT_ITERATION,1e-5,1e-12,true,1,0,1,0,.5,contextptr); 1682 } _chisquare_icdf(const gen & g,GIAC_CONTEXT)1683 gen _chisquare_icdf(const gen & g,GIAC_CONTEXT){ 1684 if ( g.type==_STRNG && g.subtype==-1) return g; 1685 if (g.type!=_VECT) 1686 return gensizeerr(contextptr); 1687 vecteur & v=*g._VECTptr; 1688 int s=int(v.size()); 1689 if (s==2) 1690 return chisquare_icdf(v[0],v[1],contextptr); 1691 if (s==3){ 1692 if (v[2]==at_left) 1693 return chisquare_icdf(v[0],v[1],contextptr); 1694 if (v[2]==at_right) 1695 return chisquare_icdf(v[0],1-v[1],contextptr); 1696 if (v[2]==at_centre) 1697 return makevecteur(chisquare_icdf(v[0],(1-v[1])/2,contextptr),chisquare_icdf(v[0],(1+v[1])/2,contextptr)); 1698 if (v[2]==at_tail) 1699 return makevecteur(chisquare_icdf(v[0],(v[1])/2,contextptr),chisquare_icdf(v[0],1-v[1]/2,contextptr)); 1700 return chisquare_icdf(v[0],v[2],contextptr)-chisquare_icdf(v[0],v[1],contextptr); 1701 } 1702 return gensizeerr(contextptr); 1703 } 1704 static const char _chisquare_icdf_s []="chisquare_icdf"; 1705 static define_unary_function_eval (__chisquare_icdf,&_chisquare_icdf,_chisquare_icdf_s); 1706 define_unary_function_ptr5( at_chisquare_icdf ,alias_at_chisquare_icdf,&__chisquare_icdf,0,true); 1707 1708 static const char _chisquared_icdf_s []="chisquared_icdf"; 1709 static define_unary_function_eval (__chisquared_icdf,&_chisquare_icdf,_chisquared_icdf_s); 1710 define_unary_function_ptr5( at_chisquared_icdf ,alias_at_chisquared_icdf,&__chisquared_icdf,0,true); 1711 snedecor(const gen & a,const gen & b,const gen & x,GIAC_CONTEXT)1712 gen snedecor(const gen & a,const gen & b,const gen & x,GIAC_CONTEXT){ 1713 if (x.type==_VECT) 1714 return apply3rd(a,b,x,contextptr,snedecor); 1715 if (is_positive(-x,contextptr)) 1716 return zero; 1717 return pow(a/b,a/2,contextptr)/Beta(a/2,b/2,contextptr) * pow(x,a/2-1,contextptr) * pow(1+a/b*x,-(a+b)/2,contextptr); 1718 } _snedecor(const gen & g,GIAC_CONTEXT)1719 gen _snedecor(const gen & g,GIAC_CONTEXT){ 1720 if ( g.type==_STRNG && g.subtype==-1) return g; 1721 if (g.type!=_VECT) 1722 return gensizeerr(contextptr); 1723 vecteur & v=*g._VECTptr; 1724 int s=int(v.size()); 1725 if (s==2){ 1726 if (abs_calc_mode(contextptr)==38) 1727 return symbolic(at_fisher,g); 1728 return symbolic(at_fisherd,g); 1729 } 1730 if (s==3) 1731 return snedecor(v[0],v[1],v[2],contextptr); 1732 return gensizeerr(contextptr); 1733 } 1734 static const char _snedecor_s []="snedecor"; 1735 static define_unary_function_eval (__snedecor,&_snedecor,_snedecor_s); 1736 define_unary_function_ptr5( at_snedecor ,alias_at_snedecor,&__snedecor,0,true); 1737 1738 static const char _snedecord_s []="snedecord"; 1739 static define_unary_function_eval (__snedecord,&_snedecor,_snedecord_s); 1740 define_unary_function_ptr5( at_snedecord ,alias_at_snedecord,&__snedecord,0,true); 1741 1742 static const char _fisher_s []="fisher"; 1743 static define_unary_function_eval (__fisher,&_snedecor,_fisher_s); 1744 define_unary_function_ptr5( at_fisher ,alias_at_fisher,&__fisher,0,true); 1745 1746 static const char _fisherd_s []="fisherd"; 1747 static define_unary_function_eval (__fisherd,&_snedecor,_fisherd_s); 1748 define_unary_function_ptr5( at_fisherd ,alias_at_fisherd,&__fisherd,0,true); 1749 1750 /* statically handled in derive.cc 1751 static gen d2_UTPF(const gen & e,GIAC_CONTEXT){ 1752 return -_snedecor(e,contextptr); 1753 } 1754 static const partial_derivative_onearg D_at_UTPF(d2_UTPF); 1755 static define_unary_function_eval (d2_UTPF_eval,&d2_UTPF,""); 1756 define_unary_function_ptr( D2_UTPF,alias_D2_UTPF,&d2_UTPF_eval); 1757 static unary_function_ptr d_UTPF(int i){ 1758 if (i<3) 1759 return at_zero; 1760 if (i==3) 1761 return D2_UTPF; 1762 return gensizeerr(contextptr); 1763 } 1764 static const partial_derivative_multiargs D_UTPF(&d_UTPF); 1765 */ UTPF(const gen & num,const gen & den,const gen & x0,GIAC_CONTEXT)1766 gen UTPF(const gen & num,const gen & den,const gen & x0,GIAC_CONTEXT){ 1767 gen gndf=num,gddf=den,gx=evalf_double(x0,1,contextptr); 1768 if (!is_integral(gndf) || !is_integral(gddf) || gx.type!=_DOUBLE_) 1769 return symbolic(at_UTPF,gen(makevecteur(num,den,x0),_SEQ__VECT)); // gensizeerr(contextptr); 1770 if (gx._DOUBLE_val<=0) 1771 return plus_one; 1772 int ndf=gndf.val,ddf=gddf.val; 1773 double x=gx._DOUBLE_val; 1774 if (ndf<1 || ddf <1 || ndf>300 || ddf>300) 1775 return gendimerr(contextptr); 1776 if (ndf==1) 1777 return 2*UTPT(ddf,std::sqrt(x),contextptr); 1778 double y= (x*ndf)/ddf, sn2= y/(1+y), cs2= 1/(1+y),sum; 1779 y=std::sqrt(y); 1780 if (ndf%2){ 1781 if (ddf%2){ // ndf && ddf odd 1782 if (y<1){ 1783 if (ddf==1) sum=0; else sum=1; 1784 sum=y*cs2*FCS(cs2,1,ddf,sum)+std::atan(y); 1785 sum=1-2/M_PI*sum; 1786 double sumB=ddf*FSS(ddf,1,sn2,ndf,1); 1787 sumB=FCS(cs2,0,ddf+2,sumB); 1788 return sum+2/M_PI*cs2*y*sumB; 1789 } 1790 else { // y>=1 1791 sum= FTS(ndf,cs2,1,ddf,1); 1792 sum= FSS2(ddf,0,sn2,ndf,sum); 1793 sum= FCS(cs2,0,ddf+2,sum); 1794 return 2/M_PI*y*cs2*sum; 1795 } 1796 } // end ndf odd , ddf odd 1797 else { // ndf odd, ddf even 1798 if (y<1) 1799 return 1-FSS(ndf,1,cs2,ddf,1)*std::pow(sn2,ndf/2.0); 1800 else { 1801 sum= FTS(ndf,cs2,1,ddf,1); 1802 sum= FSS(ndf,0,cs2,ddf+2,sum); 1803 return sum*std::pow(sn2,ndf/2.0); 1804 } 1805 } 1806 } // end ndf odd 1807 else { // ndf even 1808 return FSS(ddf,1,sn2,ndf,1)*std::pow(cs2,ddf/2.0); 1809 } 1810 } _UTPF(const gen & args,GIAC_CONTEXT)1811 gen _UTPF(const gen & args,GIAC_CONTEXT){ 1812 if ( args.type==_STRNG && args.subtype==-1) return args; 1813 if (args.type!=_VECT) 1814 return gensizeerr(contextptr); 1815 vecteur & v=*args._VECTptr; 1816 int s=int(v.size()); 1817 if (s!=3) 1818 return gensizeerr(contextptr); 1819 return UTPF(v[0],v[1],v[2],contextptr); 1820 } 1821 static const char _UTPF_s []="UTPF"; 1822 static define_unary_function_eval (__UTPF,&_UTPF,_UTPF_s); 1823 define_unary_function_ptr5( at_UTPF ,alias_at_UTPF,&__UTPF,0,true); 1824 snedecor_cdf(const gen & ndof,const gen & ddof,const gen & x,GIAC_CONTEXT)1825 gen snedecor_cdf(const gen & ndof,const gen & ddof,const gen & x,GIAC_CONTEXT){ 1826 gen gndf(ndof),gddf(ddof),gx(x); 1827 if (!is_integral(gndf) || !is_integral(gddf)) 1828 return gentypeerr(contextptr); 1829 if (gx.type!=_DOUBLE_){ 1830 if (1) {// (calc_mode(contextptr)==1) 1831 if (is_inf(x)) 1832 return symbolic(at_Beta,makesequence(ndof/2,ddof/2,1,1)); 1833 return symbolic(at_Beta,makesequence(ndof/2,ddof/2,ndof*x/(ndof*x+ddof),1)); 1834 } 1835 else 1836 return symbolic(at_snedecor_cdf,makesequence(ndof,ddof,x)); 1837 } 1838 return 1-UTPF(ndof,ddof,x,contextptr); 1839 } _snedecor_cdf(const gen & g,GIAC_CONTEXT)1840 gen _snedecor_cdf(const gen & g,GIAC_CONTEXT){ 1841 if ( g.type==_STRNG && g.subtype==-1) return g; 1842 if (g.type!=_VECT) 1843 return gensizeerr(contextptr); 1844 vecteur & v=*g._VECTptr; 1845 int s=int(v.size()); 1846 if (s==3) 1847 return snedecor_cdf(v[0],v[1],v[2],contextptr); 1848 if (s==4) 1849 return snedecor_cdf(v[0],v[1],v[3],contextptr)-snedecor_cdf(v[0],v[1],v[2],contextptr); 1850 return gensizeerr(contextptr); 1851 } 1852 static const char _snedecor_cdf_s []="snedecor_cdf"; 1853 static define_unary_function_eval (__snedecor_cdf,&_snedecor_cdf,_snedecor_cdf_s); 1854 define_unary_function_ptr5( at_snedecor_cdf ,alias_at_snedecor_cdf,&__snedecor_cdf,0,true); 1855 static const char _snedecord_cdf_s []="snedecord_cdf"; 1856 static define_unary_function_eval (__snedecord_cdf,&_snedecor_cdf,_snedecord_cdf_s); 1857 define_unary_function_ptr5( at_snedecord_cdf ,alias_at_snedecord_cdf,&__snedecord_cdf,0,true); 1858 static const char _fisher_cdf_s []="fisher_cdf"; 1859 static define_unary_function_eval (__fisher_cdf,&_snedecor_cdf,_fisher_cdf_s); 1860 define_unary_function_ptr5( at_fisher_cdf ,alias_at_fisher_cdf,&__fisher_cdf,0,true); 1861 static const char _fisherd_cdf_s []="fisherd_cdf"; 1862 static define_unary_function_eval (__fisherd_cdf,&_snedecor_cdf,_fisherd_cdf_s); 1863 define_unary_function_ptr5( at_fisherd_cdf ,alias_at_fisherd_cdf,&__fisherd_cdf,0,true); 1864 1865 // Abramowitz & Stegun 26.6.16 utpf_initial_guess(int num,int den,double y,GIAC_CONTEXT)1866 static double utpf_initial_guess(int num,int den,double y,GIAC_CONTEXT){ 1867 if (num==1){ 1868 double xp=utpt_initial_guess(den,y/2,contextptr); 1869 return xp*xp; 1870 } 1871 if (den==1){ 1872 return y-0.5; 1873 } 1874 double xp=utpn_inverse(y,contextptr)._DOUBLE_val; 1875 double lambda=(xp*xp-3)/6; 1876 double h=2/fabs(1.0/(num-1)+1.0/(den-1)); // harmonic 1877 double w=xp*std::sqrt(h+lambda)/h-(lambda+5.0/6.0-2/(3*h))*fabs(1.0/(num-1)-1.0/(den-1)); 1878 return std::exp(2*w); 1879 } 1880 snedecor_icdf(const gen & num0,const gen & den0,const gen & t_orig,GIAC_CONTEXT)1881 gen snedecor_icdf(const gen & num0,const gen & den0,const gen & t_orig,GIAC_CONTEXT){ 1882 gen t=evalf_double(t_orig,1,contextptr); 1883 gen num(num0),den(den0); 1884 if (!is_integral(num) || !is_integral(den) || num.val<0 || den.val<0 || t.type!=_DOUBLE_ || t._DOUBLE_val<0 || t._DOUBLE_val>1) 1885 return gensizeerr(contextptr); 1886 if (t._DOUBLE_val==0) 1887 return zero; 1888 if (t._DOUBLE_val==1) 1889 return plus_inf; 1890 // return utpf_initial_guess(num.val,den.val,1-t._DOUBLE_val); 1891 double x0=utpf_initial_guess(num.val,den.val,1-t._DOUBLE_val,contextptr); 1892 // FIXME 1893 identificateur x(" z"); 1894 return newton(1-UTPF(num,den,x,contextptr)-t,x,x0,NEWTON_DEFAULT_ITERATION,1e-5,1e-12,true,0,1.79769313486e+308,1,0,.5,contextptr); 1895 } _snedecor_icdf(const gen & g,GIAC_CONTEXT)1896 gen _snedecor_icdf(const gen & g,GIAC_CONTEXT){ 1897 if ( g.type==_STRNG && g.subtype==-1) return g; 1898 if (g.type!=_VECT) 1899 return gensizeerr(contextptr); 1900 vecteur & v=*g._VECTptr; 1901 int s=int(v.size()); 1902 if (s==3) 1903 return snedecor_icdf(v[0],v[1],v[2],contextptr); 1904 if (s==4){ 1905 if (v[3]==at_left) 1906 return snedecor_icdf(v[0],v[1],v[2],contextptr); 1907 if (v[3]==at_right) 1908 return snedecor_icdf(v[0],v[1],1-v[2],contextptr); 1909 if (v[3]==at_centre) 1910 return makevecteur(snedecor_icdf(v[0],v[1],(1-v[2])/2,contextptr),snedecor_icdf(v[0],v[1],(1+v[2])/2,contextptr)); 1911 if (v[3]==at_tail) 1912 return makevecteur(snedecor_icdf(v[0],v[1],(v[2])/2,contextptr),snedecor_icdf(v[0],v[1],1-v[2]/2,contextptr)); 1913 return snedecor_icdf(v[0],v[1],v[3],contextptr)-snedecor_icdf(v[0],v[1],v[2],contextptr); 1914 } 1915 return gensizeerr(contextptr); 1916 } 1917 static const char _snedecor_icdf_s []="snedecor_icdf"; 1918 static define_unary_function_eval (__snedecor_icdf,&_snedecor_icdf,_snedecor_icdf_s); 1919 define_unary_function_ptr5( at_snedecor_icdf ,alias_at_snedecor_icdf,&__snedecor_icdf,0,true); 1920 static const char _snedecord_icdf_s []="snedecord_icdf"; 1921 static define_unary_function_eval (__snedecord_icdf,&_snedecor_icdf,_snedecord_icdf_s); 1922 define_unary_function_ptr5( at_snedecord_icdf ,alias_at_snedecord_icdf,&__snedecord_icdf,0,true); 1923 1924 static const char _fisher_icdf_s []="fisher_icdf"; 1925 static define_unary_function_eval (__fisher_icdf,&_snedecor_icdf,_fisher_icdf_s); 1926 define_unary_function_ptr5( at_fisher_icdf ,alias_at_fisher_icdf,&__fisher_icdf,0,true); 1927 static const char _fisherd_icdf_s []="fisherd_icdf"; 1928 static define_unary_function_eval (__fisherd_icdf,&_snedecor_icdf,_fisherd_icdf_s); 1929 define_unary_function_ptr5( at_fisherd_icdf ,alias_at_fisherd_icdf,&__fisherd_icdf,0,true); 1930 cauchy(const gen & x0,const gen & a,const gen & x,GIAC_CONTEXT)1931 gen cauchy(const gen & x0,const gen & a,const gen & x,GIAC_CONTEXT){ 1932 if (x.type==_VECT) 1933 return apply3rd(x0,a,x,contextptr,cauchy); 1934 return plus_one/cst_pi*a/(pow(x-x0,2,contextptr)+pow(a,2,contextptr)); 1935 } _cauchy(const gen & g,GIAC_CONTEXT)1936 gen _cauchy(const gen & g,GIAC_CONTEXT){ 1937 if ( g.type==_STRNG && g.subtype==-1) return g; 1938 if (g.type!=_VECT) 1939 return cauchy(0,1,g,contextptr); 1940 vecteur & v=*g._VECTptr; 1941 int s=int(v.size()); 1942 if (s==2) 1943 return symbolic(at_cauchyd,g); 1944 if (s==3) 1945 return cauchy(v[0],v[1],v[2],contextptr); 1946 return gensizeerr(contextptr); 1947 } 1948 static const char _cauchy_s []="cauchy"; 1949 static define_unary_function_eval (__cauchy,&_cauchy,_cauchy_s); 1950 define_unary_function_ptr5( at_cauchy ,alias_at_cauchy,&__cauchy,0,true); 1951 1952 static const char _cauchyd_s []="cauchyd"; 1953 static define_unary_function_eval (__cauchyd,&_cauchy,_cauchyd_s); 1954 define_unary_function_ptr5( at_cauchyd ,alias_at_cauchyd,&__cauchyd,0,true); 1955 cauchy_cdf(const gen & x0,const gen & a,const gen & x,GIAC_CONTEXT)1956 gen cauchy_cdf(const gen &x0,const gen &a,const gen & x,GIAC_CONTEXT){ 1957 return plus_one_half+atan((x-x0)/a,contextptr)/cst_pi; 1958 } _cauchy_cdf(const gen & g,GIAC_CONTEXT)1959 gen _cauchy_cdf(const gen & g,GIAC_CONTEXT){ 1960 if ( g.type==_STRNG && g.subtype==-1) return g; 1961 if (g.type!=_VECT) 1962 return cauchy_cdf(0,1,g,contextptr); 1963 vecteur & v=*g._VECTptr; 1964 int s=int(v.size()); 1965 if (s==3) 1966 return cauchy_cdf(v[0],v[1],v[2],contextptr); 1967 if (s==4) 1968 return cauchy_cdf(v[0],v[1],v[3],contextptr)-cauchy_cdf(v[0],v[1],v[2],contextptr); 1969 return gensizeerr(contextptr); 1970 } 1971 static const char _cauchy_cdf_s []="cauchy_cdf"; 1972 static define_unary_function_eval (__cauchy_cdf,&_cauchy_cdf,_cauchy_cdf_s); 1973 define_unary_function_ptr5( at_cauchy_cdf ,alias_at_cauchy_cdf,&__cauchy_cdf,0,true); 1974 1975 static const char _cauchyd_cdf_s []="cauchyd_cdf"; 1976 static define_unary_function_eval (__cauchyd_cdf,&_cauchy_cdf,_cauchyd_cdf_s); 1977 define_unary_function_ptr5( at_cauchyd_cdf ,alias_at_cauchyd_cdf,&__cauchyd_cdf,0,true); 1978 cauchy_icdf(const gen & x0,const gen & a,const gen & x,GIAC_CONTEXT)1979 gen cauchy_icdf(const gen &x0,const gen &a,const gen & x,GIAC_CONTEXT){ 1980 return tan(cst_pi*(x-plus_one_half),contextptr)*a+x0; 1981 } _cauchy_icdf(const gen & g,GIAC_CONTEXT)1982 gen _cauchy_icdf(const gen & g,GIAC_CONTEXT){ 1983 if ( g.type==_STRNG && g.subtype==-1) return g; 1984 if (g.type!=_VECT) 1985 return cauchy_icdf(0,1,g,contextptr); 1986 vecteur & v=*g._VECTptr; 1987 int s=int(v.size()); 1988 if (s==3) 1989 return cauchy_icdf(v[0],v[1],v[2],contextptr); 1990 if (s==4) 1991 return cauchy_icdf(v[0],v[1],v[3],contextptr)-cauchy_icdf(v[0],v[1],v[2],contextptr); 1992 return gensizeerr(contextptr); 1993 } 1994 static const char _cauchy_icdf_s []="cauchy_icdf"; 1995 static define_unary_function_eval (__cauchy_icdf,&_cauchy_icdf,_cauchy_icdf_s); 1996 define_unary_function_ptr5( at_cauchy_icdf ,alias_at_cauchy_icdf,&__cauchy_icdf,0,true); 1997 1998 static const char _cauchyd_icdf_s []="cauchyd_icdf"; 1999 static define_unary_function_eval (__cauchyd_icdf,&_cauchy_icdf,_cauchyd_icdf_s); 2000 define_unary_function_ptr5( at_cauchyd_icdf ,alias_at_cauchyd_icdf,&__cauchyd_icdf,0,true); 2001 weibull(const gen & k,const gen & lambda,const gen & theta,const gen & x,GIAC_CONTEXT)2002 gen weibull(const gen & k,const gen & lambda,const gen & theta,const gen & x,GIAC_CONTEXT){ 2003 gen tmp=(x-theta)/lambda; 2004 return k/lambda*pow(tmp,k-1,contextptr)*exp(-pow(tmp,k,contextptr),contextptr); 2005 } _weibull(const gen & g,GIAC_CONTEXT)2006 gen _weibull(const gen & g,GIAC_CONTEXT){ 2007 if ( g.type==_STRNG && g.subtype==-1) return g; 2008 if (g.type!=_VECT) 2009 return gensizeerr(contextptr); 2010 vecteur & v=*g._VECTptr; 2011 int s=int(v.size()); 2012 if (s==2) 2013 return symbolic(at_weibulld,g); 2014 if (s==3) 2015 return weibull(v[0],v[1],0,v[2],contextptr); 2016 if (s==4) 2017 return weibull(v[0],v[1],v[2],v[3],contextptr); 2018 return gensizeerr(contextptr); 2019 } 2020 static const char _weibull_s []="weibull"; 2021 static define_unary_function_eval (__weibull,&_weibull,_weibull_s); 2022 define_unary_function_ptr5( at_weibull ,alias_at_weibull,&__weibull,0,true); 2023 2024 static const char _weibulld_s []="weibulld"; 2025 static define_unary_function_eval (__weibulld,&_weibull,_weibulld_s); 2026 define_unary_function_ptr5( at_weibulld ,alias_at_weibulld,&__weibulld,0,true); 2027 weibull_cdf(const gen & k,const gen & lambda,const gen & theta,const gen & x,GIAC_CONTEXT)2028 gen weibull_cdf(const gen & k,const gen & lambda,const gen & theta,const gen & x,GIAC_CONTEXT){ 2029 gen tmp=(x-theta)/lambda; 2030 return 1-exp(-pow(tmp,k,contextptr),contextptr); 2031 } _weibull_cdf(const gen & g,GIAC_CONTEXT)2032 gen _weibull_cdf(const gen & g,GIAC_CONTEXT){ 2033 if ( g.type==_STRNG && g.subtype==-1) return g; 2034 if (g.type!=_VECT) 2035 return gensizeerr(contextptr); 2036 vecteur & v=*g._VECTptr; 2037 int s=int(v.size()); 2038 if (s==3) 2039 return weibull_cdf(v[0],v[1],0,v[2],contextptr); 2040 if (s==4) 2041 return weibull_cdf(v[0],v[1],v[2],v[3],contextptr); 2042 if (s==5) 2043 return weibull_cdf(v[0],v[1],v[2],v[4],contextptr)-weibull_cdf(v[0],v[1],v[2],v[3],contextptr); 2044 return gensizeerr(contextptr); 2045 } 2046 static const char _weibull_cdf_s []="weibull_cdf"; 2047 static define_unary_function_eval (__weibull_cdf,&_weibull_cdf,_weibull_cdf_s); 2048 define_unary_function_ptr5( at_weibull_cdf ,alias_at_weibull_cdf,&__weibull_cdf,0,true); 2049 2050 static const char _weibulld_cdf_s []="weibulld_cdf"; 2051 static define_unary_function_eval (__weibulld_cdf,&_weibull_cdf,_weibulld_cdf_s); 2052 define_unary_function_ptr5( at_weibulld_cdf ,alias_at_weibulld_cdf,&__weibulld_cdf,0,true); 2053 weibull_icdf(const gen & k,const gen & lambda,const gen & theta,const gen & y,GIAC_CONTEXT)2054 gen weibull_icdf(const gen & k,const gen & lambda,const gen & theta,const gen & y,GIAC_CONTEXT){ 2055 // solve(1-exp(-((x-theta)/lambda)^k)=y,x) 2056 // x=lambda*(-ln(-(y-1)))^(1/k)+theta 2057 return lambda*pow(-ln(1-y,contextptr),plus_one/k,contextptr)+theta; 2058 } _weibull_icdf(const gen & g,GIAC_CONTEXT)2059 gen _weibull_icdf(const gen & g,GIAC_CONTEXT){ 2060 if ( g.type==_STRNG && g.subtype==-1) return g; 2061 if (g.type!=_VECT) 2062 return gensizeerr(contextptr); 2063 vecteur & v=*g._VECTptr; 2064 int s=int(v.size()); 2065 if (s==3) 2066 return weibull_icdf(v[0],v[1],0,v[2],contextptr); 2067 if (s==4) 2068 return weibull_icdf(v[0],v[1],v[2],v[3],contextptr); 2069 if (s==5) 2070 return weibull_icdf(v[0],v[1],v[2],v[4],contextptr)-weibull_icdf(v[0],v[1],v[2],v[3],contextptr); 2071 return gensizeerr(contextptr); 2072 } 2073 static const char _weibull_icdf_s []="weibull_icdf"; 2074 static define_unary_function_eval (__weibull_icdf,&_weibull_icdf,_weibull_icdf_s); 2075 define_unary_function_ptr5( at_weibull_icdf ,alias_at_weibull_icdf,&__weibull_icdf,0,true); 2076 2077 static const char _weibulld_icdf_s []="weibulld_icdf"; 2078 static define_unary_function_eval (__weibulld_icdf,&_weibull_icdf,_weibulld_icdf_s); 2079 define_unary_function_ptr5( at_weibulld_icdf ,alias_at_weibulld_icdf,&__weibulld_icdf,0,true); 2080 _randweibulld(const gen & args,GIAC_CONTEXT)2081 gen _randweibulld(const gen & args,GIAC_CONTEXT){ 2082 if (args.type==_STRNG && args.subtype==-1) return args; 2083 if (args.type!=_VECT || args._VECTptr->size()!=2) 2084 return gensizeerr(contextptr); 2085 gen k=args._VECTptr->front(); 2086 gen lambda=args._VECTptr->back(); 2087 k=evalf_double(k,1,contextptr); 2088 lambda=evalf_double(lambda,1,contextptr); 2089 if (is_positive(-k,contextptr) || is_positive(-lambda,contextptr) || k.type!=_DOUBLE_ || lambda.type!=_DOUBLE_) 2090 return gensizeerr(contextptr); 2091 return std::pow(exp_rand(contextptr),1.0/k._DOUBLE_val)*lambda; 2092 } 2093 static const char _randweibulld_s []="randweibulld"; 2094 static define_unary_function_eval (__randweibulld,&_randweibulld,_randweibulld_s); 2095 define_unary_function_ptr5( at_randweibulld ,alias_at_randweibulld,&__randweibulld,0,true); 2096 2097 static const char _weibullvariate_s []="weibullvariate"; 2098 static define_unary_function_eval (__weibullvariate,&_randweibulld,_weibullvariate_s); 2099 define_unary_function_ptr5( at_weibullvariate ,alias_at_weibullvariate,&__weibullvariate,0,true); 2100 betad(const gen & alpha,const gen & beta,const gen & x,GIAC_CONTEXT)2101 gen betad(const gen &alpha,const gen & beta,const gen & x,GIAC_CONTEXT){ 2102 if ( (x==0 && alpha==1) || (x==1 && beta==1)) 2103 return plus_one/Beta(alpha,beta,contextptr); 2104 return pow(x,alpha-1,contextptr)*pow(1-x,beta-1,contextptr)/Beta(alpha,beta,contextptr); 2105 } _betad(const gen & g,GIAC_CONTEXT)2106 gen _betad(const gen & g,GIAC_CONTEXT){ 2107 if ( g.type==_STRNG && g.subtype==-1) return g; 2108 if (g.type!=_VECT) 2109 return betad(0,1,g,contextptr); 2110 vecteur & v=*g._VECTptr; 2111 int s=int(v.size()); 2112 if (s==2) 2113 return symbolic(at_betad,g); 2114 if (s==3) 2115 return betad(v[0],v[1],v[2],contextptr); 2116 return gensizeerr(contextptr); 2117 } 2118 static const char _betad_s []="betad"; 2119 static define_unary_function_eval (__betad,&_betad,_betad_s); 2120 define_unary_function_ptr5( at_betad ,alias_at_betad,&__betad,0,true); 2121 2122 // beta_cdf=Beta regularized _betad_cdf(const gen & g,GIAC_CONTEXT)2123 gen _betad_cdf(const gen & g,GIAC_CONTEXT){ 2124 if ( g.type==_STRNG && g.subtype==-1) return g; 2125 if (g.type!=_VECT) 2126 return gensizeerr(contextptr); 2127 vecteur & v=*g._VECTptr; 2128 int s=int(v.size()); 2129 if (s==3) 2130 return _Beta(makesequence(v[0],v[1],v[2],1),contextptr); 2131 if (s==4) 2132 return _Beta(makesequence(v[0],v[1],v[3],1),contextptr)-_Beta(makesequence(v[0],v[1],v[2],1),contextptr); 2133 return gensizeerr(contextptr); 2134 } 2135 static const char _betad_cdf_s []="betad_cdf"; 2136 static define_unary_function_eval (__betad_cdf,&_betad_cdf,_betad_cdf_s); 2137 define_unary_function_ptr5( at_betad_cdf ,alias_at_betad_cdf,&__betad_cdf,0,true); 2138 betad_icdf(const gen & alpha_orig,const gen & beta_orig,const gen & t_orig,GIAC_CONTEXT)2139 gen betad_icdf(const gen &alpha_orig,const gen & beta_orig,const gen & t_orig,GIAC_CONTEXT){ 2140 if (is_zero(t_orig)|| is_one(t_orig)) 2141 return t_orig; 2142 gen t=evalf_double(t_orig,1,contextptr); 2143 gen alpha=evalf_double(alpha_orig,1,contextptr); 2144 gen beta=evalf_double(beta_orig,1,contextptr); 2145 if (alpha.type!=_DOUBLE_ || beta.type!=_DOUBLE_ || t.type!=_DOUBLE_ || alpha._DOUBLE_val<=0 || beta._DOUBLE_val<=0 || t._DOUBLE_val<0 || t._DOUBLE_val>1) 2146 return gensizeerr(contextptr); // symbolic(at_betad_icdf,makesequence(alpha_orig,beta_orig,t_orig)); 2147 double y=t._DOUBLE_val; 2148 if (y<=1e-13){ 2149 *logptr(contextptr) << "Underflow to 0" << "\n"; 2150 return 0; 2151 } 2152 if (y>=1-1e-13){ 2153 *logptr(contextptr) << "Overflow to 1" << "\n"; 2154 return 1; 2155 } 2156 // Initial guess 2157 double x0=.5; 2158 double prefactor=.5; 2159 if (alpha._DOUBLE_val>1){ 2160 if (beta._DOUBLE_val>1){ 2161 x0=(alpha._DOUBLE_val-1)/(alpha._DOUBLE_val+beta._DOUBLE_val-2); 2162 prefactor=1.; 2163 } 2164 else 2165 return 1-betad_icdf(beta,alpha,1-y,contextptr); 2166 } 2167 else { 2168 gen tmp; 2169 if (beta._DOUBLE_val<1 && y>.5) 2170 return 1-betad_icdf(beta,alpha,1-y,contextptr); 2171 double Y=y*Beta(alpha,beta,contextptr)._DOUBLE_val; 2172 tmp=exp(ln(alpha*gen(Y),contextptr)/alpha,contextptr); 2173 tmp=tmp*(1+tmp*gen(beta._DOUBLE_val-1)/(alpha._DOUBLE_val+1)); 2174 if (tmp.type==_DOUBLE_ && tmp._DOUBLE_val>0) 2175 x0=tmp._DOUBLE_val; 2176 if (x0<1e-4) 2177 return x0; 2178 } 2179 identificateur x(" x"); 2180 return newton(symbolic(at_Beta,makesequence(alpha,beta,x,1))-y,x,x0,NEWTON_DEFAULT_ITERATION,1e-5,1e-12,true,1,0,1,0,prefactor,contextptr); 2181 } _betad_icdf(const gen & g,GIAC_CONTEXT)2182 gen _betad_icdf(const gen & g,GIAC_CONTEXT){ 2183 if ( g.type==_STRNG && g.subtype==-1) return g; 2184 if (g.type!=_VECT) 2185 return gensizeerr(contextptr); 2186 vecteur & v=*g._VECTptr; 2187 int s=int(v.size()); 2188 if (s==3) 2189 return betad_icdf(v[0],v[1],v[2],contextptr); 2190 if (s==4) 2191 return betad_icdf(v[0],v[1],v[3],contextptr)-betad_icdf(v[0],v[1],v[2],contextptr); 2192 return gensizeerr(contextptr); 2193 } 2194 static const char _betad_icdf_s []="betad_icdf"; 2195 static define_unary_function_eval (__betad_icdf,&_betad_icdf,_betad_icdf_s); 2196 define_unary_function_ptr5( at_betad_icdf ,alias_at_betad_icdf,&__betad_icdf,0,true); 2197 _randbetad(const gen & args,GIAC_CONTEXT)2198 gen _randbetad(const gen & args,GIAC_CONTEXT){ 2199 if (args.type==_STRNG && args.subtype==-1) return args; 2200 if (args.type!=_VECT || args._VECTptr->size()!=2) 2201 return gensizeerr(contextptr); 2202 gen alpha=args._VECTptr->front(); 2203 gen beta=args._VECTptr->back(); 2204 alpha=evalf_double(alpha,1,contextptr); 2205 beta=evalf_double(beta,1,contextptr); 2206 if (is_positive(-alpha,contextptr) || is_positive(-beta,contextptr) || alpha.type!=_DOUBLE_ || beta.type!=_DOUBLE_) 2207 return gensizeerr(contextptr); 2208 double X=rgamma(alpha._DOUBLE_val,1.0,contextptr); 2209 double Y=rgamma(beta._DOUBLE_val,1.0,contextptr); 2210 return X/(X+Y); 2211 } 2212 static const char _randbetad_s []="randbetad"; 2213 static define_unary_function_eval (__randbetad,&_randbetad,_randbetad_s); 2214 define_unary_function_ptr5( at_randbetad ,alias_at_randbetad,&__randbetad,0,true); 2215 2216 static const char _betavariate_s []="betavariate"; 2217 static define_unary_function_eval (__betavariate,&_randbetad,_betavariate_s); 2218 define_unary_function_ptr5( at_betavariate ,alias_at_betavariate,&__betavariate,0,true); 2219 gammad(const gen & alpha,const gen & beta,const gen & x,GIAC_CONTEXT)2220 gen gammad(const gen &alpha,const gen & beta,const gen & x,GIAC_CONTEXT){ 2221 if (is_zero(x) && alpha==1) 2222 return beta; 2223 if (x==plus_inf) 2224 return 0; 2225 return pow(x,alpha-1,contextptr)*exp(-beta*x,contextptr)*pow(beta,alpha,contextptr)/Gamma(alpha,contextptr); 2226 } _gammad(const gen & g,GIAC_CONTEXT)2227 gen _gammad(const gen & g,GIAC_CONTEXT){ 2228 if ( g.type==_STRNG && g.subtype==-1) return g; 2229 if (g.type!=_VECT) 2230 return gammad(0,1,g,contextptr); 2231 vecteur & v=*g._VECTptr; 2232 int s=int(v.size()); 2233 if (s==2) 2234 return symbolic(at_gammad,g); 2235 if (s==3) 2236 return gammad(v[0],v[1],v[2],contextptr); 2237 return gensizeerr(contextptr); 2238 } 2239 static const char _gammad_s []="gammad"; 2240 static define_unary_function_eval (__gammad,&_gammad,_gammad_s); 2241 define_unary_function_ptr5( at_gammad ,alias_at_gammad,&__gammad,0,true); 2242 2243 // beta_cdf=Gamma regularized _gammad_cdf(const gen & g,GIAC_CONTEXT)2244 gen _gammad_cdf(const gen & g,GIAC_CONTEXT){ 2245 if ( g.type==_STRNG && g.subtype==-1) return g; 2246 if (g.type!=_VECT) 2247 return gensizeerr(contextptr); 2248 vecteur & v=*g._VECTptr; 2249 int s=int(v.size()); 2250 if (s==3) 2251 return _lower_incomplete_gamma(makesequence(v[0],v[1]*v[2],1),contextptr); 2252 if (s==4) 2253 return _lower_incomplete_gamma(makesequence(v[0],v[1]*v[3],1),contextptr)-_lower_incomplete_gamma(makesequence(v[0],v[1]*v[2],1),contextptr); 2254 return gensizeerr(contextptr); 2255 } 2256 static const char _gammad_cdf_s []="gammad_cdf"; 2257 static define_unary_function_eval (__gammad_cdf,&_gammad_cdf,_gammad_cdf_s); 2258 define_unary_function_ptr5( at_gammad_cdf ,alias_at_gammad_cdf,&__gammad_cdf,0,true); 2259 gammad_icdf(const gen & alpha_orig,const gen & beta_orig,const gen & t_orig,GIAC_CONTEXT)2260 gen gammad_icdf(const gen &alpha_orig,const gen & beta_orig,const gen & t_orig,GIAC_CONTEXT){ 2261 if (is_zero(t_orig)|| is_one(t_orig)) 2262 return t_orig; 2263 gen t=evalf_double(t_orig,1,contextptr); 2264 gen alpha=evalf_double(alpha_orig,1,contextptr); 2265 gen beta=evalf_double(beta_orig,1,contextptr); 2266 if (alpha.type!=_DOUBLE_ || beta.type!=_DOUBLE_ || t.type!=_DOUBLE_ || alpha._DOUBLE_val<=0 || beta._DOUBLE_val<=0 || t._DOUBLE_val<0 || t._DOUBLE_val>1) 2267 return gensizeerr(contextptr); // symbolic(at_gammad_icdf,makesequence(alpha_orig,beta_orig,t_orig)); 2268 double y=t._DOUBLE_val; 2269 if (y<=1e-13){ 2270 *logptr(contextptr) << "Underflow" << "\n"; 2271 return 0; 2272 } 2273 if (y>=1-1e-13){ 2274 *logptr(contextptr) << "Overflow" << "\n"; 2275 return plus_inf; 2276 } 2277 identificateur x(" x"); 2278 double x0=.5; // FIXME improve for y near boundaries! 2279 double prefactor=.5; 2280 if (alpha._DOUBLE_val>1){ 2281 x0=alpha._DOUBLE_val-1; 2282 prefactor=1.; 2283 } 2284 else { 2285 gen tmp=exp(ln(alpha*gen(y)*Gamma(alpha,contextptr),contextptr)/alpha,contextptr); 2286 tmp=tmp*(1-tmp/(alpha._DOUBLE_val+1)); 2287 if (tmp.type==_DOUBLE_ && tmp._DOUBLE_val>0) 2288 x0=tmp._DOUBLE_val; 2289 if (x0<1e-4) 2290 return x0; 2291 } 2292 return newton(symbolic(at_lower_incomplete_gamma,makesequence(alpha,x))-y*Gamma(alpha,contextptr),x,x0,NEWTON_DEFAULT_ITERATION,1e-5,1e-12,true,1,0,1,0,prefactor,contextptr)/beta; 2293 } _gammad_icdf(const gen & g,GIAC_CONTEXT)2294 gen _gammad_icdf(const gen & g,GIAC_CONTEXT){ 2295 if ( g.type==_STRNG && g.subtype==-1) return g; 2296 if (g.type!=_VECT) 2297 return gensizeerr(contextptr); 2298 vecteur & v=*g._VECTptr; 2299 int s=int(v.size()); 2300 if (s==3) 2301 return gammad_icdf(v[0],v[1],v[2],contextptr); 2302 if (s==4) 2303 return gammad_icdf(v[0],v[1],v[3],contextptr)-gammad_icdf(v[0],v[1],v[2],contextptr); 2304 return gensizeerr(contextptr); 2305 } 2306 static const char _gammad_icdf_s []="gammad_icdf"; 2307 static define_unary_function_eval (__gammad_icdf,&_gammad_icdf,_gammad_icdf_s); 2308 define_unary_function_ptr5( at_gammad_icdf ,alias_at_gammad_icdf,&__gammad_icdf,0,true); 2309 2310 #ifdef USE_GMP_REPLACEMENTS 2311 // a must be an integer for non GPL projects (https://en.wikipedia.org/wiki/Gamma_distribution#Generating_gamma-distributed_random_variables) rgamma(double a,double scale,GIAC_CONTEXT)2312 double rgamma(double a, double scale,GIAC_CONTEXT){ 2313 int n=a; 2314 if (a!=n) 2315 return 0.0/(n-n); // 0.0/0.0 chokes on visualc 2316 double res=0; 2317 for (int i=0;i<n;++i) 2318 res -= std::log(unif_rand(contextptr)); 2319 return res*scale; 2320 } 2321 2322 #else 2323 // modified from https://svn.r-project.org/R/trunk/src/nmath/rgamma.c, license GPL 2 or >= 2324 // assumes a>0, scale>0 2325 // thanks to R Rossignol for the reference rgamma(double a,double scale,GIAC_CONTEXT)2326 double rgamma(double a, double scale,GIAC_CONTEXT){ 2327 /* Constants : */ 2328 const static double sqrt32 = 5.656854; 2329 const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */ 2330 2331 /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k)) 2332 * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k) 2333 * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k) 2334 */ 2335 const static double q1 = 0.04166669; 2336 const static double q2 = 0.02083148; 2337 const static double q3 = 0.00801191; 2338 const static double q4 = 0.00144121; 2339 const static double q5 = -7.388e-5; 2340 const static double q6 = 2.4511e-4; 2341 const static double q7 = 2.424e-4; 2342 2343 const static double a1 = 0.3333333; 2344 const static double a2 = -0.250003; 2345 const static double a3 = 0.2000062; 2346 const static double a4 = -0.1662921; 2347 const static double a5 = 0.1423657; 2348 const static double a6 = -0.1367177; 2349 const static double a7 = 0.1233795; 2350 2351 /* State variables [FIXME for threading!] :*/ 2352 static double aa = 0.; 2353 static double aaa = 0.; 2354 static double s, s2, d; /* no. 1 (step 1) */ 2355 static double q0, b, si, c;/* no. 2 (step 4) */ 2356 2357 double e, p, q, r, t, u, v, w, x, ret_val; 2358 2359 // if(!R_FINITE(a) || !R_FINITE(scale)) return ML_POSINF; 2360 2361 if (a < 1.) { /* GS algorithm for parameters a < 1 */ 2362 e = 1.0 + exp_m1 * a; 2363 for (;;) { 2364 p = e * unif_rand(contextptr); 2365 if (p >= 1.0) { 2366 x = -std::log((e - p) / a); 2367 if (exp_rand(contextptr) >= (1.0 - a) * std::log(x)) 2368 break; 2369 } else { 2370 x = std::exp(std::log(p) / a); 2371 if (exp_rand(contextptr) >= x) 2372 break; 2373 } 2374 } 2375 return scale * x; 2376 } 2377 2378 /* --- a >= 1 : GD algorithm --- */ 2379 2380 /* Step 1: Recalculations of s2, s, d if a has changed */ 2381 if (a != aa) { 2382 aa = a; 2383 s2 = a - 0.5; 2384 s = std::sqrt(s2); 2385 d = sqrt32 - s * 12.0; 2386 } 2387 /* Step 2: t = standard normal deviate, 2388 x = (s,1/2) -normal deviate. */ 2389 2390 /* immediate acceptance (i) */ 2391 t = randNorm(contextptr); // norm_rand(); 2392 x = s + 0.5 * t; 2393 ret_val = x * x; 2394 if (t >= 0.0) 2395 return scale * ret_val; 2396 2397 /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */ 2398 u = unif_rand(contextptr); 2399 if (d * u <= t * t * t) 2400 return scale * ret_val; 2401 2402 /* Step 4: recalculations of q0, b, si, c if necessary */ 2403 2404 if (a != aaa) { 2405 aaa = a; 2406 r = 1.0 / a; 2407 q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r 2408 + q2) * r + q1) * r; 2409 2410 /* Approximation depending on size of parameter a */ 2411 /* The constants in the expressions for b, si and c */ 2412 /* were established by numerical experiments */ 2413 2414 if (a <= 3.686) { 2415 b = 0.463 + s + 0.178 * s2; 2416 si = 1.235; 2417 c = 0.195 / s - 0.079 + 0.16 * s; 2418 } else if (a <= 13.022) { 2419 b = 1.654 + 0.0076 * s2; 2420 si = 1.68 / s + 0.275; 2421 c = 0.062 / s + 0.024; 2422 } else { 2423 b = 1.77; 2424 si = 0.75; 2425 c = 0.1515 / s; 2426 } 2427 } 2428 /* Step 5: no quotient test if x not positive */ 2429 2430 if (x > 0.0) { 2431 /* Step 6: calculation of v and quotient q */ 2432 v = t / (s + s); 2433 if (fabs(v) <= 0.25) 2434 q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v 2435 + a3) * v + a2) * v + a1) * v; 2436 else 2437 q = q0 - s * t + 0.25 * t * t + (s2 + s2) * std::log(1.0 + v); 2438 2439 2440 /* Step 7: quotient acceptance (q) */ 2441 if (std::log(1.0 - u) <= q) 2442 return scale * ret_val; 2443 } 2444 2445 for(;;) { 2446 /* Step 8: e = standard exponential deviate 2447 * u = 0,1 -uniform deviate 2448 * t = (b,si)-double exponential (laplace) sample */ 2449 e = exp_rand(contextptr); 2450 u = unif_rand(contextptr); 2451 u = u + u - 1.0; 2452 if (u < 0.0) 2453 t = b - si * e; 2454 else 2455 t = b + si * e; 2456 /* Step 9: rejection if t < tau(1) = -0.71874483771719 */ 2457 if (t >= -0.71874483771719) { 2458 /* Step 10: calculation of v and quotient q */ 2459 v = t / (s + s); 2460 if (fabs(v) <= 0.25) 2461 q = q0 + 0.5 * t * t * 2462 ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v 2463 + a2) * v + a1) * v; 2464 else 2465 q = q0 - s * t + 0.25 * t * t + (s2 + s2) * std::log(1.0 + v); 2466 /* Step 11: hat acceptance (h) */ 2467 /* (if q not positive go to step 8) */ 2468 if (q > 0.0) { 2469 w = expm1(q); 2470 /* ^^^^^ original code had approximation with rel.err < 2e-7 */ 2471 /* if t is rejected sample again at step 8 */ 2472 if (c * fabs(u) <= w * std::exp(e - 0.5 * t * t)) 2473 break; 2474 } 2475 } 2476 } /* repeat .. until `t' is accepted */ 2477 x = s + 0.5 * t; 2478 return scale * x * x; 2479 } 2480 #endif 2481 2482 // args=(alpha,beta) _randgammad(const gen & args,GIAC_CONTEXT)2483 gen _randgammad(const gen & args,GIAC_CONTEXT){ 2484 if (args.type==_STRNG && args.subtype==-1) return args; 2485 if (args.type!=_VECT || args._VECTptr->size()!=2) 2486 return gensizeerr(contextptr); 2487 gen alpha=args._VECTptr->front(); 2488 gen beta=args._VECTptr->back(); 2489 alpha=evalf_double(alpha,1,contextptr); 2490 beta=evalf_double(beta,1,contextptr); 2491 if (is_positive(-alpha,contextptr) || is_positive(-beta,contextptr) || alpha.type!=_DOUBLE_ || beta.type!=_DOUBLE_) 2492 return gensizeerr(contextptr); 2493 return rgamma(alpha._DOUBLE_val,1.0/beta._DOUBLE_val,contextptr); 2494 } 2495 static const char _randgammad_s []="randgammad"; 2496 static define_unary_function_eval (__randgammad,&_randgammad,_randgammad_s); 2497 define_unary_function_ptr5( at_randgammad ,alias_at_randgammad,&__randgammad,0,true); 2498 2499 static const char _gammavariate_s []="gammavariate"; 2500 static define_unary_function_eval (__gammavariate,&_randgammad,_gammavariate_s); 2501 define_unary_function_ptr5( at_gammavariate ,alias_at_gammavariate,&__gammavariate,0,true); 2502 uniform(const gen & g,bool ckpython,GIAC_CONTEXT)2503 gen uniform(const gen & g,bool ckpython,GIAC_CONTEXT){ 2504 if ( g.type==_STRNG && g.subtype==-1) return g; 2505 if (g.type!=_VECT) 2506 return 1; 2507 vecteur & v=*g._VECTptr; 2508 int s=int(v.size()); 2509 if (s==0) 2510 return symbolic(at_uniformd,makesequence(0,1)); 2511 if (s==2){ 2512 if (ckpython) 2513 return v[0]+(giac_rand(contextptr)/(rand_max2+1.0))*(v[1]-v[0]); 2514 return symbolic(at_uniformd,makesequence(v[0],v[1])); 2515 } 2516 if (s==3) 2517 return inv(v[1]-v[0],contextptr); 2518 return gensizeerr(contextptr); 2519 } _uniform(const gen & g,GIAC_CONTEXT)2520 gen _uniform(const gen & g,GIAC_CONTEXT){ 2521 return uniform(g,true,contextptr); 2522 } 2523 static const char _uniform_s []="uniform"; 2524 static define_unary_function_eval (__uniform,&_uniform,_uniform_s); 2525 define_unary_function_ptr5( at_uniform ,alias_at_uniform,&__uniform,0,true); _uniformd(const gen & g,GIAC_CONTEXT)2526 gen _uniformd(const gen & g,GIAC_CONTEXT){ 2527 return uniform(g,false,contextptr); 2528 } 2529 static const char _uniformd_s []="uniformd"; 2530 static define_unary_function_eval (__uniformd,&_uniformd,_uniformd_s); 2531 define_unary_function_ptr5( at_uniformd ,alias_at_uniformd,&__uniformd,0,true); 2532 _uniform_cdf(const gen & g,GIAC_CONTEXT)2533 gen _uniform_cdf(const gen & g,GIAC_CONTEXT){ 2534 if ( g.type==_STRNG && g.subtype==-1) return g; 2535 if (g.type!=_VECT) 2536 return g; 2537 vecteur & v=*g._VECTptr; 2538 int s=int(v.size()); 2539 if (s==3) 2540 return (v[2]-v[0])/(v[1]-v[0]); 2541 if (s==4) 2542 return (v[3]-v[2])/(v[1]-v[0]); 2543 return gensizeerr(contextptr); 2544 } 2545 static const char _uniform_cdf_s []="uniform_cdf"; 2546 static define_unary_function_eval (__uniform_cdf,&_uniform_cdf,_uniform_cdf_s); 2547 define_unary_function_ptr5( at_uniform_cdf ,alias_at_uniform_cdf,&__uniform_cdf,0,true); 2548 static const char _uniformd_cdf_s []="uniformd_cdf"; 2549 static define_unary_function_eval (__uniformd_cdf,&_uniform_cdf,_uniformd_cdf_s); 2550 define_unary_function_ptr5( at_uniformd_cdf ,alias_at_uniformd_cdf,&__uniformd_cdf,0,true); 2551 _uniform_icdf(const gen & g,GIAC_CONTEXT)2552 gen _uniform_icdf(const gen & g,GIAC_CONTEXT){ 2553 if ( g.type==_STRNG && g.subtype==-1) return g; 2554 if (g.type!=_VECT) 2555 return g; 2556 vecteur & v=*g._VECTptr; 2557 int s=int(v.size()); 2558 if (s==3) 2559 return v[0]+v[2]*(v[1]-v[0]); 2560 if (s==4) 2561 return (v[3]-v[2])*(v[1]-v[0]); 2562 return gensizeerr(contextptr); 2563 } 2564 static const char _uniform_icdf_s []="uniform_icdf"; 2565 static define_unary_function_eval (__uniform_icdf,&_uniform_icdf,_uniform_icdf_s); 2566 define_unary_function_ptr5( at_uniform_icdf ,alias_at_uniform_icdf,&__uniform_icdf,0,true); 2567 static const char _uniformd_icdf_s []="uniformd_icdf"; 2568 static define_unary_function_eval (__uniformd_icdf,&_uniform_icdf,_uniformd_icdf_s); 2569 define_unary_function_ptr5( at_uniformd_icdf ,alias_at_uniformd_icdf,&__uniformd_icdf,0,true); 2570 _exponential(const gen & g,GIAC_CONTEXT)2571 gen _exponential(const gen & g,GIAC_CONTEXT){ 2572 if ( g.type==_STRNG && g.subtype==-1) return g; 2573 if (g.type!=_VECT) 2574 return symbolic(at_exponentiald,g); 2575 vecteur & v=*g._VECTptr; 2576 int s=int(v.size()); 2577 if (s==2) 2578 return v[0]*exp(-v[0]*v[1],contextptr); 2579 return gensizeerr(contextptr); 2580 } 2581 static const char _exponential_s []="exponential"; 2582 static define_unary_function_eval (__exponential,&_exponential,_exponential_s); 2583 define_unary_function_ptr5( at_exponential ,alias_at_exponential,&__exponential,0,true); 2584 static const char _exponentiald_s []="exponentiald"; 2585 static define_unary_function_eval (__exponentiald,&_exponential,_exponentiald_s); 2586 define_unary_function_ptr5( at_exponentiald ,alias_at_exponentiald,&__exponentiald,0,true); 2587 _exponential_cdf(const gen & g,GIAC_CONTEXT)2588 gen _exponential_cdf(const gen & g,GIAC_CONTEXT){ 2589 if ( g.type==_STRNG && g.subtype==-1) return g; 2590 if (g.type!=_VECT) 2591 return gensizeerr(contextptr); 2592 vecteur & v=*g._VECTptr; 2593 int s=int(v.size()); 2594 if (s==2) 2595 return 1-exp(-v[0]*v[1],contextptr); 2596 if (s==3) 2597 return exp(-v[0]*v[1],contextptr)-exp(-v[0]*v[2],contextptr); 2598 return gensizeerr(contextptr); 2599 } 2600 static const char _exponential_cdf_s []="exponential_cdf"; 2601 static define_unary_function_eval (__exponential_cdf,&_exponential_cdf,_exponential_cdf_s); 2602 define_unary_function_ptr5( at_exponential_cdf ,alias_at_exponential_cdf,&__exponential_cdf,0,true); 2603 static const char _exponentiald_cdf_s []="exponentiald_cdf"; 2604 static define_unary_function_eval (__exponentiald_cdf,&_exponential_cdf,_exponentiald_cdf_s); 2605 define_unary_function_ptr5( at_exponentiald_cdf ,alias_at_exponentiald_cdf,&__exponentiald_cdf,0,true); 2606 2607 // 1 - exp(-l*x)=p => exp(-l*x)=1-p => x=-1/l*ln(1-p) _exponential_icdf(const gen & g,GIAC_CONTEXT)2608 gen _exponential_icdf(const gen & g,GIAC_CONTEXT){ 2609 if ( g.type==_STRNG && g.subtype==-1) return g; 2610 if (g.type!=_VECT) 2611 return gensizeerr(contextptr); 2612 vecteur & v=*g._VECTptr; 2613 int s=int(v.size()); 2614 if (s==2) 2615 return -ln(1-v[1],contextptr)/v[0]; 2616 if (s==3) 2617 return (ln(1-v[1],contextptr)-ln(1-v[2],contextptr))/v[0]; 2618 return gensizeerr(contextptr); 2619 } 2620 static const char _exponential_icdf_s []="exponential_icdf"; 2621 static define_unary_function_eval (__exponential_icdf,&_exponential_icdf,_exponential_icdf_s); 2622 define_unary_function_ptr5( at_exponential_icdf ,alias_at_exponential_icdf,&__exponential_icdf,0,true); 2623 static const char _exponentiald_icdf_s []="exponentiald_icdf"; 2624 static define_unary_function_eval (__exponentiald_icdf,&_exponential_icdf,_exponentiald_icdf_s); 2625 define_unary_function_ptr5( at_exponentiald_icdf ,alias_at_exponentiald_icdf,&__exponentiald_icdf,0,true); 2626 2627 // geometric(p,k)=(1-p)^(k-1)*p geometric(const gen & p,const gen & k,GIAC_CONTEXT)2628 gen geometric(const gen & p,const gen & k,GIAC_CONTEXT){ 2629 gen K(k); 2630 if (is_positive(-k,contextptr)) 2631 return gensizeerr(contextptr); 2632 return pow(1-p,k-1,contextptr)*p; 2633 } _geometric(const gen & g,GIAC_CONTEXT)2634 gen _geometric(const gen & g,GIAC_CONTEXT){ 2635 if ( g.type==_STRNG && g.subtype==-1) return g; 2636 if (g.type!=_VECT) 2637 return symbolic(at_geometric,g); 2638 vecteur & v=*g._VECTptr; 2639 int s=int(v.size()); 2640 if (s==2) 2641 return geometric(v[0],v[1],contextptr); 2642 return gensizeerr(contextptr); 2643 } 2644 static const char _geometric_s []="geometric"; 2645 static define_unary_function_eval (__geometric,&_geometric,_geometric_s); 2646 define_unary_function_ptr5( at_geometric ,alias_at_geometric,&__geometric,0,true); 2647 2648 // geometric_cdf(p,k)=1-(1-p)^k geometric_cdf(const gen & p,const gen & k,GIAC_CONTEXT)2649 gen geometric_cdf(const gen & p,const gen & k,GIAC_CONTEXT){ 2650 gen K(k); 2651 if (is_strictly_positive(-k,contextptr)) 2652 return gensizeerr(contextptr); 2653 return 1-pow(1-p,k,contextptr); 2654 } _geometric_cdf(const gen & g,GIAC_CONTEXT)2655 gen _geometric_cdf(const gen & g,GIAC_CONTEXT){ 2656 if ( g.type==_STRNG && g.subtype==-1) return g; 2657 if (g.type!=_VECT) 2658 return symbolic(at_geometric_cdf,g); 2659 vecteur & v=*g._VECTptr; 2660 int s=int(v.size()); 2661 if (s==2) 2662 return geometric_cdf(v[0],v[1],contextptr); 2663 if (s==3) 2664 return geometric_cdf(v[0],v[2],contextptr)-geometric_cdf(v[0],v[1]-1,contextptr); 2665 return gensizeerr(contextptr); 2666 } 2667 static const char _geometric_cdf_s []="geometric_cdf"; 2668 static define_unary_function_eval (__geometric_cdf,&_geometric_cdf,_geometric_cdf_s); 2669 define_unary_function_ptr5( at_geometric_cdf ,alias_at_geometric_cdf,&__geometric_cdf,0,true); 2670 2671 // k=geometric_icdf(p,P) if 1-(1-p)^k>=P hence geometric_icdf(const gen & p,const gen & P,GIAC_CONTEXT)2672 gen geometric_icdf(const gen & p,const gen & P,GIAC_CONTEXT){ 2673 return _ceil(ln(1-P,contextptr)/ln(1-p,contextptr),contextptr); 2674 } _geometric_icdf(const gen & g,GIAC_CONTEXT)2675 gen _geometric_icdf(const gen & g,GIAC_CONTEXT){ 2676 if ( g.type==_STRNG && g.subtype==-1) return g; 2677 if (g.type!=_VECT) 2678 return symbolic(at_geometric_icdf,g); 2679 vecteur & v=*g._VECTptr; 2680 int s=int(v.size()); 2681 if (s==2) 2682 return geometric_icdf(v[0],v[1],contextptr); 2683 if (s==3) 2684 return geometric_icdf(v[0],v[2],contextptr)-geometric_icdf(v[0],v[1],contextptr); 2685 return gensizeerr(contextptr); 2686 } 2687 static const char _geometric_icdf_s []="geometric_icdf"; 2688 static define_unary_function_eval (__geometric_icdf,&_geometric_icdf,_geometric_icdf_s); 2689 define_unary_function_ptr5( at_geometric_icdf ,alias_at_geometric_icdf,&__geometric_icdf,0,true); 2690 _randgeometric(const gen & g,GIAC_CONTEXT)2691 gen _randgeometric(const gen & g,GIAC_CONTEXT){ 2692 if ( g.type==_STRNG && g.subtype==-1) return g; 2693 return _ceil(gen(std::log(1-giac_rand(contextptr)/(rand_max2+1.0)))/ln(1-g,contextptr),contextptr); 2694 } 2695 static const char _randgeometric_s []="randgeometric"; 2696 static define_unary_function_eval (__randgeometric,&_randgeometric,_randgeometric_s); 2697 define_unary_function_ptr5( at_randgeometric ,alias_at_randgeometric,&__randgeometric,0,true); 2698 kolmogorovd(double c2)2699 double kolmogorovd(double c2){ 2700 c2=c2*c2; 2701 // 2*sum((-1)^(r-1)*exp(-2*r^2*c^2),r,1,inf) 2702 long_double cumul=0; 2703 for (int r=1;;++r){ 2704 long_double current=std::exp(-2*r*r*c2); 2705 if (cumul==cumul+current) 2706 return 1-double(2*cumul); 2707 if (r%2) 2708 cumul += current; 2709 else 2710 cumul -= current; 2711 } 2712 } 2713 _kolmogorovd(const gen & g,GIAC_CONTEXT)2714 gen _kolmogorovd(const gen & g,GIAC_CONTEXT){ 2715 if ( g.type==_STRNG && g.subtype==-1) return g; 2716 if (g.type==_VECT) 2717 return apply(g,_kolmogorovd,contextptr); 2718 gen tmp=evalf_double(g,1,contextptr); 2719 if (tmp.type!=_DOUBLE_) 2720 return symbolic(at_kolmogorovd,g); 2721 if (is_positive(-g,contextptr)) 2722 return undef; 2723 double c2=tmp._DOUBLE_val; 2724 return kolmogorovd(c2); 2725 } 2726 static const char _kolmogorovd_s []="kolmogorovd"; 2727 static define_unary_function_eval (__kolmogorovd,&_kolmogorovd,_kolmogorovd_s); 2728 define_unary_function_ptr5( at_kolmogorovd ,alias_at_kolmogorovd,&__kolmogorovd,0,true); 2729 2730 /* 2731 gen _kolmogorov_cdf(const gen & g,GIAC_CONTEXT){ 2732 if ( g.type==_STRNG && g.subtype==-1) return g; 2733 if (g.type==_VECT) 2734 return apply(g,_kolmogorov_cdf,contextptr); 2735 gen tmp=evalf_double(g,1,contextptr); 2736 if (tmp.type!=_DOUBLE_) 2737 return symbolic(at_kolmogorov_cdf,g); 2738 if (is_positive(-g,contextptr)) 2739 return undef; 2740 double c=tmp._DOUBLE_val,sqrt2=std::sqrt(2.0); 2741 // 1-1/ln(2)*sum((-1)^(r-1)/r*erfc(sqrt(2)*r*c),r,1,inf) 2742 double cumul=0; 2743 for (int r=1;;++r){ 2744 double current=_erfc(sqrt2*r*c,contextptr)._DOUBLE_val/r; 2745 if (cumul==cumul+current) 2746 return 1-cumul/std::log(2.0); 2747 if (r%2) 2748 cumul += current; 2749 else 2750 cumul -= current; 2751 } 2752 } 2753 static const char _kolmogorov_cdf_s []="kolmogorov_cdf"; 2754 static define_unary_function_eval (__kolmogorov_cdf,&_kolmogorov_cdf,_kolmogorov_cdf_s); 2755 define_unary_function_ptr5( at_kolmogorov_cdf ,alias_at_kolmogorov_cdf,&__kolmogorov_cdf,0,true); 2756 */ is_discrete_distribution(int nd)2757 bool is_discrete_distribution(int nd){ 2758 return nd==2 || nd==3 || nd==4 || nd==12; 2759 } 2760 _kolmogorovt(const gen & g,GIAC_CONTEXT)2761 gen _kolmogorovt(const gen & g,GIAC_CONTEXT){ 2762 if ( g.type==_STRNG && g.subtype==-1) return g; 2763 if (g.type!=_VECT || g._VECTptr->size()<2) 2764 return gensizeerr(contextptr); 2765 vecteur & v = *g._VECTptr; 2766 gen x=v[0],y=v[1]; 2767 if (y==at_exp) 2768 y=at_exponential; 2769 if (v.size()==3) 2770 return _kolmogorovt(makesequence(x,y(v[2],contextptr)),contextptr); 2771 if (v.size()>2) 2772 return _kolmogorovt(makesequence(x,y(vecteur(v.begin()+2,v.end()),contextptr)),contextptr); 2773 if (is_distribution(x)){ 2774 if (is_distribution(y)) 2775 return gensizeerr(contextptr); 2776 swapgen(x,y); 2777 } 2778 int nd=is_distribution(y); 2779 if (nd){ 2780 if (y==at_normald || y==at_normal || y==at_cauchyd || y==at_cauchy) 2781 y=symbolic(*y._FUNCptr,makesequence(0,1)); 2782 if (y.type!=_SYMB) 2783 return gensizeerr(contextptr); 2784 } 2785 vector<giac_double> X; 2786 x=evalf_double(x,1,contextptr); 2787 if (x.type!=_VECT || !convert(*x._VECTptr,X,true)) 2788 return gensizeerr(contextptr); 2789 sort(X.begin(),X.end()); 2790 int n=int(X.size()); 2791 double cumulx=0,d=0,dcur,invn=1./n,ks; 2792 if (nd){ 2793 gen f=cdf(nd); 2794 if (f.type!=_FUNC) 2795 return gensizeerr(contextptr); 2796 // add parameters from y 2797 vecteur yf=gen2vecteur(y._SYMBptr->feuille); 2798 if (int(yf.size())!=distrib_nargs(nd)) 2799 return gensizeerr(contextptr); 2800 yf.push_back(0); 2801 f=symbolic(*f._FUNCptr,gen(yf,_SEQ__VECT)); 2802 gen & fback=f._SYMBptr->feuille._VECTptr->back(); 2803 if (is_discrete_distribution(nd)){ 2804 for (vector<giac_double>::const_iterator it=X.begin();it!=X.end();){ 2805 fback=double(*it)-1; 2806 gen prevtmp=evalf_double(f,1,contextptr); 2807 dcur=absdouble(prevtmp._DOUBLE_val-cumulx); 2808 if (dcur>d) 2809 d=dcur; 2810 fback=double(*it); 2811 gen tmp=evalf_double(f,1,contextptr); 2812 int i=1; 2813 for (++it;it!=X.end() && double(*it)==fback;++it) 2814 ++i; 2815 cumulx += i*invn; 2816 dcur=absdouble(tmp._DOUBLE_val-cumulx); 2817 if (dcur>d) 2818 d=dcur; 2819 } 2820 ks=d*std::sqrt(double(n)); 2821 return makevecteur(string2gen("D=",false),d,string2gen("K=",false),ks,string2gen("1-kolmogorovd(K)=",false),1-kolmogorovd(ks)); 2822 } 2823 for (int i=0;i<n;++i){ 2824 fback=double(X[i]); 2825 gen tmp=evalf_double(f,1,contextptr); 2826 if (tmp.type!=_DOUBLE_) 2827 return gensizeerr(contextptr); 2828 dcur=absdouble(tmp._DOUBLE_val-cumulx); 2829 if (dcur>d) 2830 d=dcur; 2831 cumulx += invn; 2832 dcur=absdouble(tmp._DOUBLE_val-cumulx); 2833 if (dcur>d) 2834 d=dcur; 2835 } 2836 ks=d*std::sqrt(double(n)); 2837 return makevecteur(string2gen("D=",false),d,string2gen("K=",false),ks,string2gen("1-kolmogorovd(K)=",false),1-kolmogorovd(ks)); 2838 } 2839 // 2 lists 2840 vector<giac_double> Y; 2841 y=evalf_double(y,1,contextptr); 2842 if (y.type!=_VECT || !convert(*y._VECTptr,Y,true)) 2843 return gensizeerr(contextptr); 2844 sort(Y.begin(),Y.end()); 2845 int m=int(Y.size()); 2846 double cumuly=0,invm=1./m; 2847 int i=0,j=0; 2848 while (i<n && j<m){ 2849 if (X[i]==Y[j]){ 2850 cumulx += invn; 2851 cumuly += invm; 2852 ++i; ++j; 2853 } 2854 else { 2855 if (X[i]>Y[j]){ 2856 cumuly += invm; 2857 ++j; 2858 } 2859 else { 2860 cumulx += invn; 2861 ++i; 2862 } 2863 } 2864 dcur=absdouble(cumulx-cumuly); 2865 if (dcur>d) 2866 d=dcur; 2867 } 2868 ks=d*std::sqrt((n*m)/double(n+m)); 2869 return makevecteur(string2gen("D=",false),d,string2gen("K=",false),ks,string2gen("1-kolmogorovd(K)=",false),1-kolmogorovd(ks)); 2870 } 2871 static const char _kolmogorovt_s []="kolmogorovt"; 2872 static define_unary_function_eval (__kolmogorovt,&_kolmogorovt,_kolmogorovt_s); 2873 define_unary_function_ptr5( at_kolmogorovt ,alias_at_kolmogorovt,&__kolmogorovt,0,true); 2874 2875 2876 // computes wilcoxon test value for sample x and median m_ or samples x and m_ wilcoxons(const vecteur & x,const gen & m_,GIAC_CONTEXT)2877 gen wilcoxons(const vecteur & x,const gen & m_,GIAC_CONTEXT){ 2878 if (m_.type==_VECT){ 2879 vecteur & y=*m_._VECTptr; 2880 int n=int(y.size()); 2881 int m=int(x.size()); 2882 vecteur xm; 2883 for (int i=0;i<n;++i){ 2884 xm.push_back(makevecteur(y[i],i)); 2885 } 2886 for (int i=0;i<m;++i){ 2887 xm.push_back(makevecteur(x[i],n+i)); 2888 } 2889 gen_sort_f(xm.begin(),xm.end(),first_ascend_sort); 2890 vector<double> stat(xm.size()); 2891 for (unsigned i=1;i<=xm.size();++i){ 2892 unsigned j=1; // number of ties 2893 for (;i<xm.size() && xm[i]._VECTptr->front()==xm[i-1]._VECTptr->front();){ 2894 ++j; ++i; 2895 } 2896 // xm[i]!=xm[i-1] (or i==xm.size()) 2897 // xm[i-1]==...=xm[i-j] 2898 double value=(i-1+i-j)/2.+1; 2899 for (unsigned k=1;k<=j;++k) 2900 stat[i-k]=value; 2901 } 2902 gen res=0; 2903 for (unsigned i=0;i<xm.size();++i){ 2904 if (is_strictly_greater(n,xm[i]._VECTptr->back().val,contextptr)) 2905 res += stat[i]; // int(i+1); 2906 } 2907 return res; 2908 } 2909 gen m=evalf_double(m_,1,contextptr); 2910 if (m.type!=_DOUBLE_) 2911 return gensizeerr(gettext("Invalid median")); 2912 vecteur xm(x); 2913 for (unsigned i=0;i<xm.size();++i){ 2914 xm[i]=makevecteur(abs(xm[i]-m),int(i)); 2915 } 2916 gen_sort_f(xm.begin(),xm.end(),first_ascend_sort); 2917 gen res=0; 2918 for (unsigned i=0;i<xm.size();++i){ 2919 if (is_greater(x[xm[i]._VECTptr->back().val],m,contextptr)) 2920 res += int(i+1); 2921 } 2922 return res; 2923 } 2924 _wilcoxons(const gen & args,GIAC_CONTEXT)2925 gen _wilcoxons(const gen & args,GIAC_CONTEXT){ 2926 if (args.type!=_VECT || args._VECTptr->size()!=2) 2927 return gensizeerr(contextptr); 2928 gen x=args._VECTptr->front(),m=args._VECTptr->back(); 2929 if (x.type!=_VECT || x._VECTptr->empty()) 2930 return gendimerr(contextptr); 2931 return wilcoxons(*x._VECTptr,m,contextptr); 2932 } 2933 static const char _wilcoxons_s []="wilcoxons"; 2934 static define_unary_function_eval (__wilcoxons,&_wilcoxons,_wilcoxons_s); 2935 define_unary_function_ptr5( at_wilcoxons ,alias_at_wilcoxons,&__wilcoxons,0,true); 2936 2937 // returns product_i=1^n (1+t^i) wilcoxonp(int n)2938 gen wilcoxonp(int n){ 2939 if (n<=0) 2940 return vecteur(0); 2941 gen res(gen(vecteur(1,1),_POLY1__VECT)); 2942 for (int i=1;i<=n;++i){ 2943 vecteur tmp(i+1); 2944 tmp[i]=tmp[0]=1; 2945 res=res*gen(tmp,_POLY1__VECT); 2946 } 2947 return res; 2948 } wilcoxonp(int m,int n,GIAC_CONTEXT)2949 gen wilcoxonp(int m,int n,GIAC_CONTEXT){ 2950 // sum_k Pmn(k)*x^k = 1/comb(m+n,n)*prod_{i=n+1}^{m+n}(1-x^i)/prod_{j=1}^m(1-x^j) 2951 if (n<=0 || m<=0) 2952 return vecteur(0); 2953 gen num(gen(vecteur(1,1),_POLY1__VECT)); 2954 for (int i=n+1;i<=m+n;++i){ 2955 vecteur tmp(i+1); 2956 tmp[i]=-1; tmp[0]=1; 2957 num=num*gen(tmp,_POLY1__VECT); 2958 } 2959 gen den(gen(vecteur(1,1),_POLY1__VECT)); 2960 for (int i=1;i<=m;++i){ 2961 vecteur tmp(i+1); 2962 tmp[i]=-1; tmp[0]=1; 2963 den=den*gen(tmp,_POLY1__VECT); 2964 } 2965 gen q=_quo(makesequence(num,den),contextptr); 2966 return q; 2967 } _wilcoxonp(const gen & args,GIAC_CONTEXT)2968 gen _wilcoxonp(const gen & args,GIAC_CONTEXT){ 2969 gen n(args); 2970 if (n.type==_VECT && n._VECTptr->size()==2){ 2971 gen M(n._VECTptr->front()),N(n._VECTptr->back()); 2972 if (!is_integral(M) || M.type!=_INT_ || M.val < 1 || 2973 !is_integral(N) || N.type!=_INT_ || N.val < 1 || M.val+N.val > 400 2974 ) 2975 return gendimerr(contextptr); 2976 return wilcoxonp(M.val,N.val,contextptr)/comb(M.val+N.val,N.val); 2977 } 2978 if (!is_integral(n) || n.type!=_INT_ || n.val<1 || n.val>1000) 2979 return gendimerr(contextptr); 2980 return wilcoxonp(n.val)/pow(plus_two,n,contextptr); 2981 } 2982 static const char _wilcoxonp_s []="wilcoxonp"; 2983 static define_unary_function_eval (__wilcoxonp,&_wilcoxonp,_wilcoxonp_s); 2984 define_unary_function_ptr5( at_wilcoxonp ,alias_at_wilcoxonp,&__wilcoxonp,0,true); 2985 2986 // a faire wilcoxont(echantillon,mediane,seuil alpha) renvoyant true (accepte 2987 // test 2-sided) si W est dans l'intervalle symetrique de borne inferieure 2988 // le plus grand k tel que wilcoxonp[0]+...+wilcoxonp[k-1]<alpha/2 _wilcoxont(const gen & g,GIAC_CONTEXT)2989 gen _wilcoxont(const gen & g,GIAC_CONTEXT){ 2990 if ( g.type==_STRNG && g.subtype==-1) return g; 2991 if (g.type!=_VECT || g._VECTptr->size()<2 || g._VECTptr->size()>4) 2992 return gensizeerr(contextptr); 2993 vecteur & v = *g._VECTptr; 2994 gen x=v[0],m=v[1],alpha=0.05; 2995 int typetest=0; // 1 >, -1 < 2996 if (x.type!=_VECT || x._VECTptr->empty()) 2997 return gensizeerr(gettext("Invalid sample")); 2998 if (v.size()>=3 && v[2].type!=_FUNC) 2999 alpha=evalf_double(v[2],1,contextptr); 3000 if (v.size()>=4 && v[3].type!=_FUNC) 3001 alpha=evalf_double(v[3],1,contextptr); 3002 if (v.size()>=3 && v[2].type==_FUNC){ 3003 if (v[2]==at_superieur_strict || v[2]==at_superieur_egal) 3004 typetest=1; 3005 if (v[2]==at_inferieur_strict || v[2]==at_inferieur_egal) 3006 typetest=-1; 3007 } 3008 if (v.size()>=4 && v[3].type==_FUNC){ 3009 if (v[3]==at_superieur_strict || v[3]==at_superieur_egal) 3010 typetest=1; 3011 if (v[3]==at_inferieur_strict || v[3]==at_inferieur_egal) 3012 typetest=-1; 3013 } 3014 if (alpha.type!=_DOUBLE_ || alpha._DOUBLE_val<=0 || alpha._DOUBLE_val>=1) 3015 return gensizeerr(gettext("Invalid confidence level")); 3016 int n=int(x._VECTptr->size()); 3017 if (m.type==_VECT){ 3018 // if (typetest!=0) return gensizeerr(gettext("H1 must be <> for Mann-Whitney test")); 3019 gen w=wilcoxons(*m._VECTptr,x,contextptr); 3020 if (w.type!=_DOUBLE_) 3021 return gensizeerr(contextptr); 3022 int N=int(m._VECTptr->size()),M=int(x._VECTptr->size()); 3023 gen combMN=comb(M+N,N); 3024 *logptr(contextptr) << "Mann-Whitney 2-sample test, H0 same Median, H1 "; 3025 if (typetest==0) *logptr(contextptr) << "<>"; 3026 if (typetest==1) *logptr(contextptr) << ">"; 3027 if (typetest==-1) *logptr(contextptr) << "<"; 3028 *logptr(contextptr) << "\nranksum "<< w << ", shifted ranksum " << w-M*(M+1.)/2 << "\n"; 3029 double W=M*N+M*(M+1.)/2-w._DOUBLE_val; 3030 if (typetest==0){ 3031 *logptr(contextptr) << "u1=" << W << " ,u2=" << M*N-W ; 3032 if (W>M*(N/2.)) 3033 W=M*N-W; 3034 *logptr(contextptr) << ", u=min(u1,u2)=" << W << "\n"; 3035 } 3036 gen p=wilcoxonp(M,N,contextptr); 3037 if (p.type!=_VECT || p._VECTptr->size()<double(M)*N) return gensizeerr(contextptr); 3038 vecteur & v = *p._VECTptr; 3039 gen P=0,Q=0; 3040 for (double i=0;i<=W;++i){ 3041 P += v[int(i)]; 3042 } 3043 P=P/combMN; 3044 if (typetest==0){ 3045 gen seuil=combMN*alpha/2; 3046 int um(-1); 3047 for (double i=0;i<M*N;++i){ 3048 Q += v[int(i)]; 3049 if (um==-1 && typetest==0 && is_greater(Q,seuil,contextptr)) 3050 um=int(i)-1; 3051 } 3052 *logptr(contextptr) << "Limit value to reject H0 " << um << "\n"; 3053 P=2*P; 3054 } 3055 if (typetest==1) 3056 P=1-P; 3057 *logptr(contextptr) << "P-value " << P << " (" << evalf_double(P,1,contextptr) << "), alpha=" << alpha; 3058 bool ok=is_greater(P,alpha,contextptr); 3059 *logptr(contextptr) << (ok?" H0 not rejected":" H0 rejected") << "\n"; 3060 return ok?1:0; 3061 } 3062 gen w=wilcoxons(*x._VECTptr,m,contextptr); 3063 if (w.type!=_INT_) 3064 return gensizeerr(contextptr); 3065 *logptr(contextptr) << "Wilcoxon 1-sample test, H0 Median=" << m << ", H1 M"; 3066 if (typetest==0) *logptr(contextptr) << "<>"; 3067 if (typetest==1) *logptr(contextptr) << ">"; 3068 if (typetest==-1) *logptr(contextptr) << "<"; 3069 *logptr(contextptr) << m << "\n"; 3070 gen p=wilcoxonp(n); 3071 if (p.type!=_VECT) 3072 return gensizeerr(contextptr); 3073 vecteur & pv=*p._VECTptr; 3074 gen total=0; 3075 unsigned k=0,kmax=-1; 3076 if (typetest==0){ 3077 int wbar=(n*(n+1))/2-w.val; 3078 kmax=giacmin(w.val,wbar); 3079 } 3080 if (typetest==1){ 3081 k=w.val; 3082 kmax=(n*(n+1))/2; 3083 } 3084 if (typetest==-1){ 3085 k=0; 3086 kmax=w.val; 3087 } 3088 for (;k<=kmax;k++){ 3089 total += pv[k]; 3090 } 3091 total=evalf_double(total/pow(plus_two,n,contextptr),1,contextptr); 3092 if (typetest==0) 3093 total=2*total; 3094 *logptr(contextptr) << gettext("Wilcoxon statistic: ") << w << gettext(", p-value: ") << total << gettext(", confidence level: ") << alpha << "\n"; 3095 if (is_greater(total,alpha,contextptr)) 3096 return 1; 3097 return 0; 3098 } 3099 static const char _wilcoxont_s []="wilcoxont"; 3100 static define_unary_function_eval (__wilcoxont,&_wilcoxont,_wilcoxont_s); 3101 define_unary_function_ptr5( at_wilcoxont ,alias_at_wilcoxont,&__wilcoxont,0,true); 3102 giacmin(const std::vector<int> & X)3103 int giacmin(const std::vector<int> & X){ 3104 vector<int>::const_iterator it=X.begin(),itend=X.end(); 3105 int r=RAND_MAX; 3106 for (;it!=itend;++it){ 3107 if (*it<r) 3108 r=*it; 3109 } 3110 return r; 3111 } 3112 giacmax(const std::vector<int> & X)3113 int giacmax(const std::vector<int> & X){ 3114 vector<int>::const_iterator it=X.begin(),itend=X.end(); 3115 int r=-RAND_MAX; 3116 for (;it!=itend;++it){ 3117 if (*it>r) 3118 r=*it; 3119 } 3120 return r; 3121 } 3122 effectif(const std::vector<int> & x,std::vector<int> & eff,int m)3123 void effectif(const std::vector<int> & x,std::vector<int> & eff,int m){ 3124 vector<int>::const_iterator it=x.begin(),itend=x.end(); 3125 for (;it!=itend;++it){ 3126 ++eff[*it-m]; 3127 } 3128 } 3129 somme(const vector<int> & x,const vector<int> & y,vector<int> & z)3130 void somme(const vector<int> & x,const vector<int> &y,vector<int> & z){ 3131 if (&x==&z){ 3132 vector<int>::const_iterator jt=y.begin(),jtend=y.end(); 3133 vector<int>::iterator it=z.begin(),itend=z.end(); 3134 for (;it!=itend&& jt!=jtend;++it,++jt) 3135 *it+=*jt; 3136 for (;jt!=jtend;++jt) 3137 z.push_back(*jt); 3138 return; 3139 } 3140 if (&y==&z){ 3141 somme(y,x,z); 3142 return; 3143 } 3144 z.clear(); 3145 z.reserve(giacmax(int(x.size()),int(y.size()))); 3146 vector<int>::const_iterator it=x.begin(),itend=x.end(),jt=y.begin(),jtend=y.end(); 3147 for (;it!=itend&& jt!=jtend;++it,++jt) 3148 z.push_back(*it+*jt); 3149 for (;it!=itend;++it) 3150 z.push_back(*it); 3151 for (;jt!=jtend;++jt) 3152 z.push_back(*jt); 3153 } 3154 somme(const vector<int> & x)3155 int somme(const vector<int> & x){ 3156 int s=0; 3157 vector<int>::const_iterator it=x.begin(),itend=x.end(); 3158 for (;it!=itend;++it){ 3159 s+=*it; 3160 } 3161 return s; 3162 } 3163 3164 // Chi2 test: valid inputs 3165 // Arg2: A list of probabilities (multinomial adequation test) or a distribution 3166 // or a list of integer values (two samples test) 3167 // Arg1: A list of integer values that is either effectifs or data 3168 // or a list of double values (adequation to density only) 3169 // or a matrix of classes/data (like obtained by classes) _chisquaret(const gen & g,GIAC_CONTEXT)3170 gen _chisquaret(const gen & g,GIAC_CONTEXT){ 3171 if ( g.type==_STRNG && g.subtype==-1) return g; 3172 if (g.type!=_VECT || g._VECTptr->size()<2) 3173 return gensizeerr(contextptr); 3174 vecteur v = *g._VECTptr; 3175 if (g.subtype!=_SEQ__VECT && is_integer_vecteur(v)){ 3176 vector<int> X=vecteur_2_vector_int(v); 3177 int m=giacmin(X),M=giacmax(X); 3178 // guess if data=effectifs or data=values 3179 if (X.size()>=50 && int(X.size())>5*(M-m)){ 3180 *logptr(contextptr) << gettext("Guessing data is a list of values, adequation to uniform discret distribution between ")<<m << gettext(" and ") << M << "\n"; 3181 return _chisquaret(makesequence(g,vecteur(M-m+1,1./(M-m))),contextptr); 3182 } 3183 *logptr(contextptr) << gettext("Guessing data is the list of number of elements in each class, adequation to uniform distribution")<<"\n"; 3184 return _chisquaret(makesequence(g,vecteur(X.size(),1./X.size())),contextptr); 3185 } 3186 // parse arguments for keyword classes 3187 double classmin=class_minimum,classsize=class_size; 3188 bool classdef=false; 3189 for (unsigned i=0;i<v.size();++i){ 3190 if (v[i]==at_classes){ 3191 if (i+1<v.size()){ 3192 gen tmp=evalf_double(v[i+1],1,contextptr); 3193 if (tmp.type!=_DOUBLE_) 3194 return gensizeerr(contextptr); 3195 classmin=tmp._DOUBLE_val; 3196 if (i+2<v.size()){ 3197 gen tmp=evalf_double(v[i+2],1,contextptr); 3198 if (tmp.type!=_DOUBLE_) 3199 return gensizeerr(contextptr); 3200 classsize=tmp._DOUBLE_val; 3201 classdef=true; 3202 } 3203 } 3204 v.erase(v.begin()+i,v.end()); 3205 break; 3206 } 3207 } 3208 gen x=v[0],y=v[1]; 3209 if (x.type!=_VECT) 3210 return gensizeerr(contextptr); 3211 gen ytotal=_plus(y,contextptr); 3212 bool yproba=is_zero(1-ytotal,contextptr); 3213 if (!yproba && y.type==_VECT){ 3214 // >= 2 samples, same law? 3215 // x and y are either classes/effectifs or effectifs (integer values) 3216 if (is_integer_vecteur(*x._VECTptr)){ 3217 if (!is_integer_vecteur(*y._VECTptr)) 3218 return gensizeerr(contextptr); 3219 vector< vector<int> > M(2); 3220 M[0]=vecteur_2_vector_int(*x._VECTptr); 3221 M[1]=vecteur_2_vector_int(*y._VECTptr); 3222 unsigned s=unsigned(v.size()); 3223 for (int j=2;j<int(s);++j){ 3224 if (v[j].type!=_VECT || !is_integer_vecteur(*v[j]._VECTptr)) 3225 return gensizeerr(contextptr); 3226 M.push_back(vecteur_2_vector_int(*v[j]._VECTptr)); 3227 } 3228 vector< vector<int> > effX(s); 3229 int k=int(M[0].size()); 3230 if (k!=int(M[1].size())){ 3231 // build effectifs for x and y 3232 int m=giacmin(M[0]),ma=giacmax(M[0]); 3233 for (unsigned j=1;j<s;++j){ 3234 m=giacmin(m,giacmin(M[j])); 3235 ma=giacmax(ma,giacmax(M[j])); 3236 } 3237 k=ma-m+1; 3238 for (unsigned j=0;j<s;++j){ 3239 effX[j]=vector<int>(k); 3240 effectif(M[j],effX[j],m); 3241 } 3242 } 3243 else 3244 effX=M; 3245 vector<int> eff(k); 3246 for (unsigned j=0;j<s;++j){ 3247 somme(eff,effX[j],eff); 3248 } 3249 double N=somme(eff); 3250 vector<double> NX(s); 3251 for (unsigned j=0;j<s;++j){ 3252 NX[j]=somme(effX[j]); 3253 } 3254 double res=0; 3255 for (unsigned i=0;i<eff.size();++i){ 3256 // theoric proba of class i is eff[i]/N 3257 double p_i=eff[i]/N; 3258 for (unsigned j=0;j<s;++j){ 3259 double tmp1=NX[j]*p_i; 3260 double tmp2=effX[j][i]-tmp1; 3261 res += tmp2*tmp2/tmp1; 3262 } 3263 } 3264 *logptr(contextptr) << s << gettext(" samples Chi2 test result ") << res << gettext(",\nreject adequation if superior to chisquare_icdf(") << k-1 << ",0.95)=" << chisquare_icdf(k-1,0.95,contextptr) << " or chisquare_icdf(" << k-1 <<",1-alpha) if alpha!=5%" << "\n"; 3265 return res; 3266 } 3267 } 3268 if (yproba && x.type==_VECT && !x._VECTptr->empty() && is_integer_vecteur(*x._VECTptr)){ 3269 // if x is a vector of N integers in 1..J and y a vector of J probabilities 3270 // multinomial adequation test: returns sum( (effectif[x==j]-N*y[j])^2/(N*y[j]),j=1..J) 3271 // check y 3272 if (y.type!=_VECT || y._VECTptr->size()<2) 3273 return gensizeerr(contextptr); 3274 int J=int(y._VECTptr->size()); 3275 vector<int> X=vecteur_2_vector_int(*x._VECTptr); 3276 int N=int(X.size()); 3277 gen res=0; 3278 if (N==J){ // X is directly the effectifs 3279 N=0; 3280 for (int i=0;i<J;++i){ 3281 if (X[i]<0) 3282 return gensizeerr(contextptr); 3283 N+=X[i]; 3284 } 3285 for (int j=0;j<J;++j){ 3286 gen tmp= N*y[j]; 3287 res += (X[j]-tmp)*(X[j]-tmp)/tmp; 3288 } 3289 } 3290 else { 3291 vector<int> eff(J+1); 3292 int shift=1; 3293 if (equalposcomp(X,0)) 3294 shift=0; 3295 for (int i=0;i<N;++i){ 3296 if (X[i]<shift || X[i]>=J+shift) 3297 return gensizeerr(gettext("Data should be a class number between 0 and ")+print_INT_(J)); 3298 ++eff[X[i]]; 3299 } 3300 for (int j=shift;j<J+shift;++j){ 3301 gen tmp= N*y[j-shift]; 3302 res += (eff[j]-tmp)*(eff[j]-tmp)/tmp; 3303 } 3304 } 3305 *logptr(contextptr) << gettext("Sample adequation to a finite discrete probability distribution\nChi2 test result ") << res << gettext(",\nreject adequation if superior to chisquare_icdf(") << J-1 << ",0.95)=" << chisquare_icdf(J-1,0.95,contextptr) << " or chisquare_icdf(" << J-1 <<",1-alpha) if alpha!=5%" << "\n"; 3306 return res; 3307 } 3308 int nd; 3309 if ((nd=is_distribution(y))){ 3310 // adequation to a distribution, parameters are estimated or given, depends on the 3311 // number of arguments 3312 gen xorig=x,moyenne=undef,ecart; 3313 vector<int> X,eff; 3314 int minX,maxX,Xclasses; 3315 double efftotal; 3316 if (is_integer_vecteur(*xorig._VECTptr)){ 3317 X=vecteur_2_vector_int(*xorig._VECTptr); 3318 minX=giacmin(X),maxX=giacmax(X); 3319 Xclasses=maxX-minX+1; 3320 if (X.size()>=50 && int(X.size())>=5*Xclasses){ 3321 eff.resize(Xclasses); 3322 effectif(X,eff,minX); 3323 efftotal=somme(eff); 3324 } 3325 else { 3326 minX=0; 3327 eff=X; 3328 efftotal=somme(eff); 3329 Xclasses=int(X.size()); 3330 gen m1=0,m2=0; 3331 for (unsigned i=1;i<X.size();++i){ 3332 m1 += double(i)*X[i]; 3333 m2 += (i*double(i))*X[i]; 3334 } 3335 moyenne=m1/efftotal; 3336 ecart=std::sqrt(efftotal/(efftotal-1))*sqrt(m2/efftotal-moyenne*moyenne,contextptr); 3337 } 3338 } 3339 if (is_undef(moyenne)){ 3340 x=evalf_double(x,1,contextptr); 3341 moyenne=_mean(x,contextptr); 3342 ecart=_stdDev(x,contextptr); 3343 } 3344 if (x.type!=_VECT) 3345 return gensizeerr(contextptr); 3346 int nargs=distrib_nargs(nd); 3347 vecteur w(v.begin()+2,v.end()); 3348 if (y.type==_SYMB) 3349 w=gen2vecteur(y._SYMBptr->feuille); 3350 int s=int(w.size()); 3351 if (s>nargs || s<nargs-2) 3352 return gendimerr(contextptr); 3353 int dof = nargs-s; // adjust number of degree of freedom 3354 if (s==nargs-2){ // 2 params estimation 3355 switch (nd){ 3356 case 1: 3357 w.push_back(moyenne); 3358 w.push_back(ecart); 3359 *logptr(contextptr) << gettext("Normal density, estimating mean and stddev from data ") << moyenne << " " << ecart << "\n"; 3360 break; 3361 case 2:{ 3362 // moyenne=n*p, ecart^2=n*p*(1-p) 3363 gen p=1-ecart*ecart/moyenne; 3364 if (is_greater(0,p,contextptr) || is_greater(p,1,contextptr)) 3365 return gensizeerr(contextptr); 3366 gen n=_round(moyenne/p,contextptr); 3367 p=moyenne/n; 3368 if (is_greater(0,p,contextptr) || is_greater(p,1,contextptr)) 3369 return gensizeerr(contextptr); 3370 *logptr(contextptr) << gettext("Binomial: estimating n and p from data ") << n << " " << p << "\n"; 3371 w.push_back(n); 3372 w.push_back(p); 3373 break; 3374 } 3375 case 3:{ 3376 gen p=moyenne/(ecart*ecart); 3377 if (is_greater(0,p,contextptr) || is_greater(p,1,contextptr)) 3378 return gensizeerr(contextptr); 3379 gen n=_round(moyenne*p/(1-p),contextptr); 3380 *logptr(contextptr) << gettext("Negative binomial: estimating n and p from data ") << n << " " << p << "\n"; 3381 w.push_back(n); 3382 w.push_back(p); 3383 break; 3384 } 3385 case 8: 3386 // weibull: 2 parameters to estimate k, lambda, theta assumed to be 0 3387 // lambda*Gamma(1+1/k)=moyenne, lambda^2*Gamma(1+2/k)=ecart^2-moyenne^2 3388 return gensizeerr(contextptr); 3389 break; 3390 case 13:{ // moyenne +/- sqrt(3)*ecart 3391 gen a=moyenne-std::sqrt(3.0)*ecart; 3392 gen b=moyenne+std::sqrt(3.0)*ecart; 3393 w.push_back(a); 3394 w.push_back(b); 3395 *logptr(contextptr) << gettext("Uniform density, estimating interval boundaries from data ") << a << " " << b << "\n"; 3396 break; 3397 } 3398 default: 3399 return gendimerr(contextptr); 3400 } // end switch 3401 s=int(w.size()); 3402 } // end 2 paramters to estimate 3403 if (s==nargs-1){ // 1 params estimation 3404 switch (nd){ 3405 case 1: {// normald with 1 parameter known, w[0] default is stddev 3406 if (w.empty()) 3407 return gendimerr(contextptr); 3408 if (w[0].is_symb_of_sommet(at_equal)){ 3409 if (w[0][1]==at_mean){ 3410 *logptr(contextptr) << gettext("Normal density, estimating std deviation from data ") << ecart << "\n"; 3411 w.push_back(ecart); 3412 } 3413 w[0]=w[0][2]; 3414 } 3415 if (w.size()==1){ 3416 *logptr(contextptr) << gettext("Normal density, estimating mean from data ") << moyenne << "\n"; 3417 w.insert(w.begin(),moyenne); 3418 } 3419 break; 3420 } 3421 case 4: 3422 w.push_back(moyenne); 3423 *logptr(contextptr) << gettext("Poisson distribution, estimating parameter from data mean ") << moyenne << "\n"; 3424 break; 3425 default: 3426 return gensizeerr(contextptr); 3427 } 3428 s=int(w.size()); 3429 } 3430 w.push_back(0); // last argument of the distribution 3431 gen loi; 3432 bool discret=is_discrete_distribution(nd); 3433 if (discret) 3434 loi=symbolic(y.type==_FUNC?*y._FUNCptr:y._SYMBptr->sommet,gen(w,_SEQ__VECT)); 3435 else { 3436 loi=cdf(nd); 3437 if (loi.type!=_FUNC) 3438 return gensizeerr(contextptr); 3439 loi=symbolic(*loi._FUNCptr,gen(w,_SEQ__VECT)); 3440 } 3441 // now we have the law, do the test! 3442 // if classdef=false, for discrete distributions use classsize=1, classmin=0 3443 // for densities guess from data? 3444 if (is_integer_vecteur(*xorig._VECTptr)){ 3445 if (classdef){ 3446 // not yet supported... 3447 *logptr(contextptr) << "Warning, using default class minimum=0 and class_size=1 for discrete distribution" << "\n"; 3448 } 3449 dof=Xclasses-1-dof; 3450 gen res=0; 3451 for (unsigned i=0;i<eff.size();++i){ 3452 // theoric proba of class i from law 3453 loi._SYMBptr->feuille._VECTptr->back()=minX+int(i)+(discret?0:-.5); 3454 gen tmp=evalf_double(loi,1,contextptr); 3455 if (!discret){ // could be improved, we compute twice the cdf 3456 loi._SYMBptr->feuille._VECTptr->back()=minX+int(i+1)-.5; 3457 tmp=evalf_double(loi,1,contextptr)-tmp; 3458 } 3459 double p_i=tmp._DOUBLE_val; 3460 double tmp1=efftotal*p_i; 3461 double tmp2=eff[i]-tmp1; 3462 res += tmp2*tmp2/tmp1; 3463 } 3464 loi._SYMBptr->feuille._VECTptr->back()=identificateur("."); 3465 *logptr(contextptr) << gettext("Sample adequation to ")<< loi <<gettext(", Chi2 test result ") << res << gettext(",\nreject adequation if superior to chisquare_icdf(") << dof << ",0.95)=" << chisquare_icdf(dof,0.95,contextptr) << gettext(" or chisquare_icdf(") << dof <<gettext(",1-alpha) if alpha!=5%") << "\n"; 3466 return res; 3467 } // end if is_integer_vecteur(x) 3468 if (discret) 3469 return gensizeerr(gettext("Non integer values for adequation to a discrete distribution")); 3470 // compute classes and recall itself 3471 matrice m; 3472 if (ckmatrix(x)) 3473 m=*x._VECTptr; 3474 else 3475 m=effectifs(*x._VECTptr,classmin,classsize,contextptr); 3476 dof=int(m.size())-dof; 3477 if (dof<2) 3478 return gensizeerr(gettext("Not enough degree of freedom with default values. Add classes,classmin,classsize parameters")); 3479 efftotal=0.0; 3480 for (unsigned i=0;i<m.size();++i){ 3481 gen effi=m[i][1]; 3482 if (effi.type!=_INT_) 3483 return gensizeerr(contextptr); 3484 efftotal += effi.val; 3485 } 3486 gen res=0; 3487 for (unsigned i=0;i<m.size();++i){ 3488 gen tmp=m[i][0]; 3489 if (!tmp.is_symb_of_sommet(at_interval)) 3490 return gensizeerr(contextptr); 3491 gen tmp0=evalf_double(tmp._SYMBptr->feuille[0],1,contextptr), 3492 tmp1=evalf_double(tmp._SYMBptr->feuille[1],1,contextptr), 3493 tmp2=m[i][1]; 3494 if (tmp0.type!=_DOUBLE_ || tmp1.type!=_DOUBLE_ || tmp1._DOUBLE_val<=tmp0._DOUBLE_val) 3495 return gensizeerr(contextptr); 3496 loi._SYMBptr->feuille._VECTptr->back()=tmp1; 3497 gen theoric=evalf_double(loi,1,contextptr); 3498 loi._SYMBptr->feuille._VECTptr->back()=tmp0; 3499 theoric=theoric-evalf_double(loi,1,contextptr); 3500 theoric=efftotal*theoric; 3501 res += (tmp2-theoric)*(tmp2-theoric)/theoric; 3502 } 3503 loi._SYMBptr->feuille._VECTptr->back()=identificateur("."); 3504 *logptr(contextptr) << gettext("Sample adequation to ") << loi <<gettext(", Chi2 test result ") << res << gettext(",\nreject adequation if superior to chisquare_icdf(") << dof << ",0.95)=" << chisquare_icdf(dof,0.95,contextptr) << " or chisquare_icdf(" << dof <<",1-alpha) if alpha!=5%" << "\n"; 3505 return res; 3506 } 3507 return undef; 3508 } 3509 static const char _chisquaret_s []="chisquaret"; 3510 static define_unary_function_eval (__chisquaret,&_chisquaret,_chisquaret_s); 3511 define_unary_function_ptr5( at_chisquaret ,alias_at_chisquaret,&__chisquaret,0,true); 3512 3513 // normalt arguments: 3514 // arg1 = 1/ [x,n] number of success, number of trials 3515 // = 2/ [mu1,n] mean, sample size 3516 // = 3/ data (list of values) 3517 // arg2 = 1/ p proportion 3518 // = 2, 3/ mu population mean, or data 3519 // arg3 optionnal for 2, 3/: sigma. If not given and 3520 // arg1=[int,int] and arg2=p in ]0,1[, sigma is derived from p 3521 // nextarg < > or != alternative hypothesis 3522 // nextarg optional alpha level of confidence, default value 0.05 zttest(const gen & g,bool ztest,GIAC_CONTEXT)3523 gen zttest(const gen & g,bool ztest,GIAC_CONTEXT){ 3524 if ( g.type==_STRNG && g.subtype==-1) return g; 3525 if (g.type!=_VECT || g._VECTptr->size()<3 || g._VECTptr->front().type!=_VECT) 3526 return gensizeerr(contextptr); 3527 vecteur v =*g._VECTptr; 3528 vecteur v0=*v[0]._VECTptr; 3529 gen MU1=evalf_double(v[1],1,contextptr),test=undef; 3530 double mu0,mu1,alpha=0.05,sigma=1; 3531 double n0,n1=0; // sample size for data0 and data1, 0 for data1 if not a sample 3532 int dof=0; 3533 if (v0.size()<2) 3534 return gensizeerr(contextptr); 3535 bool proportion=false,data0=false,data1=MU1.type==_VECT; 3536 if (data1){ 3537 n1=double(MU1._VECTptr->size()); 3538 dof = int(n1)-1; // mean estimated from data 3539 gen tmp=_mean(MU1,contextptr); 3540 if (tmp.type!=_DOUBLE_) 3541 return gensizeerr(contextptr); 3542 mu1=tmp._DOUBLE_val; 3543 *logptr(contextptr) << gettext("Estimated mean from sample(s)") << mu1 << "\n"; 3544 } 3545 else { 3546 if (MU1.type!=_DOUBLE_) 3547 return gensizeerr(contextptr); 3548 mu1=MU1._DOUBLE_val; 3549 } 3550 if ( (n0=double(v0.size()))==2 && is_integer(v0[1])){ 3551 gen v00=evalf_double(v0[0],1,contextptr); 3552 if (v00.type!=_DOUBLE_) 3553 return gensizeerr(contextptr); 3554 // arg1 =[mu1,n] or [x,n] 3555 n0=evalf_double(v0[1],1,contextptr)._DOUBLE_val; 3556 dof += int(n0); 3557 proportion=ztest && is_integer(v0[0]) && MU1.type==_DOUBLE_ && MU1._DOUBLE_val>0 && MU1._DOUBLE_val<1 && v[2].type==_FUNC; 3558 if (proportion) 3559 mu0=v00._DOUBLE_val/n0; 3560 else 3561 mu0=v00._DOUBLE_val; 3562 } 3563 else { 3564 dof += int(n0); 3565 data0=true; 3566 gen tmp=evalf_double(_mean(v0,contextptr),1,contextptr); 3567 if (tmp.type!=_DOUBLE_) 3568 return gensizeerr(contextptr); 3569 mu0=tmp._DOUBLE_val; 3570 } 3571 if (v[2].type==_FUNC){ 3572 // estimate sigma value from data 3573 --dof; 3574 gen S=0; 3575 if (data0) 3576 S=n0*_variance(v[0],contextptr); 3577 if (data1) 3578 S += n1*_variance(v[1],contextptr); 3579 S=evalf_double(S/dof,1,contextptr); 3580 if (S.type!=_DOUBLE_ || S._DOUBLE_val<=0) 3581 return gensizeerr(gettext("Unable to guess sigma using data")); 3582 sigma=std::sqrt(S._DOUBLE_val); 3583 *logptr(contextptr) << gettext("Estimated sigma from sample(s)") << sigma << "\n"; 3584 if (v.size()>4) 3585 return gendimerr(contextptr); 3586 if (v.size()==4){ 3587 gen tmp=evalf_double(v[3],1,contextptr); 3588 if (tmp.type!=_DOUBLE_ || tmp._DOUBLE_val<=0 || tmp._DOUBLE_val>=1) 3589 return gensizeerr(contextptr); 3590 alpha=tmp._DOUBLE_val; 3591 } 3592 test=v[2]; 3593 } 3594 else { 3595 if (v.size()<4) 3596 return gendimerr(contextptr); 3597 gen tmp=evalf_double(v[2],1,contextptr); 3598 if (tmp.type!=_DOUBLE_ || tmp._DOUBLE_val<=0) 3599 return gensizeerr(contextptr); 3600 sigma=tmp._DOUBLE_val; 3601 test=v[3]; 3602 if (v.size()==5){ 3603 tmp=evalf_double(v[4],1,contextptr); 3604 if (tmp.type!=_DOUBLE_ || tmp._DOUBLE_val<=0 || tmp._DOUBLE_val>=1) 3605 return gensizeerr(contextptr); 3606 alpha=tmp._DOUBLE_val; 3607 } 3608 } 3609 if (test==at_inferieur_strict) 3610 test=-1; 3611 if (test==at_superieur_strict) 3612 test=1; 3613 if (test==at_different) 3614 test=0; 3615 if (test!=0 && test!=-1 && test!=1) 3616 return gensizeerr(gettext("Bad alternative hypothesis, should be '>', '<' or '!='")); 3617 // Ready to test: compare mu0 to mu1-f(alpha)*sigma 3618 gen Falpha; 3619 if (test==0) 3620 Falpha=ztest?_normal_icdf(makesequence(0,1,1-alpha/2),contextptr):_student_icdf(makesequence(dof,1-alpha),contextptr); 3621 else 3622 Falpha=ztest?_normal_icdf(makesequence(0,1,1-alpha),contextptr):_student_icdf(makesequence(dof,1-alpha),contextptr); 3623 double falpha=evalf_double(Falpha,1,contextptr)._DOUBLE_val; 3624 bool ok=true; 3625 double sqrtn0=std::sqrt(n0); 3626 if (test==-1){ 3627 if (mu0<mu1-sigma*falpha/sqrtn0) 3628 ok=false; 3629 } 3630 if (test==0){ 3631 // mu0 should be < mu1 for test to fail 3632 if (mu0<mu1-sigma*falpha/sqrtn0 || mu0>mu1+sigma*falpha/sqrtn0) 3633 ok=false; 3634 } 3635 if (test==1){ 3636 if (mu0>mu1+sigma*falpha/sqrtn0) 3637 ok=false; 3638 } 3639 *logptr(contextptr) << gettext("*** TEST RESULT ") << (ok?"1 ***":"0 ***") << "\n" << "Summary " << (ztest?"Z-Test":"T-Test") << " null hypothesis H0 " << (proportion?"p1=p2":"mu1=mu2") << ", alt. hyp. H1 mu1" << (test==-1? "<":(test==0?"!=":">")) << "mu2." << "\n"; 3640 *logptr(contextptr) << gettext("Test returns 0 if probability to observe data is less than ") << alpha << "\n"; 3641 *logptr(contextptr) << gettext("(null hyp. mu1=mu2 rejected with less than alpha probability error)") << "\n"; 3642 *logptr(contextptr) << gettext("Test returns 1 otherwise (can not reject null hypothesis)") << "\n"; 3643 *logptr(contextptr) << gettext("Data mean mu1=") << mu0 << gettext(", population mean mu2=") << mu1 ; 3644 if (!ztest) 3645 *logptr(contextptr) << gettext(", degrees of freedom ") << dof; 3646 *logptr(contextptr) << "\n" << "alpha level " << alpha << ", multiplier*stddev/sqrt(sample size)= " << falpha << "*" << sigma << "/" << sqrtn0 << "\n"; 3647 return (ok?1:0); 3648 } _normalt(const gen & g,GIAC_CONTEXT)3649 gen _normalt(const gen & g,GIAC_CONTEXT){ 3650 return zttest(g,true,contextptr); 3651 } 3652 static const char _normalt_s []="normalt"; 3653 static define_unary_function_eval (__normalt,&_normalt,_normalt_s); 3654 define_unary_function_ptr5( at_normalt ,alias_at_normalt,&__normalt,0,true); 3655 _studentt(const gen & g,GIAC_CONTEXT)3656 gen _studentt(const gen & g,GIAC_CONTEXT){ 3657 return zttest(g,false,contextptr); 3658 } 3659 static const char _studentt_s []="studentt"; 3660 static define_unary_function_eval (__studentt,&_studentt,_studentt_s); 3661 define_unary_function_ptr5( at_studentt ,alias_at_studentt,&__studentt,0,true); 3662 3663 // return 0 if not distrib 3664 // 1 normal, 2 binomial, 3 negbinomial, 4 poisson, 5 student, 3665 // 6 fisher, 7 cauchy, 8 weibull, 9 betad, 10 gammad, 11 chisquare 3666 // 12 geometric, 13 uniformd, 14 exponentiald is_distribution(const gen & args)3667 int is_distribution(const gen & args){ 3668 if (args.type==_SYMB && args._SYMBptr->sommet!=at_exp){ 3669 int res=is_distribution(args._SYMBptr->sommet) ; 3670 if (!res) 3671 return res; 3672 int s=distrib_nargs(res); 3673 if (s!=int(gen2vecteur(args._SYMBptr->feuille).size())) 3674 return 0; 3675 return res; 3676 } 3677 if (args.type==_FUNC){ 3678 if (args==at_normald || args==at_NORMALD) 3679 return 1; 3680 if (args==at_binomial || args==at_BINOMIAL) 3681 return 2; 3682 if (args==at_negbinomial) 3683 return 3; 3684 if (args==at_poisson || args==at_POISSON) 3685 return 4; 3686 if (args==at_student || args==at_studentd) 3687 return 5; 3688 if (args==at_fisher || args==at_fisherd || args==at_snedecor) 3689 return 6; 3690 if (args==at_cauchy || args==at_cauchyd) 3691 return 7; 3692 if (args==at_weibull || args==at_weibulld) 3693 return 8; 3694 if (args==at_betad) 3695 return 9; 3696 if (args==at_gammad) 3697 return 10; 3698 if (args==at_chisquared || args==at_chisquare) 3699 return 11; 3700 if (args==at_geometric) 3701 return 12; 3702 if (args==at_uniform || args==at_uniformd) 3703 return 13; 3704 if (args==at_exp || args==at_exponential || args==at_exponentiald) 3705 return 14; 3706 } 3707 return 0; 3708 } 3709 distribution(int nd)3710 gen distribution(int nd){ 3711 static vecteur * d_static=0; 3712 if (!d_static) 3713 d_static=new vecteur(makevecteur(at_normald,at_binomial,at_negbinomial,at_poisson,at_studentd,at_fisherd,at_cauchyd,at_weibulld,at_betad,at_gammad,at_chisquared,at_geometric,at_uniformd,at_exponentiald)); 3714 if (nd<=0 || nd>int(d_static->size())) 3715 return undef; 3716 return (*d_static)[nd-1]; 3717 } 3718 distrib_nargs(int nd)3719 int distrib_nargs(int nd){ 3720 switch (nd){ 3721 case 4: case 5: case 11: case 12: case 14: 3722 return 1; 3723 #if 0 // Luka Marohnić comment: Furthermore, I propose excluding case 8 in distrib_nargs routine in moyal.cc. Namely, weibulld(k,lambda,theta) is evaluated to a real number, the value of Weibull PDF with k=a and lambda=b at point x=theta (not as a three-parameter Weibull). In contrast, weibulld(k,lambda) is a function and should be detected by is_distribution, but fails on the mentioned case 8. Perhaps distrib_nargs should return 2 on weibull distribution 3724 case 8: 3725 return 3; 3726 #endif 3727 default: 3728 return 2; 3729 } 3730 } 3731 _mgf(const gen & g,GIAC_CONTEXT)3732 gen _mgf(const gen & g,GIAC_CONTEXT){ 3733 if ( g.type==_STRNG && g.subtype==-1) return g; 3734 if (g.type==_SYMB){ 3735 vecteur v(gen2vecteur(g._SYMBptr->feuille)); 3736 v.insert(v.begin(),g._SYMBptr->sommet); 3737 return _mgf(v,contextptr); 3738 } 3739 int nd; 3740 if (g.type!=_VECT || g._VECTptr->empty() || !(nd=is_distribution(g._VECTptr->front())) ) 3741 return gensizeerr(contextptr); 3742 vecteur & v=*g._VECTptr; 3743 int s=int(v.size()); 3744 if (s!=distrib_nargs(nd)+1) 3745 return gensizeerr(contextptr); 3746 gen t(identificateur("t")); 3747 if (nd==1) 3748 return exp(v[1]*t+plus_one_half*pow(v[2],2,contextptr)*pow(t,2,contextptr),contextptr); 3749 if (nd==2) 3750 return pow((1-v[2])+v[2]*exp(t,contextptr),v[1],contextptr); 3751 if (nd==3) 3752 return pow(v[2]/(1-(1-v[2])*exp(t,contextptr)),v[1],contextptr); 3753 if (nd==4) 3754 return exp(v[1]*(exp(t,contextptr)-1),contextptr); 3755 if (nd==10) 3756 return pow(1-t/v[2],-v[1],contextptr); 3757 if (nd==11) 3758 return pow(1-2*t,-v[1]/2,contextptr); 3759 if (nd==12) 3760 return v[1]*exp(t,contextptr)/(1-(1-v[1])*exp(t,contextptr)); 3761 if (nd==13) 3762 return (exp(t*v[2],contextptr)-exp(t*v[1],contextptr))/(t*(v[2]-v[1])); 3763 if (nd==14) 3764 return inv(1-t/v[1],contextptr); 3765 return undef; 3766 } 3767 static const char _mgf_s []="mgf"; 3768 static define_unary_function_eval (__mgf,&_mgf,_mgf_s); 3769 define_unary_function_ptr5( at_mgf ,alias_at_mgf,&__mgf,0,true); 3770 icdf(int n)3771 gen icdf(int n){ 3772 static vecteur * icdf_static=0; 3773 if (!icdf_static) 3774 icdf_static=new vecteur(makevecteur(at_normald_icdf,at_binomial_icdf,undef,at_poisson_icdf,at_studentd_icdf,at_fisherd_icdf,at_cauchyd_icdf,at_weibulld_icdf,at_betad_icdf,at_gammad_icdf,at_chisquared_icdf,at_geometric_icdf,at_uniformd_icdf,at_exponentiald_icdf)); 3775 if (n<=0 || n>int(icdf_static->size())) 3776 return undef; 3777 return (*icdf_static)[n-1]; 3778 } 3779 cdf(int n)3780 gen cdf(int n){ 3781 static vecteur * cdf_static=0; 3782 if (!cdf_static) 3783 cdf_static=new vecteur(makevecteur(at_normald_cdf,at_binomial_cdf,undef,at_poisson_cdf,at_studentd_cdf,at_fisherd_cdf,at_cauchyd_cdf,at_weibulld_cdf,at_betad_cdf,at_gammad_cdf,at_chisquared_cdf,at_geometric_cdf,at_uniformd_cdf,at_exponentiald_cdf)); 3784 if (n<=0 || n>int(cdf_static->size())) 3785 return undef; 3786 return (*cdf_static)[n-1]; 3787 } 3788 3789 // set a and b to the boundaries of the support of distrib number nd 3790 // truncate infinities to 100 if truncate is true distrib_support(int nd,gen & a,gen & b,bool truncate)3791 bool distrib_support(int nd,gen & a,gen &b,bool truncate){ 3792 a=truncate?gnuplot_xmin:minus_inf; 3793 b=truncate?gnuplot_xmax:plus_inf; 3794 if (nd==2 || nd==3 || nd==4 || nd==9 || nd==10 || nd==11 || nd==14) 3795 a=0; 3796 if (nd==9) 3797 b=1; 3798 if (nd==12){ 3799 a=1; 3800 b=10; 3801 } 3802 if (nd==2 || nd==3 || nd==4 || nd==12) 3803 return false; 3804 return true; 3805 } 3806 _cdf(const gen & g,GIAC_CONTEXT)3807 gen _cdf(const gen & g,GIAC_CONTEXT){ 3808 if ( g.type==_STRNG && g.subtype==-1) return g; 3809 int nd; 3810 if (g.type!=_VECT || g._VECTptr->empty()) 3811 return gensizeerr(contextptr); 3812 vecteur w=*g._VECTptr; 3813 int s=int(w.size()); 3814 if (!(nd=is_distribution(w.front()))){ 3815 if (g.subtype!=_SEQ__VECT){ 3816 w=gen2vecteur(_sort(g,contextptr)); 3817 s=int(w.size()); 3818 vecteur res; 3819 for (int i=0;i<s;++i){ 3820 res.push_back(makevecteur(w[i],gen(i)/gen(s))); 3821 } 3822 return res; 3823 } 3824 vector<giac_double> D; 3825 gen w0=evalf_double(w[0],1,contextptr); 3826 if (s==2 && ckmatrix(w0)){ // list of classes/effectifs 3827 matrice m=*w0._VECTptr; 3828 if (m.empty()) return gensizeerr(contextptr); 3829 if (m.front()._VECTptr->size()!=2) 3830 m=mtran(m); 3831 if (m.front()._VECTptr->size()!=2) 3832 return gensizeerr(contextptr); 3833 gen_sort_f(m.begin(),m.end(),first_ascend_sort); 3834 s=int(m.size()); 3835 double tot=0,cur=0; 3836 for (int i=0;i<s;++i){ 3837 if (m[i][1].type!=_DOUBLE_) 3838 return gensizeerr(contextptr); 3839 tot += m[i][1]._DOUBLE_val; 3840 } 3841 if (w[1]!=at_plot){ 3842 w[1]=evalf_double(w[1],1,contextptr); 3843 if (w[1].type!=_DOUBLE_ || w[1]._DOUBLE_val<=0 || w[1]._DOUBLE_val>=1) 3844 return gensizeerr(contextptr); 3845 tot *= w[1]._DOUBLE_val; 3846 for (int i=0;i<s;++i){ 3847 cur += m[i][1]._DOUBLE_val; 3848 if (cur>=tot) 3849 return m[i][0]; 3850 } 3851 } 3852 // plot cdf for frequencies 3853 vecteur res; res.reserve(2*s); 3854 for (int i=0;i<s;++i){ 3855 gen mi0=m[i][0]; 3856 if (mi0.is_symb_of_sommet(at_interval) && mi0._SYMBptr->feuille.type==_VECT && mi0._SYMBptr->feuille._VECTptr->size()==2) 3857 mi0=(mi0._SYMBptr->feuille._VECTptr->front()+mi0._SYMBptr->feuille._VECTptr->back())/2; 3858 res.push_back(mi0+(cur/tot)*cst_i); 3859 cur += m[i][1]._DOUBLE_val; 3860 res.push_back(mi0+(cur/tot)*cst_i); 3861 } 3862 return _polygone_ouvert(gen(res,_SEQ__VECT),contextptr); 3863 } 3864 if (s!=2 || w0.type!=_VECT || !convert(*w0._VECTptr,D,false)) 3865 return gensizeerr(contextptr); 3866 sort(D.begin(),D.end()); 3867 s=int(D.size()); 3868 if (w[1]!=at_plot){ 3869 gen tmp=evalf_double(w[1],1,contextptr); 3870 if (tmp.type!=_DOUBLE_) 3871 return gensizeerr(contextptr); 3872 return gen(dichotomy(D,tmp._DOUBLE_val))/gen(s); 3873 } 3874 vecteur res; 3875 for (int i=0;i<s;++i){ 3876 res.push_back(double(D[i])+i/double(s)*cst_i); 3877 res.push_back(double(D[i])+(i+1)/double(s)*cst_i); 3878 } 3879 return _polygone_ouvert(gen(res,_SEQ__VECT),contextptr); 3880 } 3881 if (s && w.front().type==_SYMB){ 3882 vecteur v(gen2vecteur(w.front()._SYMBptr->feuille)); 3883 v.insert(v.begin(),w.front()._SYMBptr->sommet); 3884 for (unsigned j=1;j<w.size();++j) 3885 v.push_back(w[j]); 3886 return _cdf(v,contextptr); 3887 } 3888 bool plot=false; 3889 if (w.back()==at_plot || w.back()==at_plotfunc){ 3890 if (nd==13){ 3891 gen a=w[1],b=w[2]; 3892 return _segment(makesequence(a,b+cst_i),contextptr); 3893 } 3894 w.back()=vx_var; 3895 plot=true; 3896 } 3897 if (s==distrib_nargs(nd)+1){ 3898 w.push_back(vx_var); 3899 ++s; 3900 } 3901 if (s!=distrib_nargs(nd)+2) 3902 return gensizeerr(contextptr); 3903 gen res=gen(vecteur(w.begin()+1,w.end()),_SEQ__VECT); 3904 res=cdf(nd)(res,contextptr); 3905 if (plot){ 3906 gen a,b; 3907 if (distrib_support(nd,a,b,true)) // true if density 3908 return _plotfunc(makesequence(res,symb_equal(vx_var,symb_interval(a,b))),contextptr); 3909 if (nd==2) // binomial 3910 b=w[1]; 3911 if (a.type!=_INT_ || !is_integral(b) || b.type!=_INT_ || b.val<=0) 3912 return gensizeerr(contextptr); 3913 int A=a.val,B=b.val; 3914 vecteur v; 3915 for (int i=A;i<B;++i){ 3916 gen y=subst(res,vx_var,i,false,contextptr); 3917 v.push_back(i+cst_i*y); 3918 v.push_back(i+1+cst_i*y); 3919 } 3920 gen y=subst(res,vx_var,B,false,contextptr); 3921 v.push_back(B+cst_i*y); 3922 return _polygone_ouvert(gen(v,_SEQ__VECT),contextptr); 3923 } 3924 return res; 3925 } 3926 static const char _cdf_s []="cdf"; 3927 static define_unary_function_eval (__cdf,&_cdf,_cdf_s); 3928 define_unary_function_ptr5( at_cdf ,alias_at_cdf,&__cdf,0,true); 3929 _plotcdf(const gen & g,GIAC_CONTEXT)3930 gen _plotcdf(const gen & g,GIAC_CONTEXT){ 3931 if ( g.type==_STRNG && g.subtype==-1) return g; 3932 vecteur v(makevecteur(g,at_plot)); 3933 if (g.type==_VECT && g.subtype==_SEQ__VECT){ 3934 v=*g._VECTptr; 3935 v.push_back(at_plot); 3936 } 3937 return _cdf(gen(v,_SEQ__VECT),contextptr); 3938 } 3939 static const char _plotcdf_s []="plotcdf"; 3940 static define_unary_function_eval (__plotcdf,&_plotcdf,_plotcdf_s); 3941 define_unary_function_ptr5( at_plotcdf ,alias_at_plotcdf,&__plotcdf,0,true); 3942 3943 static const char _cdfplot_s []="cdfplot"; 3944 static define_unary_function_eval (__cdfplot,&_plotcdf,_cdfplot_s); 3945 define_unary_function_ptr5( at_cdfplot ,alias_at_cdfplot,&__cdfplot,0,true); 3946 _icdf(const gen & g,GIAC_CONTEXT)3947 gen _icdf(const gen & g,GIAC_CONTEXT){ 3948 if ( g.type==_STRNG && g.subtype==-1) return g; 3949 int nd; 3950 if (g.type!=_VECT || g._VECTptr->empty()) 3951 return gensizeerr(contextptr); 3952 vecteur w=*g._VECTptr; 3953 int s=int(w.size()); 3954 if (!(nd=is_distribution(w.front()))){ 3955 if (g.subtype!=_SEQ__VECT || s!=2) 3956 return gensizeerr(contextptr); 3957 if (w[1]==at_plot) 3958 return _symetrie(makesequence(_droite(makesequence(0,1+cst_i),contextptr),_cdf(g,contextptr)),contextptr); 3959 return _quantile(g,contextptr); 3960 } 3961 if (s && w.front().type==_SYMB){ 3962 vecteur v(gen2vecteur(w.front()._SYMBptr->feuille)); 3963 v.insert(v.begin(),w.front()._SYMBptr->sommet); 3964 for (unsigned j=1;j<w.size();++j) 3965 v.push_back(w[j]); 3966 return _icdf(v,contextptr); 3967 } 3968 bool plot=false; 3969 if (w.back()==at_plot || w.back()==at_plotfunc){ 3970 if (nd==13){ 3971 gen a=w[1],b=w[2]; 3972 return _segment(makesequence(cst_i*a,1+cst_i*b),contextptr); 3973 } 3974 w.back()=vx_var; 3975 plot=true; 3976 } 3977 if (s==distrib_nargs(nd)+1){ 3978 w.push_back(vx_var); 3979 ++s; 3980 } 3981 if (s!=distrib_nargs(nd)+2) 3982 return gensizeerr(contextptr); 3983 gen res=gen(vecteur(w.begin()+1,w.end()),_SEQ__VECT); 3984 res=icdf(nd)(res,contextptr); 3985 if (plot){ 3986 gen a,b; 3987 if (distrib_support(nd,a,b,true)) // true if density 3988 return _plotfunc(makesequence(res,symb_equal(vx_var,symb_interval(0,1))),contextptr); 3989 return _symetrie(makesequence(_droite(makesequence(0,1+cst_i),contextptr),_cdf(g,contextptr)),contextptr); 3990 } 3991 return res; 3992 } 3993 static const char _icdf_s []="icdf"; 3994 static define_unary_function_eval (__icdf,&_icdf,_icdf_s); 3995 define_unary_function_ptr5( at_icdf ,alias_at_icdf,&__icdf,0,true); 3996 3997 3998 // kind=0: BesselI, =1 BesselJ, =2 BesselK, =3 BesselY Bessel(const gen & g,int kind,GIAC_CONTEXT)3999 gen Bessel(const gen & g,int kind,GIAC_CONTEXT){ 4000 #if defined BESTA_OS || defined FXCG || defined KHICAS 4001 return gensizeerr(gettext("Bessel not implemented")); 4002 #else 4003 int n; 4004 gen a,x; 4005 if (!find_n_x(g,n,x,a)) 4006 return gensizeerr(contextptr); 4007 if (has_evalf(x,a,1,contextptr) && a.type==_DOUBLE_ 4008 #ifndef POCKETCAS 4009 && j0!=NULL 4010 #endif 4011 ){ 4012 double d=a._DOUBLE_val; 4013 switch (kind){ 4014 case 1: 4015 if (n==0) return j0(d); 4016 if (n==1) return j1(d); 4017 return jn(n,d); 4018 case 3: 4019 if (n==0) return y0(d); 4020 if (n==1) return y1(d); 4021 return yn(n,d); 4022 } 4023 } 4024 gen gn=gen(makevecteur(n,x),_SEQ__VECT); 4025 switch (kind){ 4026 case 0: 4027 return symbolic(at_BesselI,gn); 4028 case 1: 4029 return symbolic(at_BesselJ,gn); 4030 case 2: 4031 return symbolic(at_BesselK,gn); 4032 case 3: 4033 return symbolic(at_BesselY,gn); 4034 } 4035 return gensizeerr(gettext("Bessel")); 4036 #endif 4037 } bessel(const gen & g,int kind,GIAC_CONTEXT)4038 gen bessel(const gen & g,int kind,GIAC_CONTEXT){ 4039 if (g.type==_VECT && g._VECTptr->size()>=2) 4040 return Bessel(makesequence(g[1],g[0]),kind,contextptr); 4041 return gensizeerr(contextptr); 4042 } 4043 4044 // numerical eval BesselJ(n,x) and BesselY(n,x) are implemented for x.type==_DOUBLE_ _BesselI(const gen & g,GIAC_CONTEXT)4045 gen _BesselI(const gen & g,GIAC_CONTEXT){ 4046 if ( g.type==_STRNG && g.subtype==-1) return g; 4047 return Bessel(g,0,contextptr); 4048 } 4049 static const char _BesselI_s []="BesselI"; 4050 static define_unary_function_eval (__BesselI,&_BesselI,_BesselI_s); 4051 define_unary_function_ptr5( at_BesselI ,alias_at_BesselI,&__BesselI,0,true); 4052 _BesselJ(const gen & g,GIAC_CONTEXT)4053 gen _BesselJ(const gen & g,GIAC_CONTEXT){ 4054 if ( g.type==_STRNG && g.subtype==-1) return g; 4055 return Bessel(g,1,contextptr); 4056 } 4057 static const char _BesselJ_s []="BesselJ"; 4058 static define_unary_function_eval (__BesselJ,&_BesselJ,_BesselJ_s); 4059 define_unary_function_ptr5( at_BesselJ ,alias_at_BesselJ,&__BesselJ,0,true); 4060 _BesselK(const gen & g,GIAC_CONTEXT)4061 gen _BesselK(const gen & g,GIAC_CONTEXT){ 4062 if ( g.type==_STRNG && g.subtype==-1) return g; 4063 return Bessel(g,2,contextptr); 4064 } 4065 static const char _BesselK_s []="BesselK"; 4066 static define_unary_function_eval (__BesselK,&_BesselK,_BesselK_s); 4067 define_unary_function_ptr5( at_BesselK ,alias_at_BesselK,&__BesselK,0,true); 4068 _BesselY(const gen & g,GIAC_CONTEXT)4069 gen _BesselY(const gen & g,GIAC_CONTEXT){ 4070 if ( g.type==_STRNG && g.subtype==-1) return g; 4071 return Bessel(g,3,contextptr); 4072 } 4073 static const char _BesselY_s []="BesselY"; 4074 static define_unary_function_eval (__BesselY,&_BesselY,_BesselY_s); 4075 define_unary_function_ptr5( at_BesselY ,alias_at_BesselY,&__BesselY,0,true); 4076 _besselI(const gen & g,GIAC_CONTEXT)4077 gen _besselI(const gen & g,GIAC_CONTEXT){ 4078 if ( g.type==_STRNG && g.subtype==-1) return g; 4079 return bessel(g,0,contextptr); 4080 } 4081 static const char _besselI_s []="besselI"; 4082 static define_unary_function_eval (__besselI,&_besselI,_besselI_s); 4083 define_unary_function_ptr5( at_besselI ,alias_at_besselI,&__besselI,0,true); 4084 _besselJ(const gen & g,GIAC_CONTEXT)4085 gen _besselJ(const gen & g,GIAC_CONTEXT){ 4086 if ( g.type==_STRNG && g.subtype==-1) return g; 4087 return bessel(g,1,contextptr); 4088 } 4089 static const char _besselJ_s []="besselJ"; 4090 static define_unary_function_eval (__besselJ,&_besselJ,_besselJ_s); 4091 define_unary_function_ptr5( at_besselJ ,alias_at_besselJ,&__besselJ,0,true); 4092 _besselK(const gen & g,GIAC_CONTEXT)4093 gen _besselK(const gen & g,GIAC_CONTEXT){ 4094 if ( g.type==_STRNG && g.subtype==-1) return g; 4095 return bessel(g,2,contextptr); 4096 } 4097 static const char _besselK_s []="besselK"; 4098 static define_unary_function_eval (__besselK,&_besselK,_besselK_s); 4099 define_unary_function_ptr5( at_besselK ,alias_at_besselK,&__besselK,0,true); 4100 _besselY(const gen & g,GIAC_CONTEXT)4101 gen _besselY(const gen & g,GIAC_CONTEXT){ 4102 if ( g.type==_STRNG && g.subtype==-1) return g; 4103 return bessel(g,3,contextptr); 4104 } 4105 static const char _besselY_s []="besselY"; 4106 static define_unary_function_eval (__besselY,&_besselY,_besselY_s); 4107 define_unary_function_ptr5( at_besselY ,alias_at_besselY,&__besselY,0,true); 4108 4109 // sum_{k=1}^n 1/k^e _harmonic(const gen & g,GIAC_CONTEXT)4110 gen _harmonic(const gen & g,GIAC_CONTEXT){ 4111 if ( g.type==_STRNG && g.subtype==-1) return g; 4112 gen e=1; 4113 gen n=g; 4114 if (g.type==_VECT && g.subtype==_SEQ__VECT && g._VECTptr->size()==2){ 4115 e=g._VECTptr->front(); 4116 n=g._VECTptr->back(); 4117 } 4118 if (n==plus_inf) 4119 return Zeta(e,contextptr); 4120 if (!is_integral(n)) 4121 return symbolic(at_harmonic,g); 4122 if (is_greater(0,n,contextptr) || is_greater(n,1e7,contextptr)) 4123 return gendimerr(contextptr); 4124 gen res=1; 4125 for (int k=2;k<=n.val;++k){ 4126 res += plus_one/pow(gen(k),e,contextptr); 4127 } 4128 return res; 4129 } 4130 static const char _harmonic_s []="harmonic"; 4131 static define_unary_function_eval (__harmonic,&_harmonic,_harmonic_s); 4132 define_unary_function_ptr5( at_harmonic ,alias_at_harmonic,&__harmonic,0,true); 4133 4134 #include "input_parser.h" 4135 _constants_catalog(const gen & g,GIAC_CONTEXT)4136 gen _constants_catalog(const gen & g,GIAC_CONTEXT){ 4137 if ( g.type==_STRNG && g.subtype==-1) return g; 4138 if (g.type!=_STRNG) 4139 return undef; 4140 const string & gs=*g._STRNGptr; 4141 gen G; 4142 if (gs=="black"){ 4143 G.val=_BLACK; 4144 G.subtype=_INT_COLOR; 4145 return G; 4146 } 4147 if (gs=="white"){ 4148 G.val=_WHITE; 4149 G.subtype=_INT_COLOR; 4150 return G; 4151 } 4152 if (gs=="red"){ 4153 G.val=_RED; 4154 G.subtype=_INT_COLOR; 4155 return G; 4156 } 4157 if (gs=="green"){ 4158 G.val=_GREEN; 4159 G.subtype=_INT_COLOR; 4160 return G; 4161 } 4162 if (gs=="blue"){ 4163 G.val=_BLUE; 4164 G.subtype=_INT_COLOR; 4165 return G; 4166 } 4167 if (gs=="yellow"){ 4168 G.val=_YELLOW; 4169 G.subtype=_INT_COLOR; 4170 return G; 4171 } 4172 if (gs=="magenta"){ 4173 G.val=_MAGENTA; 4174 G.subtype=_INT_COLOR; 4175 return G; 4176 } 4177 if (gs=="cyan"){ 4178 G.val=_CYAN; 4179 G.subtype=_INT_COLOR; 4180 return G; 4181 } 4182 if (gs=="filled"){ 4183 G.val=_FILL_POLYGON; 4184 G.subtype=_INT_COLOR; 4185 return G; 4186 } 4187 if (gs=="hidden_name"){ 4188 G.val=_HIDDEN_NAME; 4189 G.subtype=_INT_COLOR; 4190 return G; 4191 } 4192 return -1; 4193 } 4194 static const char _black_s []="black"; 4195 static define_unary_function_eval (__black,&_constants_catalog,_black_s); 4196 define_unary_function_ptr5( at_black ,alias_at_black,&__black,0,T_NUMBER); 4197 4198 static const char _white_s []="white"; 4199 static define_unary_function_eval (__white,&_constants_catalog,_white_s); 4200 define_unary_function_ptr5( at_white ,alias_at_white,&__white,0,T_NUMBER); 4201 4202 static const char _red_s []="red"; 4203 static define_unary_function_eval (__red,&_constants_catalog,_red_s); 4204 define_unary_function_ptr5( at_red ,alias_at_red,&__red,0,T_NUMBER); 4205 4206 static const char _green_s []="green"; 4207 static define_unary_function_eval (__green,&_constants_catalog,_green_s); 4208 define_unary_function_ptr5( at_green ,alias_at_green,&__green,0,T_NUMBER); 4209 4210 static const char _blue_s []="blue"; 4211 static define_unary_function_eval (__blue,&_constants_catalog,_blue_s); 4212 define_unary_function_ptr5( at_blue ,alias_at_blue,&__blue,0,T_NUMBER); 4213 4214 static const char _cyan_s []="cyan"; 4215 static define_unary_function_eval (__cyan,&_constants_catalog,_cyan_s); 4216 define_unary_function_ptr5( at_cyan ,alias_at_cyan,&__cyan,0,T_NUMBER); 4217 4218 static const char _yellow_s []="yellow"; 4219 static define_unary_function_eval (__yellow,&_constants_catalog,_yellow_s); 4220 define_unary_function_ptr5( at_yellow ,alias_at_yellow,&__yellow,0,T_NUMBER); 4221 4222 static const char _magenta_s []="magenta"; 4223 static define_unary_function_eval (__magenta,&_constants_catalog,_magenta_s); 4224 define_unary_function_ptr5( at_magenta ,alias_at_magenta,&__magenta,0,T_NUMBER); 4225 4226 static const char _filled_s []="filled"; 4227 static define_unary_function_eval (__filled,&_constants_catalog,_filled_s); 4228 define_unary_function_ptr5( at_filled ,alias_at_filled,&__filled,0,T_NUMBER); 4229 4230 static const char _hidden_name_s []="hidden_name"; 4231 static define_unary_function_eval (__hidden_name,&_constants_catalog,_hidden_name_s); 4232 define_unary_function_ptr5( at_hidden_name ,alias_at_hidden_name,&__hidden_name,0,T_NUMBER); 4233 4234 #if defined(VISUALC) || defined(BESTA_OS) 4235 const double M_E=2.7182818284590452; 4236 #endif 4237 LambertW(complex<double> z,int n)4238 complex<double> LambertW(complex<double> z,int n){ 4239 // n!=0 is not implemented yet 4240 if (z==0) return z; 4241 complex<double> w; 4242 // initial guess 4243 w=2.0*(M_E*z+1.0); 4244 if (std::abs(w)<0.1 && (n==0 || ( n==1 && z.imag()<0) || (n==-1 && z.imag()>0))){ 4245 // near -1/e, set p=sqrt(2(ez+1)), -1+p-1/3*p^2+11/72*p^3+... 4246 w=std::sqrt(w); 4247 if (n==0) w=-1.0+w*(1.0+w*(-1./3.+w*11./72.)); 4248 if (n==1 || n==-1) w=-1.0+w*(-1.0+w*(-1./3.-w*11./72.)); 4249 } 4250 else { 4251 if (z.imag()==0 && z.real()<1 && w.real()>0 && (n==0 || n==-1)){ 4252 w=1; 4253 if (n==-1 && z.real()<0){ 4254 double lnw=std::log(-z.real()); 4255 double lnlnw=std::log(-lnw); 4256 w=lnw-lnlnw-lnlnw/lnw; 4257 } 4258 } 4259 else { 4260 // almost everywhere Log(z)-ln(Log(z)) 4261 w=std::log(z)+2.0*n*complex<double>(0,M_PI); 4262 if (std::abs(z)>=3) 4263 w=w-std::log(w); 4264 } 4265 } 4266 if (n==0 && std::abs(z - .5)<=.5) 4267 w = (0.35173371 * (0.1237166 + 7.061302897 * z)) / (2. + 0.827184 * (1. + 2. * z));// (1,1) Pade approximant for W(z,0) 4268 if (n==-1 && std::abs(z - .5)<=.5) 4269 w = -((complex<double>(2.2591588985 ,4.22096) * (complex<double>(-14.073271 ,-33.767687754) * z - complex<double>(12.7127,-19.071643) * (1. + 2.*z))) / (2. - complex<double>(17.23103,-10.629721) * (1. + 2.*z)));// (1,1) Pade 4270 if (z.imag()==0 && w.imag()==0){ 4271 double Z=z.real(),W=w.real(); 4272 while (1){ 4273 // wnext=w-(w*exp(w)-z)/(exp(w)*(w+1)-(w+2)*(w*exp(w)-z)/(2*w+2)) 4274 double expw(std::exp(W)),wexpwz(W*expw-Z),w1(W+1.0); 4275 double wnext(W-wexpwz/(w1*expw-(W+2.0)*wexpwz/w1/2.0)); 4276 if (abs(wnext-W)<1e-13*std::abs(W)) 4277 return wnext; 4278 W=wnext; 4279 } 4280 } 4281 while (1){ 4282 // wnext=w-(w*exp(w)-z)/(exp(w)*(w+1)-(w+2)*(w*exp(w)-z)/(2*w+2)) 4283 complex<double> expw(std::exp(w)),wexpwz(w*expw-z),w1(w+1.0); 4284 complex<double> wnext(w-wexpwz/(w1*expw-(w+2.0)*wexpwz/w1/2.0)); 4285 if (abs(wnext-w)<1e-13*(1.0+std::abs(w))) 4286 return wnext; 4287 w=wnext; 4288 } 4289 } 4290 4291 #ifdef HAVE_LIBMPFR LambertW(const gen & Z,int n)4292 gen LambertW(const gen & Z,int n){ 4293 gen z(Z); 4294 if (z==0) return z; 4295 int nbits=45; 4296 if (z.type==_REAL) 4297 nbits=mpfr_get_prec(z._REALptr->inf); 4298 if (z.type==_CPLX && z._CPLXptr->type==_REAL) 4299 nbits=mpfr_get_prec(z._CPLXptr->_REALptr->inf); 4300 // initial guess 4301 gen w=evalf_double(z,1,context0); 4302 if (w.type==_DOUBLE_) 4303 w=LambertW(complex<double>(w._DOUBLE_val,0),n); 4304 else { 4305 if (w.type!=_CPLX || w.subtype!=3) 4306 return gensizeerr("Unable to convert to float"); 4307 w=LambertW(complex<double>(w._CPLXptr->_DOUBLE_val,(w._CPLXptr+1)->_DOUBLE_val)); 4308 } 4309 if (nbits<=45) 4310 return w; 4311 int addprec=10; 4312 gen tmp=abs(w,context0); 4313 if (is_greater(tmp,1,context0)) 4314 addprec += int(std::floor(evalf_double(ln(tmp,context0),1,context0)._DOUBLE_val)); 4315 w=accurate_evalf(w,nbits+addprec); 4316 z=accurate_evalf(z,nbits+addprec); 4317 gen eps=accurate_evalf(pow(inv(2,context0),nbits,context0),100);//std::pow(.5,nbits); 4318 while (1){ 4319 // wnext=w-(w*exp(w)-z)/(exp(w)*(w+1)-(w+2)*(w*exp(w)-z)/(2*w+2)) 4320 gen expw(exp(w,context0)),wexpwz(w*expw-z),w1(w+1); 4321 gen wnext(w-wexpwz/(w1*expw-(w+2)*wexpwz/w1/2)); 4322 if (is_greater(eps*(1+abs(w,context0)),abs(wnext-w,context0),context0)) 4323 return accurate_evalf(wnext,nbits); 4324 w=wnext; 4325 } 4326 } 4327 #endif 4328 4329 #ifndef NO_NAMESPACE_GIAC 4330 } // namespace giac 4331 #endif // ndef NO_NAMESPACE_GIAC 4332