1 /* -*- mode:C++ ; compile-command: "g++-3.4 -I.. -g -c permu.cc -DHAVE_CONFIG_H -DIN_GIAC" -*- */ 2 #include "giacPCH.h" 3 /* 4 * Copyright (C) 2005, 2007 R. De Graeve & B. Parisse, 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 20 using namespace std; 21 #include "permu.h" 22 #include "usual.h" 23 #include "sym2poly.h" 24 #include "rpn.h" 25 #include "prog.h" 26 #include "derive.h" 27 #include "subst.h" 28 #include "misc.h" 29 #include "plot.h" 30 #include "intg.h" 31 #include "ifactor.h" 32 #include "lin.h" 33 #include "modpoly.h" 34 #include "giacintl.h" 35 36 #ifndef NO_NAMESPACE_GIAC 37 namespace giac { 38 #endif // ndef NO_NAMESPACE_GIAC 39 vector_int_2_vecteur(const vector<int> & v,GIAC_CONTEXT)40 vecteur vector_int_2_vecteur(const vector<int> & v,GIAC_CONTEXT){ 41 //transforme un vector<int> en vecteur 42 vector<int>::const_iterator it=v.begin(),itend=v.end(); 43 vecteur res; 44 res.reserve(itend-it); 45 if (array_start(contextptr)){ //(xcas_mode(contextptr) || abs_calc_mode(contextptr)==38){ 46 for (;it!=itend;++it) 47 res.push_back(*it+1); 48 } 49 else { 50 for (;it!=itend;++it) 51 res.push_back(*it); 52 } 53 return res; 54 } 55 vector_int_2_vecteur(const vector<int> & v)56 vecteur vector_int_2_vecteur(const vector<int> & v){ 57 //transforme un vector<int> en vecteur 58 vector<int>::const_iterator it=v.begin(),itend=v.end(); 59 vecteur res; 60 res.reserve(itend-it); 61 for (;it!=itend;++it) 62 res.push_back(*it); 63 return res; 64 } 65 vecteur_2_vector_int(const vecteur & v,GIAC_CONTEXT)66 static vector<int> vecteur_2_vector_int(const vecteur & v,GIAC_CONTEXT){ 67 //transforme un vecteur en vector<int> -> empty vector on error 68 vecteur::const_iterator it=v.begin(),itend=v.end(); 69 vector<int> res; 70 res.reserve(itend-it); 71 if (array_start(contextptr)){ //(xcas_mode(contextptr) || abs_calc_mode(contextptr)==38){ 72 for (;it!=itend;++it) 73 if ((*it).type==_INT_) 74 res.push_back((*it).val-1); 75 else 76 return vector<int>(0); 77 } 78 else { 79 for (;it!=itend;++it) 80 if ((*it).type==_INT_) 81 res.push_back((*it).val); 82 else 83 return vector<int>(0); 84 } 85 return res; 86 } 87 vecteur_2_vector_int(const vecteur & v)88 vector<int> vecteur_2_vector_int(const vecteur & v){ 89 //transforme un vecteur en vector<int> -> empty vector on error 90 vecteur::const_iterator it=v.begin(),itend=v.end(); 91 vector<int> res; 92 res.reserve(itend-it); 93 for (;it!=itend;++it){ 94 if ((*it).type==_INT_) 95 res.push_back((*it).val); 96 else 97 return vector<int>(0); 98 } 99 return res; 100 } 101 vecteur_2_vectvector_int(const vecteur & v,GIAC_CONTEXT)102 static vector< vector<int> > vecteur_2_vectvector_int(const vecteur & v,GIAC_CONTEXT){ 103 //transforme un vecteur en vector< vector<int> > -> empty vector on error 104 vecteur::const_iterator it=v.begin(),itend=v.end(); 105 vector< vector<int> > res; 106 res.reserve(itend-it); 107 for (;it!=itend;++it){ 108 if (it->type!=_VECT) 109 return vector< vector<int> >(0); 110 res.push_back(vecteur_2_vector_int(*it->_VECTptr,contextptr)); 111 } 112 return res; 113 } 114 vecteur_2_vectvector_int(const vecteur & v)115 vector< vector<int> > vecteur_2_vectvector_int(const vecteur & v){ 116 //transforme un vecteur en vector< vector<int> > -> empty vector on error 117 vecteur::const_iterator it=v.begin(),itend=v.end(); 118 vector< vector<int> > res; 119 res.reserve(itend-it); 120 for (;it!=itend;++it){ 121 if (it->type!=_VECT) 122 return vector< vector<int> >(0); 123 res.push_back(vecteur_2_vector_int(*it->_VECTptr)); 124 } 125 return res; 126 } 127 vectvector_int_2_vecteur(const vector<vector<int>> & v,GIAC_CONTEXT)128 static vecteur vectvector_int_2_vecteur(const vector< vector<int> > & v,GIAC_CONTEXT){ 129 //transforme un vector< vector<int> > en vecteur 130 int s=int(v.size()); 131 vecteur res; 132 res.reserve(s); 133 for (int i=0;i<s;++i) 134 res.push_back(vector_int_2_vecteur(v[i],contextptr)); 135 return res; 136 } 137 vectvector_int_2_vecteur(const vector<vector<int>> & v)138 vecteur vectvector_int_2_vecteur(const vector< vector<int> > & v){ 139 //transforme un vector< vector<int> > en vecteur 140 int s=int(v.size()); 141 vecteur res; 142 res.reserve(s); 143 for (int i=0;i<s;++i) 144 res.push_back(vector_int_2_vecteur(v[i])); 145 return res; 146 } 147 sizes(const vector<vector<int>> & v)148 vector<int> sizes(const vector< vector<int> > & v){ 149 //donne la liste des tailles des vecteurs qui forment v 150 int s=int(v.size()); 151 vector<int> res(s); 152 for (int i=0;i<s;i++){ 153 vector<int> vi; 154 vi=v[i]; 155 res[i]=int(vi.size()); 156 //res.push_back(vi.size());pourqoi? 157 } 158 return res; 159 } 160 _sizes(const gen & args,GIAC_CONTEXT)161 gen _sizes(const gen & args,GIAC_CONTEXT){ 162 if ( args.type==_STRNG && args.subtype==-1) return args; 163 if (args.type!=_VECT) 164 return gentypeerr(contextptr); 165 vecteur v(*args._VECTptr); 166 vecteur res; 167 vecteur::const_iterator it=v.begin(),itend=v.end(); 168 res.reserve(itend-it); 169 for (;it!=itend;++it){ 170 if (it->type!=_VECT) 171 return gensizeerr(contextptr); 172 res.push_back(int(it->_VECTptr->size())); 173 } 174 return res; 175 } 176 static const char _sizes_s[]="sizes"; 177 static define_unary_function_eval (__sizes,&_sizes,_sizes_s); 178 define_unary_function_ptr5( at_sizes ,alias_at_sizes,&__sizes,0,true); 179 lcm(int a,int b)180 longlong lcm(int a,int b){ 181 int d=gcd(a,b); 182 return (a/d)*longlong(b); 183 } 184 lcm(const vector<int> & v)185 int lcm(const vector<int> & v){ 186 vector<int>::const_iterator it=v.begin(),itend=v.end(); 187 if (itend==it) return 0; 188 if (itend==it+1) return *it; 189 int r=lcm(*it,*(it+1)); 190 for (it+=2;it<itend;++it) 191 r=lcm(r,*it); 192 return r; 193 } 194 cyclesorder(const vector<vector<int>> & v,GIAC_CONTEXT)195 static int cyclesorder(const vector< vector<int> > & v,GIAC_CONTEXT){ 196 vector<int> ss; 197 ss=sizes(v); 198 return lcm(ss); 199 } permuorder(const vector<int> & p,GIAC_CONTEXT)200 static int permuorder(const vector<int> & p,GIAC_CONTEXT){ 201 vector< vector<int> > c; 202 c=permu2cycles(p); 203 return lcm(sizes(c)); 204 } 205 _permuorder(const gen & args,GIAC_CONTEXT)206 gen _permuorder(const gen & args,GIAC_CONTEXT){ 207 if ( args.type==_STRNG && args.subtype==-1) return args; 208 if (args.type!=_VECT) 209 return gensizeerr(contextptr); 210 vecteur v(*args._VECTptr); 211 vector<int> p; 212 if (!is_permu(v,p,contextptr)) 213 return gensizeerr(contextptr); 214 return (permuorder(p,contextptr)); 215 } 216 static const char _permuorder_s[]="permuorder"; 217 static define_unary_function_eval (__permuorder,&_permuorder,_permuorder_s); 218 define_unary_function_ptr5( at_permuorder ,alias_at_permuorder,&__permuorder,0,true); 219 shuffle(vector<int> & temp)220 void shuffle(vector<int> & temp){ 221 int n=int(temp.size()); 222 // source wikipedia Fisher-Yates shuffle article 223 for (int i=0;i<n-1;++i){ 224 // j ← random integer such that i ≤ j < n 225 // exchange a[i] and a[j] 226 int j=i+int((std_rand()/(rand_max2+1.0))*(n-i)); 227 std::swap(temp[i],temp[j]); 228 } 229 } 230 rand_k_n(int k,int n,bool sorted)231 vector<int> rand_k_n(int k,int n,bool sorted){ 232 if (k<=0 || n<=0) 233 return vector<int>(0); 234 if (//n>=65536 && 235 k*double(k)<=n/4){ 236 vector<int> t(k),ts(k); 237 for (int essai=20;essai>=0;--essai){ 238 int i; 239 for (i=0;i<k;++i) 240 ts[i]=t[i]=int(std_rand()/(rand_max2+1.0)*n); 241 sort(ts.begin(),ts.end()); 242 for (i=1;i<k;++i){ 243 if (ts[i]==ts[i-1]) 244 break; 245 } 246 if (i==k) 247 return sorted?ts:t; 248 } 249 } 250 if (k>=n/3 || (sorted && k*std::log(double(k))>n) ){ 251 vector<int> t; t.reserve(k); 252 // (algorithm suggested by O. Garet) 253 while (n>0){ 254 int r=int(std_rand()/(rand_max2+1.0)*n); 255 if (r<n-k) // (n-k)/n=proba that the current n is not in the list 256 --n; 257 else { 258 --n; 259 t.push_back(n); 260 --k; 261 } 262 } 263 if (sorted) 264 reverse(t.begin(),t.end()); 265 else 266 shuffle(t); 267 return t; 268 } 269 vector<bool> tab(n,true); 270 vector<int> v(k); 271 for (int j=0;j<k;++j){ 272 int r=-1; 273 for (;;){ 274 r=int(std_rand()/(rand_max2+1.0)*n); 275 if (tab[r]) break; 276 } 277 v[j]=r; 278 } 279 if (sorted) 280 sort(v.begin(),v.end()); 281 return v; 282 } 283 randperm(const int & n)284 vector<int> randperm(const int & n){ 285 //renvoie une permutation au hasard de long n 286 vector<int> temp(n); 287 for (int k=0;k<n;k++) temp[k]=k; 288 #if 1 289 shuffle(temp); 290 return temp; 291 #else 292 vector<int> p(n); 293 //on chosit au hasard h et alors p[k]=temp[h] 294 int m; 295 m=n; 296 for (int k=0;k<n;k++) { 297 int h=int(std::floor(m*(std_rand()/(rand_max2+1.0)))); 298 p[k]=temp[h]; m=m-1; 299 //mise a jour de temp :il faut supprimer temp[h] 300 for (int j=h;j<m;j++) { 301 temp[j]=temp[j+1]; 302 } 303 } 304 return(p); 305 #endif 306 } _randperm(const gen & args,GIAC_CONTEXT)307 gen _randperm(const gen & args,GIAC_CONTEXT){ 308 if ( args.type==_STRNG && args.subtype==-1) return args; 309 if (args.type==_VECT){ 310 vecteur & v = *args._VECTptr; 311 vector<int> p=randperm(int(v.size())); 312 vecteur w(v); 313 for (unsigned i=0;i<v.size();++i){ 314 w[i]=v[p[i]]; 315 } 316 return w; 317 } 318 gen n=args; 319 if (!is_integral(n) || n.type!=_INT_ || n.val<=0) 320 return gensizeerr(contextptr); 321 return vector_int_2_vecteur(randperm(n.val),contextptr); 322 } 323 static const char _randperm_s[]="randperm"; 324 static define_unary_function_eval (__randperm,&_randperm,_randperm_s); 325 define_unary_function_ptr5( at_randperm ,alias_at_randperm,&__randperm,0,true); 326 is_permu(const vecteur & p,vector<int> & p1,GIAC_CONTEXT)327 bool is_permu(const vecteur &p,vector<int> & p1,GIAC_CONTEXT) { 328 //renvoie true si p est une perm et transforme p en le vector<int> p1 329 int n; 330 n=int(p.size()); 331 vector<int> p2(n); 332 p1=p2; 333 vector<int> temp(n); 334 335 for (int j=0;j<n;j++){ if (p[j].type!=_INT_){return(false);}} 336 337 for (int j=0;j<n;j++){ 338 if (array_start(contextptr)) //xcas_mode(contextptr)>0 || abs_calc_mode(contextptr)==38) 339 p1[j]=p[j].val-1; 340 else 341 p1[j]=p[j].val; 342 if ((n<=p1[j])|| (p1[j])<0) { 343 return(false); 344 } 345 } 346 int k; 347 k=0; 348 while (k<n) { 349 int p1k=p1[k]; 350 if (p1k<0 || p1k>=n) {return(false);} 351 if (temp[p1k]) { 352 return(false);} 353 else {temp[p1k]=1;} 354 k=k+1; 355 } 356 return(true); 357 } _is_permu(const gen & args,GIAC_CONTEXT)358 gen _is_permu(const gen & args,GIAC_CONTEXT){ 359 if ( args.type==_STRNG && args.subtype==-1) return args; 360 if (args.type!=_VECT) 361 return gensizeerr(contextptr); 362 vecteur v(*args._VECTptr); 363 vector<int> p1; 364 return is_permu(v,p1,contextptr); 365 } 366 static const char _is_permu_s[]="is_permu"; 367 static define_unary_function_eval (__is_permu,&_is_permu,_is_permu_s); 368 define_unary_function_ptr5( at_is_permu ,alias_at_is_permu,&__is_permu,0,true); 369 370 is_cycle(const vecteur & c,vector<int> & c1,GIAC_CONTEXT)371 bool is_cycle(const vecteur & c,vector<int> & c1,GIAC_CONTEXT) { 372 //renvoie true si c est un cycle et transf c vecteur en le cycle c1 vector<int> 373 int n1; 374 n1=int(c.size()); 375 //vector<int> p; 376 c1.resize(n1); 377 int b=array_start(contextptr); 378 for (int j=0;j<n1;j++){ 379 int tmp=c[j].val-b; 380 if (tmp<0) 381 return false; 382 c1[j]=tmp; 383 } 384 vector<int> c2(c1); 385 sort(c2.begin(),c2.end()); 386 for (int k=1;k<n1;k++) { 387 if (c2[k]==c2[k-1]) 388 return false; 389 } 390 return true; 391 } _is_cycle(const gen & args,GIAC_CONTEXT)392 gen _is_cycle(const gen & args,GIAC_CONTEXT){ 393 if ( args.type==_STRNG && args.subtype==-1) return args; 394 if (args.type!=_VECT) 395 return gensizeerr(contextptr); 396 vecteur v(*args._VECTptr); 397 vector<int> c1; 398 return is_cycle(v,c1,contextptr); 399 } 400 static const char _is_cycle_s[]="is_cycle"; 401 static define_unary_function_eval (__is_cycle,&_is_cycle,_is_cycle_s); 402 define_unary_function_ptr5( at_is_cycle ,alias_at_is_cycle,&__is_cycle,0,true); 403 cycle2perm(const vector<int> & c)404 vector<int> cycle2perm(const vector<int> & c) { 405 //transforme c en la permu p et renvoie p 406 int n1; 407 n1=int(c.size()); 408 //vector<int> c1(n1); 409 int n; 410 n=c[0]; 411 for (int k=1;k<n1;k++) { 412 if (n<c[k]) { 413 n=c[k]; 414 } 415 } 416 n=n+1; 417 vector<int> p(n); 418 for (int k=0;k<n;k++) { 419 p[k]=k;} 420 for (int k=0;k<n1-1;k++) {p[c[k]]=c[k+1];} 421 p[c[n1-1]]=c[0]; 422 return(p); 423 } 424 _cycle2perm(const gen & args,GIAC_CONTEXT)425 gen _cycle2perm(const gen & args,GIAC_CONTEXT){ 426 if ( args.type==_STRNG && args.subtype==-1) return args; 427 if (args.type!=_VECT) 428 return gensizeerr(contextptr); 429 vecteur v(*args._VECTptr); 430 vector<int> c; 431 if (!is_cycle(v,c,contextptr)) 432 return gensizeerr(contextptr); 433 return vector_int_2_vecteur(cycle2perm(c),contextptr); 434 } 435 static const char _cycle2perm_s[]="cycle2perm"; 436 static define_unary_function_eval (__cycle2perm,&_cycle2perm,_cycle2perm_s); 437 define_unary_function_ptr5( at_cycle2perm ,alias_at_cycle2perm,&__cycle2perm,0,true); 438 439 p1op2(const vector<int> & p1,const vector<int> & p2)440 vector<int> p1op2(const vector<int> & p1,const vector<int> & p2) { 441 //composition de 2 perm p1 et p2 de long n1 et n2 : a pour long max(n1,n2) 442 int n1; 443 n1=int(p1.size()); 444 int n2; 445 n2=int(p2.size()); 446 vector<int> p3; 447 vector<int> p4; 448 p3=p1; 449 p4=p2; 450 if (n1>n2) { 451 for (int k=n2;k<n1;k++) {p4.push_back(k);n2=n1;} 452 } else { 453 for (int k=n1;k<n2;k++) {p3.push_back(k);n1=n2;} 454 } 455 ; 456 vector<int> p(n1); 457 for (int k=0;k<n1;k++) { 458 p[k]=p3[p4[k]]; 459 } 460 return(p); 461 } 462 _p1op2(const gen & args,GIAC_CONTEXT)463 gen _p1op2(const gen & args,GIAC_CONTEXT){ 464 if ( args.type==_STRNG && args.subtype==-1) return args; 465 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 466 return gentypeerr(contextptr); 467 vecteur v(*args._VECTptr); 468 gen v1=v.front(),v2=v.back(); 469 if ( (v1.type!=_VECT) || (v2.type!=_VECT)) 470 return gentypeerr(contextptr); 471 vector<int> p1,p2; 472 if (!is_permu(*v1._VECTptr,p1,contextptr) || !is_permu(*v2._VECTptr,p2,contextptr)) 473 return gensizeerr(contextptr); 474 return vector_int_2_vecteur(p1op2(p1,p2),contextptr); 475 } 476 static const char _p1op2_s[]="p1op2"; 477 static define_unary_function_eval (__p1op2,&_p1op2,_p1op2_s); 478 define_unary_function_ptr5( at_p1op2 ,alias_at_p1op2,&__p1op2,0,true); 479 c1oc2(const vector<int> & c1,const vector<int> & c2)480 vector<int> c1oc2(const vector<int> & c1,const vector<int> & c2) { 481 //composition de 2 cycles en une perm de long min 482 vector<int> p1; 483 p1=cycle2perm(c1); 484 vector<int> p2; 485 p2=cycle2perm(c2); 486 int n1; 487 n1=int(p1.size()); 488 int n2; 489 n2=int(p2.size()); 490 int n; 491 if (n1>n2) { 492 n=n1; 493 for (int k=n2;k<n;k++) {p2.push_back(k);} 494 } 495 else { 496 n=n2; 497 for (int k=n1;k<n;k++) {p1.push_back(k);} 498 } 499 vector<int> p(n); 500 for (int k=0;k<n;k++) { 501 p[k]=p1[p2[k]]; 502 } 503 return(p); 504 } 505 _c1oc2(const gen & args,GIAC_CONTEXT)506 gen _c1oc2(const gen & args,GIAC_CONTEXT){ 507 if ( args.type==_STRNG && args.subtype==-1) return args; 508 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 509 return gentypeerr(contextptr); 510 vecteur v(*args._VECTptr); 511 gen v1=v.front(),v2=v.back(); 512 if ( (v1.type!=_VECT) || (v2.type!=_VECT)) 513 return gentypeerr(contextptr); 514 vector<int> c1,c2; 515 //c1=vecteur_2_vector_int(*v1._VECTptr); 516 //c2=vecteur_2_vector_int(*v2._VECTptr); 517 if (!is_cycle(*v1._VECTptr,c1,contextptr) || !is_cycle(*v2._VECTptr,c2,contextptr)) 518 return gensizeerr(contextptr); 519 return vector_int_2_vecteur(c1oc2(c1,c2),contextptr); 520 } 521 static const char _c1oc2_s[]="c1oc2"; 522 static define_unary_function_eval (__c1oc2,&_c1oc2,_c1oc2_s); 523 define_unary_function_ptr5( at_c1oc2 ,alias_at_c1oc2,&__c1oc2,0,true); 524 c1op2(const vector<int> & c1,const vector<int> & p2)525 vector<int> c1op2(const vector<int> & c1, const vector<int> & p2) { 526 //composition d'un cycle et d'une perm en une perm de long min 527 vector<int> p1,p3; 528 p1=cycle2perm(c1); 529 int n1; 530 n1=int(p1.size()); 531 int n2; 532 n2=int(p2.size()); 533 p3=p2; 534 int n; 535 if (n1>n2) { 536 n=n1; 537 for (int k=n2;k<n;k++) {p3.push_back(k);} 538 } 539 else { 540 n=n2; 541 for (int k=n1;k<n;k++) {p1.push_back(k);} 542 } 543 vector<int> p(n); 544 for (int k=0;k<n;k++) { 545 p[k]=p1[p3[k]]; 546 } 547 return(p); 548 } 549 _c1op2(const gen & args,GIAC_CONTEXT)550 gen _c1op2(const gen & args,GIAC_CONTEXT){ 551 if ( args.type==_STRNG && args.subtype==-1) return args; 552 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 553 return gentypeerr(contextptr); 554 vecteur v(*args._VECTptr); 555 gen v1=v.front(),v2=v.back(); 556 if ( (v1.type!=_VECT) || (v2.type!=_VECT)) 557 return gentypeerr(contextptr); 558 vector<int> c1,p2; 559 if (!is_cycle(*v1._VECTptr,c1,contextptr) || !is_permu(*v2._VECTptr,p2,contextptr)) 560 return gensizeerr(contextptr); 561 return vector_int_2_vecteur(c1op2(c1,p2),contextptr); 562 } 563 static const char _c1op2_s[]="c1op2"; 564 static define_unary_function_eval (__c1op2,&_c1op2,_c1op2_s); 565 define_unary_function_ptr5( at_c1op2 ,alias_at_c1op2,&__c1op2,0,true); 566 p1oc2(const vector<int> & p1,const vector<int> & c2)567 vector<int> p1oc2(const vector<int> & p1, const vector<int> & c2) { 568 //composition d'une perm et d'un cycle en une perm de long min 569 vector<int> p2,p3; 570 p2=cycle2perm(c2); 571 int n2; 572 n2=int(p2.size()); 573 int n1; 574 n1=int(p1.size()); 575 p3=p1; 576 int n; 577 if (n1>n2) {n=n1; 578 for (int k=n2;k<n;k++) {p2.push_back(k);} 579 } else {n=n2; 580 for (int k=n1;k<n;k++) {p3.push_back(k);} 581 } 582 vector<int> p(n); 583 for (int k=0;k<n;k++) { 584 p[k]=p3[p2[k]]; 585 } 586 return(p); 587 } 588 _p1oc2(const gen & args,GIAC_CONTEXT)589 gen _p1oc2(const gen & args,GIAC_CONTEXT){ 590 if ( args.type==_STRNG && args.subtype==-1) return args; 591 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 592 return gentypeerr(contextptr); 593 vecteur v(*args._VECTptr); 594 gen v1=v.front(),v2=v.back(); 595 if ( (v1.type!=_VECT) || (v2.type!=_VECT)) 596 return gentypeerr(contextptr); 597 vector<int> p1,c2; 598 if (!is_cycle(*v2._VECTptr,c2,contextptr) || !is_permu(*v1._VECTptr,p1,contextptr)) 599 return gensizeerr(contextptr); 600 return vector_int_2_vecteur(p1oc2(p1,c2),contextptr); 601 } 602 static const char _p1oc2_s[]="p1oc2"; 603 static define_unary_function_eval (__p1oc2,&_p1oc2,_p1oc2_s); 604 define_unary_function_ptr5( at_p1oc2 ,alias_at_p1oc2,&__p1oc2,0,true); 605 cycles2permu(const vector<vector<int>> & c)606 vector<int> cycles2permu(const vector< vector<int> > & c) { 607 //transforme une liste de cycles en la permutation produit 608 int n; 609 n=int(c.size()); 610 vector<int> pk; 611 vector<int> p; 612 vector<int> ck; 613 vector<int> c1; 614 ck=c[n-1]; 615 vector<int> c0(1); 616 c0[0]=0; 617 p=c1oc2(ck,c0); 618 for (int k=n-2;k>=0;k--){ 619 ck=c[k]; 620 p=c1op2(ck,p); 621 } 622 return(p); 623 } 624 _cycles2permu(const gen & args,GIAC_CONTEXT)625 gen _cycles2permu(const gen & args,GIAC_CONTEXT){ 626 if ( args.type==_STRNG && args.subtype==-1) return args; 627 if (args.type!=_VECT) 628 return gentypeerr(contextptr); 629 vecteur v(*args._VECTptr); 630 vecteur::const_iterator it=v.begin(),itend=v.end(); 631 for (;it!=itend;++it){ 632 vector<int> c; 633 if (it->type!=_VECT || !is_cycle(*it->_VECTptr,c,contextptr)) 634 return gensizeerr(contextptr); 635 } 636 return vector_int_2_vecteur(cycles2permu(vecteur_2_vectvector_int(v,contextptr)),contextptr); 637 } 638 static const char _cycles2permu_s[]="cycles2permu"; 639 static define_unary_function_eval (__cycles2permu,&_cycles2permu,_cycles2permu_s); 640 define_unary_function_ptr5( at_cycles2permu ,alias_at_cycles2permu,&__cycles2permu,0,true); 641 permu2cycles(const vector<int> & p)642 vector< vector<int> > permu2cycles(const vector<int> & p) { 643 //transforme la permutation p en une liste de cycles (p= produit de ces cycles) 644 int l=int(p.size()); 645 int n=0; 646 int deb; 647 vector<int> p1(l); 648 p1=p; 649 vector<int> temp(l+1); 650 vector< vector<int> > c; 651 if (p1[l-1]==l-1) { 652 vector<int> c1; 653 c1.push_back(l-1); 654 c.push_back(c1); l=l-1; 655 } 656 temp[l]=0; 657 for (int k=0;k<l;k++) temp[k]=p1[k]; 658 while (n<l){ 659 vector<int> v; 660 v.push_back(n);deb=n; 661 while (p1[n]!=deb){ 662 v.push_back(p1[n]); 663 temp[n]=0; 664 n=p1[n]; 665 } 666 if (n!=deb) {c.push_back(v);} 667 temp[n]=0; 668 n=deb+1; 669 while ((n<l)&&(temp[n]==0)) n++; 670 } 671 return(c); 672 } 673 _permu2cycles(const gen & args,GIAC_CONTEXT)674 gen _permu2cycles(const gen & args,GIAC_CONTEXT){ 675 if ( args.type==_STRNG && args.subtype==-1) return args; 676 if (args.type!=_VECT) 677 return gensizeerr(contextptr); 678 vecteur v(*args._VECTptr); 679 vector<int> p; 680 if (!is_permu(v,p,contextptr)) 681 return gensizeerr(contextptr); 682 return gen(vectvector_int_2_vecteur(permu2cycles(p),contextptr),_LIST__VECT); 683 } 684 685 static const char _permu2cycles_s[]="permu2cycles"; 686 static define_unary_function_eval (__permu2cycles,&_permu2cycles,_permu2cycles_s); 687 define_unary_function_ptr5( at_permu2cycles ,alias_at_permu2cycles,&__permu2cycles,0,true); 688 perminv(const vector<int> & p)689 vector<int> perminv(const vector<int> & p){ 690 int n; 691 n=int(p.size()); 692 vector<int> p1(n); 693 for (int j=0;j<n;j++){ 694 p1[p[j]]=j; 695 } 696 return p1; 697 } _perminv(const gen & args,GIAC_CONTEXT)698 gen _perminv(const gen & args,GIAC_CONTEXT){ 699 if ( args.type==_STRNG && args.subtype==-1) return args; 700 if (args.type!=_VECT) 701 return gensizeerr(contextptr); 702 vecteur v(*args._VECTptr); 703 vector<int> p; 704 if (!is_permu(v,p,contextptr)) 705 return gensizeerr(contextptr); 706 return vector_int_2_vecteur(perminv(p),contextptr); 707 } 708 static const char _perminv_s[]="perminv"; 709 static define_unary_function_eval (__perminv,&_perminv,_perminv_s); 710 define_unary_function_ptr5( at_perminv ,alias_at_perminv,&__perminv,0,true); 711 cycleinv(const vector<int> & c)712 vector<int> cycleinv(const vector<int> & c){ 713 int n; 714 n=int(c.size()); 715 vector<int> c1(n); 716 for (int j=0;j<n;j++){ 717 c1[j]=c[n-j-1]; 718 } 719 return c1; 720 } _cycleinv(const gen & args,GIAC_CONTEXT)721 gen _cycleinv(const gen & args,GIAC_CONTEXT){ 722 if ( args.type==_STRNG && args.subtype==-1) return args; 723 if (args.type!=_VECT) 724 return gensizeerr(contextptr); 725 vecteur v(*args._VECTptr); 726 vector<int> c; 727 if (!is_cycle(v,c,contextptr)) 728 return gensizeerr(contextptr); 729 return vector_int_2_vecteur(cycleinv(c),contextptr); 730 } 731 static const char _cycleinv_s[]="cycleinv"; 732 static define_unary_function_eval (__cycleinv,&_cycleinv,_cycleinv_s); 733 define_unary_function_ptr5( at_cycleinv ,alias_at_cycleinv,&__cycleinv,0,true); 734 735 signature(const vector<int> & p)736 int signature(const vector<int> & p) { 737 //renvoie la signature de la permutation p 738 int s; 739 s=1; 740 vector< vector<int> > c; 741 c=permu2cycles(p); 742 int l=int(c.size()); 743 for (int k=0;k<l;k++){ 744 int lk=int(c[k].size())-1; 745 if (lk%2) s=s*-1; 746 } 747 return(s); 748 } 749 _signature(const gen & args,GIAC_CONTEXT)750 gen _signature(const gen & args,GIAC_CONTEXT){ 751 if ( args.type==_STRNG && args.subtype==-1) return args; 752 if (args.type!=_VECT) 753 return gensizeerr(contextptr); 754 vecteur v(*args._VECTptr); 755 vector<int> p; 756 if (!is_permu(v,p,contextptr)) 757 return gensizeerr(contextptr); 758 return signature(p); 759 } 760 761 static const char _signature_s[]="signature"; 762 static define_unary_function_eval (__signature,&_signature,_signature_s); 763 define_unary_function_ptr5( at_signature ,alias_at_signature,&__signature,0,true); 764 vector_giac_double_2_vecteur(const vector<giac_double> & v)765 vecteur vector_giac_double_2_vecteur(const vector<giac_double> & v){ 766 //transforme un vector<double> en vecteur 767 vector<giac_double>::const_iterator it=v.begin(),itend=v.end(); 768 vecteur res; 769 res.reserve(itend-it); 770 for (;it!=itend;++it) 771 res.push_back(double(*it)); 772 return res; 773 } 774 vectvector_giac_double_2_vecteur(const vector<vector<giac_double>> & v)775 vecteur vectvector_giac_double_2_vecteur(const vector< vector<giac_double> > & v){ 776 //transforme un vector< vector<double> > en vecteur 777 int s=int(v.size()); 778 vecteur res; 779 res.reserve(s); 780 for (int i=0;i<s;++i) 781 res.push_back(vector_giac_double_2_vecteur(v[i])); 782 return res; 783 } 784 _hilbert(const gen & args,GIAC_CONTEXT)785 gen _hilbert(const gen & args,GIAC_CONTEXT){ 786 if ( args.type==_STRNG && args.subtype==-1) return args; 787 int n,p; 788 if (args.type==_INT_) { 789 n=args.val; 790 p=args.val; 791 } 792 else { 793 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 794 return gentypeerr(contextptr); 795 vecteur v(*args._VECTptr); 796 gen v1=v.front(),v2=v.back(); 797 n=v1.val; 798 p=v2.val; 799 } 800 vecteur c; 801 for (int k=0;k<n;k++){ 802 vecteur l(p); 803 for (int j=0;j<p;j++){ 804 l[j]=rdiv(1,k+j+1,contextptr); 805 } 806 c.push_back(l); 807 } 808 return gen(c,_MATRIX__VECT); 809 } 810 811 static const char _hilbert_s[]="hilbert"; 812 static define_unary_function_eval (__hilbert,&_hilbert,_hilbert_s); 813 define_unary_function_ptr5( at_hilbert ,alias_at_hilbert,&__hilbert,0,true); 814 l2norm2(const gen & g)815 gen l2norm2(const gen & g){ 816 if (g.type!=_VECT) 817 return g*g; 818 const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 819 gen res(0); 820 mpz_t tmpz; 821 mpz_init(tmpz); 822 for (;it!=itend;++it){ 823 if (res.type==_ZINT && is_integer(*it)){ 824 #if defined(INT128) && !defined(USE_GMP_REPLACEMENTS) 825 if (it->type==_INT_) 826 mpz_add_ui(*res._ZINTptr,*res._ZINTptr,longlong(it->val)*(it->val)); 827 else { 828 mpz_mul(tmpz,*it->_ZINTptr,*it->_ZINTptr); 829 mpz_add(*res._ZINTptr,*res._ZINTptr,tmpz); 830 } 831 #else 832 if (it->type==_INT_){ 833 mpz_set_si(tmpz,it->val); 834 mpz_mul(tmpz,tmpz,tmpz); 835 } 836 else 837 mpz_mul(tmpz,*it->_ZINTptr,*it->_ZINTptr); 838 mpz_add(*res._ZINTptr,*res._ZINTptr,tmpz); 839 #endif 840 } 841 else 842 res += (*it)*(*it); 843 } 844 mpz_clear(tmpz); 845 return res; 846 } 847 square_hadamard_bound(const matrice & m)848 gen square_hadamard_bound(const matrice & m){ 849 const_iterateur it=m.begin(),itend=m.end(); 850 gen prod(1); 851 for (;it!=itend;++it){ 852 type_operator_times(prod,l2norm2(*it),prod); 853 // prod=prod*l2norm2(*it); 854 } 855 return prod; 856 } 857 _hadamard(const gen & args,GIAC_CONTEXT)858 gen _hadamard(const gen & args,GIAC_CONTEXT){ 859 if ( args.type==_STRNG && args.subtype==-1) return args; 860 if (ckmatrix(args) || args[0][0].type!=_VECT){ 861 // Hadamard bound on det(args) 862 matrice & m=*args._VECTptr; 863 if (has_num_coeff(m)){ 864 gen r=1.0; 865 for (unsigned i=0;i<m.size();++i) 866 r = r*l2norm(*m[i]._VECTptr,contextptr); 867 return r; 868 } 869 return sqrt(min(square_hadamard_bound(m),square_hadamard_bound(mtran(m)),contextptr),contextptr); 870 } 871 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 872 return gentypeerr(contextptr); 873 vecteur v(*args._VECTptr); 874 gen g1=v.front(),g2=v.back(); 875 876 if ((g1.type!=_VECT) ||(g2.type!=_VECT)) 877 return gentypeerr(contextptr); 878 vecteur v1(*g1._VECTptr); 879 vecteur v2(*g2._VECTptr); 880 if (v1.size()!=v2.size()) return gensizeerr(contextptr); 881 int n=int(v1.size()); 882 vecteur c; 883 for (int k=0;k<n;k++){ 884 if ((v1[k].type!=_VECT) ||(v2[k].type!=_VECT)) return gentypeerr(contextptr); 885 vecteur l1(*(v1[k])._VECTptr); 886 vecteur l2(*(v2[k])._VECTptr); 887 if (l1.size()!=l2.size()) return gensizeerr(contextptr); 888 int p=int(l1.size()); 889 vecteur l(p); 890 for (int j=0;j<p;j++){ 891 l[j]=l1[j]*l2[j]; 892 } 893 894 c.push_back(l); 895 } 896 return gen(c,_MATRIX__VECT); 897 } 898 static const char _hadamard_s[]="hadamard"; 899 static define_unary_function_eval (__hadamard,&_hadamard,_hadamard_s); 900 define_unary_function_ptr5( at_hadamard ,alias_at_hadamard,&__hadamard,0,true); 901 _trn(const gen & args,GIAC_CONTEXT)902 gen _trn(const gen & args,GIAC_CONTEXT){ 903 return conj(_tran(args,contextptr),contextptr); 904 } 905 906 static const char _trn_s[]="trn"; 907 static define_unary_function_eval (__trn,&_trn,_trn_s); 908 define_unary_function_ptr5( at_trn ,alias_at_trn,&__trn,0,true); 909 _syst2mat(const gen & args,GIAC_CONTEXT)910 gen _syst2mat(const gen & args,GIAC_CONTEXT){ 911 if ( args.type==_STRNG && args.subtype==-1) return args; 912 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 913 return gentypeerr(contextptr); 914 vecteur v(*args._VECTptr); 915 gen g1=_equal2diff(v.front(),contextptr),g2=v.back(); 916 917 if ((g1.type!=_VECT) ||(g2.type!=_VECT)) 918 return gentypeerr(contextptr); 919 vecteur v1(*g1._VECTptr); 920 vecteur v2(*g2._VECTptr); 921 922 int n=int(v2.size()),m=int(v1.size()); 923 vecteur c; 924 for (int k=0;k<m;k++){ 925 vecteur l(n+1); 926 gen ln=v1[k]; 927 for (int j=0;j<n;j++){ 928 l[j]=derive(v1[k],v2[j],contextptr); 929 if (is_undef(l[j])) return l[j]; 930 ln=subst(ln,v2[j],0,false,contextptr); 931 } 932 l[n]=ln; 933 c.push_back(l); 934 } 935 return gen(c,_MATRIX__VECT); 936 } 937 static const char _syst2mat_s[]="syst2mat"; 938 static define_unary_function_eval (__syst2mat,&_syst2mat,_syst2mat_s); 939 define_unary_function_ptr5( at_syst2mat ,alias_at_syst2mat,&__syst2mat,0,true); 940 _vandermonde(const gen & args,GIAC_CONTEXT)941 gen _vandermonde(const gen & args,GIAC_CONTEXT){ 942 if ( args.type==_STRNG && args.subtype==-1) return args; 943 if (args.type!=_VECT) 944 return gentypeerr(contextptr); 945 vecteur v(*args._VECTptr); 946 int n=int(v.size()),m=n; 947 if (n==2 && v[1].type==_INT_ && v[0].type==_VECT){ 948 m=v[1].val; 949 v=*v[0]._VECTptr; 950 n=int(v.size()); 951 } 952 vecteur c; 953 vecteur l(n); 954 for (int j=0;j<n;j++){ 955 l[j]=1; 956 } 957 c.push_back(l); 958 for (int k=1;k<m;k++){ 959 for (int j=0;j<n;j++){ 960 l[j]=l[j]*v[j]; 961 } 962 c.push_back(l); 963 } 964 return gen(mtran(c),_MATRIX__VECT); 965 } 966 static const char _vandermonde_s[]="vandermonde"; 967 static define_unary_function_eval (__vandermonde,&_vandermonde,_vandermonde_s); 968 define_unary_function_ptr5( at_vandermonde ,alias_at_vandermonde,&__vandermonde,0,true); 969 970 _laplacian(const gen & args,GIAC_CONTEXT)971 gen _laplacian(const gen & args,GIAC_CONTEXT){ 972 if ( args.type==_STRNG && args.subtype==-1) return args; 973 if (args.type==_DOUBLE_ && args._DOUBLE_val==int(args._DOUBLE_val)){ 974 return evalf(_laplacian(int(args._DOUBLE_val),contextptr),1,context0); 975 } 976 if (args.type==_IDNT){ 977 gen g=eval(args,1,contextptr); 978 if (g.type==_INT_) 979 return _laplacian(g,contextptr); 980 } 981 if (args.type==_INT_ && args.val>0 && args.val<int(std::sqrt(double(LIST_SIZE_LIMIT)))){ 982 // discrete 1-d laplacian matrix 983 int n=args.val; 984 vecteur res(n); 985 for (int i=0;i<n;++i){ 986 vecteur resi(n); 987 if (i) 988 resi[i-1]=-1; 989 resi[i]=2; 990 if (i<n-1) 991 resi[i+1]=-1; 992 res[i]=resi; 993 } 994 return res; 995 } 996 if (args.type==_VECT && args._VECTptr->size()==3){ 997 gen opt=args._VECTptr->back(); 998 if (opt.is_symb_of_sommet(at_equal) && (opt._SYMBptr->feuille[0]==at_coordonnees || (opt._SYMBptr->feuille[0].type==_INT_ && opt._SYMBptr->feuille[0].val==_COORDS))){ 999 gen f=eval(args._VECTptr->front(),1,contextptr); 1000 gen coord=(*args._VECTptr)[1]; 1001 if (coord.type==_VECT && coord._VECTptr->size()==3){ 1002 if (opt._SYMBptr->feuille[1]==at_sphere){ 1003 gen r=coord[0],r2=pow(r,2,contextptr),t=coord[1],s=sin(t,contextptr),p=coord[2]; 1004 gen res=derive(r2*derive(f,r,contextptr),r,contextptr); 1005 res += derive(derive(s*f,t,contextptr),t,contextptr)/s; 1006 res += derive(derive(f,p,contextptr),p,contextptr)/(s*s); 1007 return res/r2; 1008 } 1009 if (opt._SYMBptr->feuille[1]==at_cylindre){ 1010 gen r=coord[0],t=coord[1],z=coord[2]; 1011 gen res=derive(r*derive(f,r,contextptr),r,contextptr); 1012 res += derive(derive(f,t,contextptr),t,contextptr)/(r*r); 1013 res += derive(derive(f,z,contextptr),z,contextptr); 1014 return res; 1015 } 1016 } 1017 } 1018 } 1019 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 1020 return gentypeerr(contextptr); 1021 vecteur v(plotpreprocess(args,contextptr)); 1022 if (is_undef(v)) 1023 return v; 1024 gen g1=v.front(),g2=v.back(); 1025 if (g2.type!=_VECT) 1026 return gentypeerr(contextptr); 1027 vecteur v2(*g2._VECTptr); 1028 int n=int(v2.size()); 1029 gen la; 1030 la=0; 1031 for (int k=0;k<n;k++){ 1032 la=la+derive(derive(g1,v2[k],contextptr),v2[k],contextptr); 1033 } 1034 return normal(la,contextptr); 1035 } 1036 static const char _laplacian_s[]="laplacian"; 1037 static define_unary_function_eval_quoted (__laplacian,&_laplacian,_laplacian_s); 1038 define_unary_function_ptr5( at_laplacian ,alias_at_laplacian,&__laplacian,_QUOTE_ARGUMENTS,true); 1039 _hessian(const gen & args,GIAC_CONTEXT)1040 gen _hessian(const gen & args,GIAC_CONTEXT){ 1041 if ( args.type==_STRNG && args.subtype==-1) return args; 1042 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 1043 return gentypeerr(contextptr); 1044 vecteur v(plotpreprocess(args,contextptr)); 1045 if (is_undef(v)) 1046 return v; 1047 gen g1=v.front(),g2=v.back(); 1048 if (g2.type!=_VECT) 1049 return gentypeerr(contextptr); 1050 vecteur v2(*g2._VECTptr); 1051 int n=int(v2.size()); 1052 vecteur he; 1053 for (int k=0;k<n;k++){ 1054 vecteur l(n); 1055 for (int j=0;j<n;j++){ 1056 l[j]=derive(derive(g1,v2[k],contextptr),v2[j],contextptr); 1057 } 1058 he.push_back(l); 1059 } 1060 return (he); 1061 } 1062 static const char _hessian_s[]="hessian"; 1063 static define_unary_function_eval_quoted (__hessian,&_hessian,_hessian_s); 1064 define_unary_function_ptr5( at_hessian ,alias_at_hessian,&__hessian,_QUOTE_ARGUMENTS,true); 1065 _divergence(const gen & args,GIAC_CONTEXT)1066 gen _divergence(const gen & args,GIAC_CONTEXT){ 1067 if ( args.type==_STRNG && args.subtype==-1) return args; 1068 if (args.type==_VECT && args._VECTptr->size()==3){ 1069 gen opt=args._VECTptr->back(); 1070 if (opt.is_symb_of_sommet(at_equal) && (opt._SYMBptr->feuille[0]==at_coordonnees || (opt._SYMBptr->feuille[0].type==_INT_ && opt._SYMBptr->feuille[0].val==_COORDS))){ 1071 gen g=eval(args._VECTptr->front(),1,contextptr); 1072 gen coord=(*args._VECTptr)[1]; 1073 if (g.type==_VECT && coord.type==_VECT && g._VECTptr->size()==3 && coord._VECTptr->size()==3){ 1074 if (opt._SYMBptr->feuille[1]==at_sphere){ 1075 // spheric: 1/r^2*d_r(r^2*A_r)+1/r/sin(theta)*d_theta(sin(theta)*A_theta)+1/r*sin(theta)*d_phi(A_phi) 1076 gen Ar=g[0],At=g[1],Ap=g[2],r=coord[0],r2=pow(r,2,contextptr),t=coord[1],s=sin(t,contextptr),p=coord[2]; 1077 return derive(r2*Ar,r,contextptr)/r2+(derive(s*At,t,contextptr)+derive(Ap,p,contextptr))/(r*s); 1078 } 1079 if (opt._SYMBptr->feuille[1]==at_cylindre){ 1080 // cylindric: 1/r*d_r(r*A_r)+1/r*d_theta(A_theta)+d_z(A_z) 1081 gen Ar=g[0],At=g[1],Az=g[2],r=coord[0],t=coord[1],z=coord[2]; 1082 return (derive(r*Ar,r,contextptr)+derive(At,t,contextptr))/r+derive(Az,z,contextptr); 1083 } 1084 } 1085 } 1086 } 1087 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 1088 return gentypeerr(contextptr); 1089 vecteur v(plotpreprocess(args,contextptr)); 1090 if (is_undef(v)) 1091 return v; 1092 gen g1=v.front(),g2=v.back(); 1093 if ((g1.type!=_VECT) ||(g2.type!=_VECT)) 1094 return gentypeerr(contextptr); 1095 vecteur v1(*g1._VECTptr); 1096 vecteur v2(*g2._VECTptr); 1097 int n=int(v2.size()); 1098 gen di; 1099 di=0; 1100 for (int k=0;k<n;k++){ 1101 di=di+derive(v1[k],v2[k],contextptr); 1102 } 1103 return normal(di,contextptr); 1104 } 1105 static const char _divergence_s[]="divergence"; 1106 static define_unary_function_eval_quoted (__divergence,&_divergence,_divergence_s); 1107 define_unary_function_ptr5( at_divergence ,alias_at_divergence,&__divergence,_QUOTE_ARGUMENTS,true); 1108 _curl(const gen & args,GIAC_CONTEXT)1109 gen _curl(const gen & args,GIAC_CONTEXT){ 1110 if ( args.type==_STRNG && args.subtype==-1) return args; 1111 if (args.type==_VECT && args._VECTptr->size()==3){ 1112 gen opt=args._VECTptr->back(); 1113 if (opt.is_symb_of_sommet(at_equal) && (opt._SYMBptr->feuille[0]==at_coordonnees || (opt._SYMBptr->feuille[0].type==_INT_ && opt._SYMBptr->feuille[0].val==_COORDS))){ 1114 gen g=eval(args._VECTptr->front(),1,contextptr); 1115 gen coord=(*args._VECTptr)[1]; 1116 if (g.type==_VECT && coord.type==_VECT && g._VECTptr->size()==3 && coord._VECTptr->size()==3){ 1117 if (opt._SYMBptr->feuille[1]==at_sphere){ 1118 // spheric: 1/r/sin(theta)*(d_theta(sin(theta)*A_phi)-d_phi(A_theta)), 1119 // 1/r/sin(theta)*d_phi(A_r)-1/r*d_r(r*A_phi), 1/r*(d_r(r*A_theta)-d_theta(A_r)) 1120 gen Ar=g[0],At=g[1],Ap=g[2],r=coord[0],r2=pow(r,2,contextptr),t=coord[1],s=sin(t,contextptr),p=coord[2]; 1121 return makevecteur((derive(s*Ap,t,contextptr)-derive(At,p,contextptr))/(r*s),derive(Ar,p,contextptr)/(r*s)-derive(r*Ap,r,contextptr)/r,(derive(r*At,r,contextptr)-derive(Ar,t,contextptr))/r); 1122 } 1123 if (opt._SYMBptr->feuille[1]==at_cylindre){ 1124 // cylindric (1/r*d_theta(A_z)-d_z(A_theta)),d_z(A_r)-d_r(A_z),1/r*d_r(r*A_theta)-d_theta(A_r)) 1125 gen Ar=g[0],At=g[1],Az=g[2],r=coord[0],t=coord[1],z=coord[2]; 1126 return makevecteur(derive(Az,t,contextptr)/r-derive(At,z,contextptr),derive(Ar,z,contextptr)-derive(Az,r,contextptr),(derive(r*At,r,contextptr)-derive(Ar,t,contextptr))/r); 1127 } 1128 } 1129 } 1130 } 1131 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 1132 return gentypeerr(contextptr); 1133 vecteur v(plotpreprocess(args,contextptr)); 1134 if (is_undef(v)) 1135 return v; 1136 gen g1=v.front(),g2=v.back(); 1137 if ((g1.type!=_VECT) ||(g2.type!=_VECT)) 1138 return gentypeerr(contextptr); 1139 vecteur v1(*g1._VECTptr); 1140 vecteur v2(*g2._VECTptr); 1141 int n=int(v2.size()); 1142 if (n!=3) return gensizeerr(contextptr); 1143 vecteur rot(3); 1144 rot[0]=derive(v1[2],v2[1],contextptr)-derive(v1[1],v2[2],contextptr); 1145 rot[1]=derive(v1[0],v2[2],contextptr)-derive(v1[2],v2[0],contextptr); 1146 rot[2]=derive(v1[1],v2[0],contextptr)-derive(v1[0],v2[1],contextptr); 1147 return rot; 1148 } 1149 static const char _curl_s[]="curl"; 1150 static define_unary_function_eval_quoted (__curl,&_curl,_curl_s); 1151 define_unary_function_ptr5( at_curl ,alias_at_curl,&__curl,_QUOTE_ARGUMENTS,true); 1152 find_n_x(const gen & args,int & n,gen & x,gen & a)1153 bool find_n_x(const gen & args,int & n,gen & x,gen & a){ 1154 if (args.type==_INT_){ 1155 n=args.val; 1156 x=vx_var; 1157 a=a__IDNT_e; 1158 return n>=0; 1159 } 1160 if (args.type==_DOUBLE_){ 1161 n=int(args._DOUBLE_val); 1162 if (n!=args._DOUBLE_val) 1163 return false; 1164 x=vx_var; 1165 a=a__IDNT_e; 1166 return n>=0; 1167 } 1168 if (args.type!=_VECT || args._VECTptr->size()<2) 1169 return false; 1170 vecteur v=*args._VECTptr; 1171 if (v[0].type==_DOUBLE_ && v[0]._DOUBLE_val==int(v[0]._DOUBLE_val)) 1172 v[0]=int(v[0]._DOUBLE_val); 1173 if (v[0].type==_INT_){ 1174 n=v[0].val; 1175 x=v[1]; 1176 } 1177 else 1178 return false; 1179 if (v.size()>2) 1180 a=v[2]; 1181 else 1182 a=a__IDNT_e; 1183 return n>=0; 1184 } 1185 hermite(int n)1186 vecteur hermite(int n){ 1187 vecteur v(n+1); 1188 v[0]=pow(plus_two,n); 1189 for (int k=2;k<=n;k+=2){ 1190 v[k]=-((n+2-k)*(n+1-k)*v[k-2])/(2*k); 1191 if (is_undef(v[k])) 1192 return v; 1193 } 1194 return v; 1195 } 1196 poly2symbmat(matrice & A,const gen & var,GIAC_CONTEXT)1197 bool poly2symbmat(matrice & A,const gen & var,GIAC_CONTEXT){ 1198 iterateur it=A.begin(),itend=A.end(); 1199 for (;it!=itend;++it){ 1200 if (it->type!=_VECT) return false; 1201 iterateur jt=it->_VECTptr->begin(),jtend=it->_VECTptr->end(); 1202 for (;jt!=jtend;++jt){ 1203 *jt=_poly2symb(makesequence(*jt,var),contextptr); 1204 } 1205 } 1206 return true; 1207 } 1208 unmod(std_matrix<gen> & a,const gen & p)1209 bool unmod(std_matrix<gen> & a,const gen & p){ 1210 for (size_t i=0;i<a.size();++i){ 1211 dbgprint_vector<gen> & v=a[i]; 1212 for (size_t j=0;j<v.size();++j){ 1213 gen & g=v[j]; 1214 if (g.type==_MOD){ 1215 if (*(g._MODptr+1)!=p) 1216 return false; 1217 g=unmod(g); 1218 } 1219 if (g.type==_VECT){ 1220 iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 1221 for (;it!=itend;++it){ 1222 gen & h =*it; 1223 if (h.type==_MOD){ 1224 if (*(h._MODptr+1)!=p) 1225 return false; 1226 h=unmod(h); 1227 } 1228 } 1229 } 1230 } 1231 } 1232 return true; 1233 } 1234 hnf(const matrice & Aorig,const gen & var,matrice & U,matrice & A,matrice & V,bool do_smith,GIAC_CONTEXT)1235 bool hnf(const matrice & Aorig,const gen & var,matrice & U,matrice & A,matrice & V,bool do_smith,GIAC_CONTEXT){ 1236 std_matrix<gen> aorig,u,a,v; 1237 if (var.type==_VECT){ 1238 if (var._VECTptr->empty()) 1239 matrice2std_matrix_gen(Aorig,aorig); 1240 else 1241 return false; // multivariate polynomial are not yet allowed 1242 } 1243 else { 1244 gen tmp=_symb2poly(makesequence(Aorig,var),contextptr); 1245 if (!ckmatrix(tmp)) 1246 return false; 1247 matrice2std_matrix_gen(*tmp._VECTptr,aorig); 1248 } 1249 // fixme, detect modular computation 1250 environment env; 1251 gen g=aorig[0][0]; 1252 if (g.type==_VECT && g._VECTptr->back().type==_MOD){ 1253 env.modulo=*(g._VECTptr->back()._MODptr+1); 1254 env.moduloon=true; 1255 if (!unmod(aorig,env.modulo)) 1256 return false; 1257 } 1258 else 1259 env.moduloon=false; 1260 if (do_smith){ 1261 if (!smith(aorig,u,a,v,&env,contextptr)) 1262 return false; 1263 } 1264 else { 1265 if (!hermite(aorig,u,a,&env,contextptr)) 1266 return false; 1267 } 1268 std_matrix_gen2matrice_destroy(u,U); 1269 std_matrix_gen2matrice_destroy(a,A); 1270 if (do_smith) 1271 std_matrix_gen2matrice_destroy(v,V); 1272 if (var.type!=_VECT || !var._VECTptr->empty()) { 1273 poly2symbmat(U,var,contextptr); 1274 poly2symbmat(A,var,contextptr); 1275 if (do_smith) 1276 poly2symbmat(V,var,contextptr); 1277 } 1278 if (env.moduloon){ 1279 U=*makemod(U,env.modulo)._VECTptr; 1280 A=*makemod(A,env.modulo)._VECTptr; 1281 if (do_smith) 1282 V=*makemod(V,env.modulo)._VECTptr; 1283 } 1284 return true; 1285 } 1286 _hermite(const gen & args,GIAC_CONTEXT)1287 gen _hermite(const gen & args,GIAC_CONTEXT){ 1288 if ( args.type==_STRNG && args.subtype==-1) return args; 1289 int n; 1290 gen a,x; 1291 if (!find_n_x(args,n,x,a)){ 1292 matrice U,A,V; // U invertible and A such that U*args==A 1293 if (ckmatrix(args) && hnf(*args._VECTptr,ggb_var(args),U,A,V,false,contextptr)) 1294 return makesequence(U,A); 1295 if (args.type==_VECT && args._VECTptr->size()==2 && ckmatrix(args._VECTptr->front()) && hnf(*args._VECTptr->front()._VECTptr,args._VECTptr->back(),U,A,V,false,contextptr)) 1296 return makesequence(U,A); 1297 return gensizeerr(contextptr); 1298 } 1299 return r2e(hermite(n),x,contextptr); 1300 } 1301 static const char _hermite_s[]="hermite"; 1302 static define_unary_function_eval (__hermite,&_hermite,_hermite_s); 1303 define_unary_function_ptr5( at_hermite ,alias_at_hermite,&__hermite,0,true); 1304 _smith(const gen & args,GIAC_CONTEXT)1305 gen _smith(const gen & args,GIAC_CONTEXT){ 1306 if ( args.type==_STRNG && args.subtype==-1) return args; 1307 matrice U,A,V; // U invertible and A such that U*args==A 1308 if (ckmatrix(args) && hnf(*args._VECTptr,ggb_var(args),U,A,V,true,contextptr)) 1309 return makesequence(U,A,V); 1310 if (args.type==_VECT && args._VECTptr->size()==2 && ckmatrix(args._VECTptr->front()) && hnf(*args._VECTptr->front()._VECTptr,args._VECTptr->back(),U,A,V,true,contextptr)) 1311 return makesequence(U,A,V); 1312 return gensizeerr(contextptr); 1313 } 1314 static const char _smith_s[]="smith"; 1315 static define_unary_function_eval (__smith,&_smith,_smith_s); 1316 define_unary_function_ptr5( at_smith ,alias_at_smith,&__smith,0,true); 1317 laguerre(int n)1318 vecteur laguerre(int n){ 1319 // L_k(x)=v_k(x)/k! 1320 // L_{n+1}=1/(n+1)*((2n+1-x)*L_n(x)-nL_{n-1}(x)) 1321 // hence v_{n+1}=(2n+1-x)v_n-n*(n-1)*v_{n-1} 1322 vecteur v0,v1,tmp1,tmp2,tmpx; 1323 v0.reserve(n+1); 1324 v1.reserve(n+1); 1325 tmp1.reserve(n+1); 1326 tmp2.reserve(n+1); 1327 tmpx=makevecteur(-1,0); 1328 v0.push_back(1); // v0=1 1329 v1.push_back(-1); 1330 v1.push_back(1); // v1=1-x 1331 for (int k=1;k<n;k++){ 1332 tmpx[1]=2*k+1; 1333 mulmodpoly(tmpx,v1,0,tmp1); 1334 mulmodpoly(v0,gen(k*k),0,tmp2); 1335 submodpoly(tmp1,tmp2,0,tmp1); 1336 // v0, v1, tmp1 -> v1, v0, tmp1 -> v1, tmp1, v0 1337 v0.swap(v1); 1338 v1.swap(tmp1); 1339 if (is_undef(v1)) 1340 return v1; 1341 // cerr << v0 << " " << v1 << '\n'; 1342 } 1343 return v1; 1344 } 1345 _laguerre(const gen & args,GIAC_CONTEXT)1346 gen _laguerre(const gen & args,GIAC_CONTEXT){ 1347 if ( args.type==_STRNG && args.subtype==-1) return args; 1348 int n; 1349 gen a,x; 1350 if (!find_n_x(args,n,x,a)) 1351 return gensizeerr(contextptr); 1352 if (is_zero(a)) 1353 return inv(factorial(n),contextptr)*symb_horner(laguerre(n),x); 1354 gen p0,p1,p2; 1355 p0=1; 1356 p1=1+a-x; 1357 if (n==0) return p0; 1358 if (n==1) return p1; 1359 for (int k=2;k<=n;k++){ 1360 //p2=rdiv(2*k+a-1-x,k)*p1-rdiv(k+a-1,k)*p0; 1361 p2=(2*k+a-1-x)*p1-(k-1)*(k+a-1)*p0; 1362 p0=p1; 1363 p1=p2; 1364 } 1365 //return normal(p2,contextptr); 1366 return normal(rdiv(p2,factorial(n),contextptr),contextptr); 1367 } 1368 1369 static const char _laguerre_s[]="laguerre"; 1370 static define_unary_function_eval (__laguerre,&_laguerre,_laguerre_s); 1371 define_unary_function_ptr5( at_laguerre ,alias_at_laguerre,&__laguerre,0,true); 1372 1373 // Improved one tchebyshev1(int n)1374 vecteur tchebyshev1(int n){ 1375 if (n==0) return vecteur(1,1); 1376 vecteur v(n+1); 1377 v[0]=pow(gen(2),n-1); 1378 if (n==1) return v; 1379 for (int k=2;k<=n;k+=2){ 1380 v[k]=-((n-k+2)*(n-k+1)*v[k-2])/(2*k*(n-k/2)); 1381 if (is_undef(v[k])) 1382 return v; 1383 } 1384 return v; 1385 } tchebyshev_eval(const gen & n,const gen & x,const vecteur & v,GIAC_CONTEXT)1386 gen tchebyshev_eval(const gen & n,const gen &x,const vecteur & v,GIAC_CONTEXT){ 1387 if (is_zero(n)) 1388 return 1; 1389 if (!is_positive(n,contextptr)) 1390 return gensizeerr(contextptr); 1391 #if 1 1392 vecteur w(v); 1393 modpoly X(makevecteur(1,0)); 1394 modpoly B(makevecteur(1,-2*x,1)); 1395 environment env; 1396 if (x.type==_MOD){ 1397 env.moduloon=true; 1398 env.modulo=*(x._MODptr+1); 1399 B[1]=*B[1]._MODptr; 1400 if (v[1].type==_MOD) 1401 w[1]=*v[1]._MODptr; 1402 } 1403 else 1404 env.moduloon=false; 1405 modpoly res=powmod(X,n,B,&env); 1406 if (res.empty()) 1407 return gensizeerr(contextptr); 1408 if (res.size()==1) 1409 res[0]=res[0]*w[0]; 1410 else { 1411 matrice A(makevecteur(makevecteur(res[1],res[0]),makevecteur(-res[0],2*res[0]*x+res[1]))); 1412 multmatvecteur(A,w,res); 1413 } 1414 if (env.moduloon && res[0].type!=_MOD) 1415 res[0]=makemod(res[0],env.modulo); 1416 return res[0]; 1417 #else 1418 matrice A(makevecteur(makevecteur(0,1),makevecteur(-1,2*x))); 1419 gen An=pow(A,n,contextptr); 1420 if (An.type!=_VECT) 1421 return undef; 1422 An=ckmultmatvecteur(*An._VECTptr,v,context0); 1423 return An[0]; 1424 #endif 1425 } _tchebyshev1(const gen & args,GIAC_CONTEXT)1426 gen _tchebyshev1(const gen & args,GIAC_CONTEXT){ 1427 if ( args.type==_STRNG && args.subtype==-1) return args; 1428 int n; 1429 gen a,x; 1430 if (args.type==_VECT && args._VECTptr->size()==2 && is_integer(args._VECTptr->front())) 1431 return tchebyshev_eval(args._VECTptr->front(),args._VECTptr->back(),makevecteur(1,args._VECTptr->back()),contextptr); 1432 if (!find_n_x(args,n,x,a)) 1433 return gensizeerr(contextptr); 1434 if (n==0) 1435 return 1; 1436 return r2e(tchebyshev1(n),x,contextptr); 1437 } 1438 static const char _tchebyshev1_s[]="tchebyshev1"; 1439 static define_unary_function_eval (__tchebyshev1,&_tchebyshev1,_tchebyshev1_s); 1440 define_unary_function_ptr5( at_tchebyshev1 ,alias_at_tchebyshev1,&__tchebyshev1,0,true); 1441 1442 // Improved one tchebyshev2(int n)1443 vecteur tchebyshev2(int n){ 1444 vecteur v(n+1); 1445 v[0]=pow(gen(2),n); 1446 for (int k=1;k<=n/2;++k){ 1447 v[2*k]=-(n+2-2*k)*(n+1-2*k)*v[2*k-2]/(4*k*(n+1-k)); 1448 if (is_undef(v[2*k])) 1449 return v; 1450 } 1451 return v; 1452 } _tchebyshev2(const gen & args,GIAC_CONTEXT)1453 gen _tchebyshev2(const gen & args,GIAC_CONTEXT){ 1454 if ( args.type==_STRNG && args.subtype==-1) return args; 1455 if (args.type==_VECT && args._VECTptr->size()==2 && is_integer(args._VECTptr->front())) 1456 return tchebyshev_eval(args._VECTptr->front(),args._VECTptr->back(),makevecteur(0,1),contextptr); 1457 int n; 1458 gen a,x; 1459 if (!find_n_x(args,n,x,a)) 1460 return gensizeerr(contextptr); 1461 return r2e(tchebyshev2(n),x,contextptr); 1462 } 1463 static const char _tchebyshev2_s[]="tchebyshev2"; 1464 static define_unary_function_eval (__tchebyshev2,&_tchebyshev2,_tchebyshev2_s); 1465 define_unary_function_ptr5( at_tchebyshev2 ,alias_at_tchebyshev2,&__tchebyshev2,0,true); 1466 1467 // Legendre is the vecteur returned divided by n! legendre(int n)1468 vecteur legendre(int n){ 1469 vecteur v0,v1,vtmp1,vtmp2; 1470 v0.push_back(1); 1471 v1.push_back(1); 1472 v1.push_back(0); 1473 if (!n) return v0; 1474 if (n==1) return v1; 1475 for (int k=2;k<=n;k++){ 1476 multvecteur(2*k-1,v1,vtmp1); 1477 vtmp1.push_back(0); // (2k-1)*x*p1 1478 multvecteur((k-1)*(k-1),v0,vtmp2); // (k-1)^2*p0 1479 vtmp1=vtmp1-vtmp2; // p2=(2*k-1)*x*p1-(k-1)*(k-1)*p0; 1480 v0=v1; 1481 v1=vtmp1; 1482 } 1483 return v1; 1484 } 1485 _legendre(const gen & args,GIAC_CONTEXT)1486 gen _legendre(const gen & args,GIAC_CONTEXT){ 1487 if ( args.type==_STRNG && args.subtype==-1) return args; 1488 int n; 1489 gen a,x; 1490 if (!find_n_x(args,n,x,a)) 1491 return gensizeerr(contextptr); 1492 vecteur v=multvecteur(inv(factorial(n),contextptr),legendre(n)); 1493 return r2e(v,x,contextptr); 1494 } 1495 static const char _legendre_s[]="legendre"; 1496 static define_unary_function_eval (__legendre,&_legendre,_legendre_s); 1497 define_unary_function_ptr5( at_legendre ,alias_at_legendre,&__legendre,0,true); 1498 1499 // arithmetic mean column by column mean(const matrice & m,bool column)1500 vecteur mean(const matrice & m,bool column){ 1501 matrice mt; 1502 if (column) 1503 mt=mtran(m); 1504 else 1505 mt=m; 1506 vecteur res; 1507 const_iterateur it=mt.begin(),itend=mt.end(); 1508 for (;it!=itend;++it){ 1509 const gen & g =*it; 1510 if (g.type!=_VECT){ 1511 res.push_back(g); 1512 continue; 1513 } 1514 vecteur & v=*g._VECTptr; 1515 if (v.empty()){ 1516 res.push_back(undef); 1517 continue; 1518 } 1519 const_iterateur jt=v.begin(),jtend=v.end(); 1520 int s=int(jtend-jt); 1521 gen somme(0); 1522 for (;jt!=jtend;++jt){ 1523 //somme = somme + evalf(*jt); 1524 somme = somme + *jt; 1525 } 1526 res.push_back(rdiv(somme,s,context0)); 1527 } 1528 return res; 1529 } 1530 stddev(const matrice & m,bool column,int variance)1531 vecteur stddev(const matrice & m,bool column,int variance){ 1532 matrice mt; 1533 if (column) 1534 mt=mtran(m); 1535 else 1536 mt=m; 1537 vecteur moyenne(mean(mt,false)); 1538 vecteur res; 1539 const_iterateur it=mt.begin(),itend=mt.end(); 1540 for (int i=0;it!=itend;++it,++i){ 1541 const gen & g =*it; 1542 if (g.type!=_VECT){ 1543 res.push_back(0); 1544 continue; 1545 } 1546 vecteur & v=*g._VECTptr; 1547 if (v.empty()){ 1548 res.push_back(undef); 1549 continue; 1550 } 1551 const_iterateur jt=v.begin(),jtend=v.end(); 1552 int s=int(jtend-jt); 1553 gen somme(0); 1554 for (;jt!=jtend;++jt){ 1555 // somme = somme + evalf((*jt)*(*jt)); 1556 somme = somme + (*jt)*(*jt); 1557 } 1558 if (variance!=3) 1559 res.push_back(sqrt(rdiv(somme-s*moyenne[i]*moyenne[i],s-(variance==2),context0),context0)); 1560 else 1561 res.push_back(rdiv(somme,s,context0)-moyenne[i]*moyenne[i]); 1562 } 1563 return res; 1564 } 1565 ascsort(const matrice & m,bool column)1566 matrice ascsort(const matrice & m,bool column){ 1567 matrice mt; 1568 if (column) 1569 mt=mtran(m); 1570 else 1571 mt=m; 1572 iterateur it=mt.begin(),itend=mt.end(); 1573 gen tmp; 1574 for (;it!=itend;++it){ 1575 gen & g =*it; 1576 if (g.type!=_VECT) 1577 continue; 1578 vecteur v=*g._VECTptr; 1579 if (v.empty()) 1580 continue; 1581 const_iterateur jt=v.begin(),jtend=v.end(); 1582 int n=int(jtend-jt); 1583 vector<double> vv(n); 1584 for (int j=0;jt!=jtend;++jt,++j){ 1585 if ( (jt->type==_VECT) && (jt->_VECTptr->size()==3) ) 1586 tmp=(*jt->_VECTptr)[1]; 1587 else 1588 tmp=*jt; 1589 tmp=evalf(tmp,1,0); 1590 if (tmp.type!=_DOUBLE_) 1591 vv[j]=0; 1592 else 1593 vv[j]=tmp._DOUBLE_val; 1594 } 1595 sort(vv.begin(),vv.end()); 1596 for (int j=0;j<n;++j) 1597 v[j]=vv[j]; 1598 *it=v; 1599 } 1600 return mt; 1601 } 1602 _permu2mat(const gen & args,GIAC_CONTEXT)1603 gen _permu2mat(const gen & args,GIAC_CONTEXT){ 1604 if ( args.type==_STRNG && args.subtype==-1) return args; 1605 //transforme une permutation en une matrice obtenue en permutant les lignes de la matrice identite 1606 if (args.type!=_VECT) 1607 return gentypeerr(contextptr); 1608 vector<int> p1; 1609 vecteur p(*args._VECTptr); 1610 if (!(is_permu(p,p1,contextptr))) 1611 return gentypeerr(contextptr); 1612 int n=int(p.size()); 1613 vecteur c; 1614 vecteur l(n); 1615 for (int k=0;k<n;k++){ 1616 for (int j=0;j<n;j++){ 1617 if (p[k]==j+array_start(contextptr)) { 1618 l[j]=1; 1619 } else { 1620 l[j]=0; 1621 } 1622 } 1623 c.push_back(l); 1624 } 1625 return c; 1626 } 1627 1628 static const char _permu2mat_s[]="permu2mat"; 1629 static define_unary_function_eval (__permu2mat,&_permu2mat,_permu2mat_s); 1630 define_unary_function_ptr5( at_permu2mat ,alias_at_permu2mat,&__permu2mat,0,true); 1631 est_dans(const vector<int> & a,const int n,vector<vector<int>> s)1632 static bool est_dans(const vector<int> & a , const int n, vector< vector<int> > s) { 1633 //teste si a est egal a l'un des n premiers elements de s 1634 bool cont=true; 1635 int j=0; 1636 while (j<=n && cont) { 1637 if (a==s[j]) { 1638 cont=false; 1639 } 1640 j=j+1; 1641 } 1642 return (! cont); 1643 } 1644 groupermu(const vector<int> & p1,const vector<int> & p2)1645 vector< vector<int> > groupermu(const vector<int> & p1,const vector<int> & p2) { 1646 //groupe engendre par de 2 perm p1 et p2 de long n1 et n2 : a pour long max(n1,n2) 1647 int n1; 1648 n1=int(p1.size()); 1649 int n2; 1650 n2=int(p2.size()); 1651 vector<int> a; 1652 vector<int> b; 1653 a=p1; 1654 b=p2; 1655 if (n1>n2) { 1656 for (int k=n2;k<n1;k++) {b.push_back(k);n2=n1;} 1657 } else { 1658 for (int k=n1;k<n2;k++) {a.push_back(k);n1=n2;} 1659 } 1660 ; 1661 vector< vector<int> > s(2); 1662 s[0]=a; 1663 s[1]=b; 1664 int p=0; 1665 int q=1; 1666 bool pfini=true; 1667 while (pfini) { 1668 int k=0; 1669 for (int j=p;j<=q;j++){ 1670 vector<int> na; 1671 vector<int> nb; 1672 na=p1op2(a,s[j]); 1673 if (!(est_dans(na,q+k,s))){ 1674 k=k+1; 1675 s.push_back(na); 1676 } 1677 nb=p1op2(b,s[j]); 1678 if (!(est_dans(nb,q+k,s))){ 1679 k=k+1; 1680 s.push_back(nb); 1681 } 1682 } 1683 if (k!=0) { 1684 p=q+1; 1685 q=q+k; 1686 } else { 1687 pfini=false; 1688 } 1689 } 1690 //q est l'ordre du groupe = size(s) 1691 return(s); 1692 } 1693 _groupermu(const gen & args,GIAC_CONTEXT)1694 gen _groupermu(const gen & args,GIAC_CONTEXT){ 1695 if ( args.type==_STRNG && args.subtype==-1) return args; 1696 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 1697 return gentypeerr(contextptr); 1698 vecteur v(*args._VECTptr); 1699 gen v1=v.front(),v2=v.back(); 1700 if ( (v1.type!=_VECT) || (v2.type!=_VECT)) 1701 return gentypeerr(contextptr); 1702 vector<int> p1,p2; 1703 if (!is_permu(*v1._VECTptr,p1,contextptr) || !is_permu(*v2._VECTptr,p2,contextptr)) 1704 return gensizeerr(contextptr); 1705 return vectvector_int_2_vecteur(groupermu(p1,p2),contextptr); 1706 } 1707 static const char _groupermu_s[]="groupermu"; 1708 static define_unary_function_eval (__groupermu,&_groupermu,_groupermu_s); 1709 define_unary_function_ptr5( at_groupermu ,alias_at_groupermu,&__groupermu,0,true); 1710 _nextperm(const gen & args,GIAC_CONTEXT)1711 gen _nextperm(const gen & args,GIAC_CONTEXT){ 1712 if ( args.type==_STRNG && args.subtype==-1) return args; 1713 if (args.type!=_VECT) 1714 return gentypeerr(contextptr); 1715 const vecteur & v1=*args._VECTptr; 1716 vector<int> p1; 1717 if (!is_permu(v1,p1,contextptr)) 1718 return gensizeerr(contextptr); 1719 if (next_permutation(p1.begin(),p1.end())) 1720 return vector_int_2_vecteur(p1,contextptr); 1721 else 1722 return undef; 1723 } 1724 static const char _nextperm_s[]="nextperm"; 1725 static define_unary_function_eval (__nextperm,&_nextperm,_nextperm_s); 1726 define_unary_function_ptr5( at_nextperm ,alias_at_nextperm,&__nextperm,0,true); 1727 _prevperm(const gen & args,GIAC_CONTEXT)1728 gen _prevperm(const gen & args,GIAC_CONTEXT){ 1729 if ( args.type==_STRNG && args.subtype==-1) return args; 1730 if (args.type!=_VECT) 1731 return gentypeerr(contextptr); 1732 const vecteur & v1=*args._VECTptr; 1733 vector<int> p1; 1734 if (!is_permu(v1,p1,contextptr)) 1735 return gensizeerr(contextptr); 1736 if (prev_permutation(p1.begin(),p1.end())) 1737 return vector_int_2_vecteur(p1,contextptr); 1738 else 1739 return undef; 1740 } 1741 static const char _prevperm_s[]="prevperm"; 1742 static define_unary_function_eval (__prevperm,&_prevperm,_prevperm_s); 1743 define_unary_function_ptr5( at_prevperm ,alias_at_prevperm,&__prevperm,0,true); 1744 _split(const gen & args,GIAC_CONTEXT)1745 gen _split(const gen & args,GIAC_CONTEXT){ 1746 if ( args.type==_STRNG && args.subtype==-1) return args; 1747 //renvoie [ax,ay] si les arg sont g1=ax*an (sans denominateur) et g2=[x,y] 1748 //sinon renvoie [0] 1749 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 1750 return gentypeerr(contextptr); 1751 vecteur v(*args._VECTptr); 1752 gen g1=v.front(),g2=v.back(); 1753 if (g1.type==_STRNG && g2.type==_STRNG){ // Python "Ceci est une phrase".split(" ") 1754 if (g2._STRNGptr->empty()) 1755 return gendimerr(contextptr); 1756 vecteur res; 1757 int pos=0,ss=g1._STRNGptr->size(); 1758 for (;pos<ss;){ 1759 int npos=g1._STRNGptr->find(*g2._STRNGptr,pos); 1760 if (npos<0 || npos>=ss) 1761 break; 1762 res.push_back(string2gen(g1._STRNGptr->substr(pos,npos-pos),false)); 1763 pos=npos+g2._STRNGptr->size(); 1764 } 1765 if (pos<ss) 1766 res.push_back(string2gen(g1._STRNGptr->substr(pos,ss-pos),false)); 1767 return res; 1768 } 1769 if (g2.type!=_VECT) 1770 return gentypeerr(contextptr); 1771 vecteur v2(*g2._VECTptr); 1772 int n=int(v2.size()); 1773 if (n!=2) return gensizeerr(contextptr); 1774 vecteur fa; 1775 fa=factors(g1,vx_var,contextptr); 1776 int l=int(fa.size()); 1777 gen ax=1; 1778 gen ay=1; 1779 for (int k=0;k<l;k=k+2){ 1780 gen f=fa[k]; 1781 if (derive(f,v2[0],contextptr)==0) { 1782 ay=ay*pow(f,fa[k+1],contextptr); 1783 } 1784 else { 1785 if (derive(f,v2[1],contextptr)==0){ 1786 ax=ax*pow(f,fa[k+1],contextptr); 1787 } 1788 else {vecteur res(1);return (res);} 1789 } 1790 } 1791 vecteur res(2); 1792 res[0]=ax; 1793 res[1]=ay; 1794 return res; 1795 } 1796 1797 static const char _split_s[]="split"; 1798 static define_unary_function_eval (__split,&_split,_split_s); 1799 define_unary_function_ptr5( at_split ,alias_at_split,&__split,0,true); 1800 _join(const gen & args,GIAC_CONTEXT)1801 gen _join(const gen & args,GIAC_CONTEXT){ 1802 if ( args.type==_STRNG && args.subtype==-1) return args; 1803 if (args.type==_VECT && args._VECTptr->size()==2){ 1804 gen g1=args._VECTptr->front(),g2=args._VECTptr->back(); 1805 if (g1.type==_STRNG && g2.type==_VECT){ // Python 1806 const_iterateur it=g2._VECTptr->begin(),itend=g2._VECTptr->end(); 1807 string res; 1808 for (;it!=itend;){ 1809 if (it->type==_STRNG) 1810 res += *it->_STRNGptr; 1811 else 1812 res += it->print(contextptr); 1813 ++it; 1814 if (it==itend) 1815 break; 1816 res += *g1._STRNGptr; 1817 } 1818 return string2gen(res,false); 1819 } 1820 } 1821 return gensizeerr(contextptr); 1822 } 1823 static const char _join_s[]="join"; 1824 static define_unary_function_eval (__join,&_join,_join_s); 1825 define_unary_function_ptr5( at_join ,alias_at_join,&__join,0,true); 1826 _sum_riemann(const gen & args,GIAC_CONTEXT)1827 gen _sum_riemann(const gen & args,GIAC_CONTEXT){ 1828 if ( args.type==_STRNG && args.subtype==-1) return args; 1829 //renvoie l'equivalent pour n=infini de sum de k=1 a n de g1 avec g2=[n,k]) 1830 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 1831 return gentypeerr(contextptr); 1832 gen g1=args[0],g2=args[1]; 1833 if (g2.type!=_VECT) 1834 return gentypeerr(contextptr); 1835 vecteur v2(*g2._VECTptr); 1836 int sv2=int(v2.size()); 1837 if (sv2!=2) return gensizeerr(contextptr); 1838 identificateur x("rieman_sum_x"); 1839 //on pose k=n*x 1840 //mettre que v2[0]=n et v2[1]=k sont ds N 1841 //_assume(makevecteur(is_plus(v2[0]))); 1842 //_assume(makevecteur(is_plus(v2[1]))); 1843 gen a=recursive_normal(v2[0]*subst(g1,v2[1],x*v2[0],false,contextptr),contextptr); 1844 gen nda=_fxnd(a,contextptr); 1845 vecteur nada(*nda._VECTptr); 1846 //na numerateur de a et da denominateur de a ou k=n*x 1847 gen na=nada[0]; 1848 gen da=nada[1]; 1849 vecteur var(2); 1850 var[0]=na; 1851 v2[1]=x; 1852 var[1]=v2; 1853 //var=[na,[n,x]] 1854 gen b0=_split(var,contextptr); 1855 if (is_undef(b0)) return b0; 1856 //on separe les variables du numerateur 1857 vecteur vb0(*b0._VECTptr); 1858 if ((vb0.size()==1)&& (vb0[0]==0)) 1859 return string2gen(gettext("Probably not a Riemann sum"),false); 1860 var[0]=da; 1861 //var=[da,[n,x]] 1862 gen b1=_split(var,contextptr); 1863 if (is_undef(b1)) return b1; 1864 //on separe les variables du denominateur 1865 vecteur vb1(*b1._VECTptr); 1866 if ((vb1.size()==1)&& (vb1[0]==0)) 1867 return string2gen(gettext("Probably not a Riemann sum"),false); 1868 gen an=vb0[0]/vb1[0]; 1869 gen ax=vb0[1]/vb1[1]; 1870 gen tmp=_integrate(makesequence(ax,x,0,1),contextptr); 1871 if (is_undef(tmp)) return tmp; 1872 gen tmp2=_limit(makesequence(an,v2[0],plus_inf),contextptr); 1873 //tmp n'a pas ete calcule sum_riemann(pi/(2*n)*log(sin(pi*k/(2*n))),[n,k]) 1874 //if ((tmp.type==_SYMB)&&(tmp._SYMBptr->sommet==at_integrate)) 1875 //return _limit(makevecteur(tmp*an,v2[0],plus_inf)); 1876 //tmp ne doit pas etre infini qd tmp2 est nul et reciproquement 1877 if (is_inf(tmp)&& is_zero(tmp2)) 1878 return string2gen(gettext("Probably not a Riemann sum"),false); 1879 if (is_zero(tmp)&& is_inf(tmp2)) 1880 return string2gen(gettext("Probably not a Riemann sum"),false); 1881 gen res=recursive_normal(tmp*_series(makesequence(an,v2[0],plus_inf),contextptr),contextptr); 1882 res=_limit(makesequence(res,v2[0],plus_inf),contextptr); 1883 return res; 1884 } 1885 static const char _sum_riemann_s[]="sum_riemann"; 1886 static define_unary_function_eval (__sum_riemann,&_sum_riemann,_sum_riemann_s); 1887 define_unary_function_ptr5( at_sum_riemann,alias_at_sum_rieman,&__sum_riemann,0,true); 1888 1889 #ifndef NO_NAMESPACE_GIAC 1890 } // namespace giac 1891 #endif // ndef NO_NAMESPACE_GIAC 1892