1 // -*- mode:C++ ; compile-command: "g++-3.4 -I.. -g -c gauss.cc -Wall" -*- 2 #include "giacPCH.h" 3 4 /* 5 * Copyright (C) 2001,14 R. De Graeve, Institut Fourier, 38402 St Martin d'Heres 6 * 7 * This program is free software; you can redistribute it and/or modify 8 * it under the terms of the GNU General Public License as published by 9 * the Free Software Foundation; either version 3 of the License, or 10 * (at your option) any later version. 11 * 12 * This program is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 * GNU General Public License for more details. 16 * 17 * You should have received a copy of the GNU General Public License 18 * along with this program. If not, see <http://www.gnu.org/licenses/>. 19 */ 20 21 using namespace std; 22 //#include <fstream> 23 #include "gauss.h" 24 #include "vecteur.h" 25 #include "derive.h" 26 #include "subst.h" 27 #include "usual.h" 28 #include "sym2poly.h" 29 #include "solve.h" 30 #include "ti89.h" 31 #include "plot.h" 32 #include "misc.h" 33 #include "ifactor.h" 34 #include "prog.h" 35 #include "giacintl.h" 36 37 #ifndef NO_NAMESPACE_GIAC 38 namespace giac { 39 #endif // ndef NO_NAMESPACE_GIAC 40 quad(int & b,const gen & q,const vecteur & x,GIAC_CONTEXT)41 vecteur quad(int &b,const gen & q, const vecteur & x,GIAC_CONTEXT){ 42 //x=vecteur des variables, q=la fonction a tester,n=la dimension de x 43 //b=2 si q est quadratique,=0,1 ou 3 si il y des termes d'ordre 0,1 ou 3 44 //renvoie (la jacobienne de q)/2 45 gen qs; 46 gen dq; 47 gen dqs; 48 gen qdd; 49 int n=int(x.size()); 50 51 vecteur A; 52 //creation d'une matrice carree A d'ordre n 53 for (int i=0;i<n;i++){ 54 vecteur li(n); 55 A.push_back(li); 56 } 57 //A est un vecteur de vecteur=une matrice! 58 //on met ds A :(la jacobienne de q)/2 59 for (int i=0;i<n;i++){ 60 gen qxi(derive(q,x[i],contextptr)); 61 for (int j=i;j<n;j++){ 62 qdd=derive(qxi,x[j],contextptr); 63 qdd=recursive_normal(qdd,contextptr); 64 //cout<<i<<","<<j<<qdd<<'\n'; 65 (*A[j]._VECTptr)[i]=(*A[i]._VECTptr)[j]=rdiv(qdd,2,contextptr); 66 } 67 } 68 //2*A=jacobienne de q 69 //on calcule qs=q en zero 70 //cout<<A<<'\n'; 71 qs=subst(q,x,vecteur(n),false,contextptr); 72 //qs=la valeur de q en 0 73 if (qs !=0){ 74 b=0; 75 return(A); 76 } 77 //on regarde si il y des termes lineaires 78 for (int j=0;j<n;j++){ 79 dq=derive(q,x[j],contextptr); 80 dqs=subst(dq,x,vecteur(n),false,contextptr); 81 //dqs=la diff de q en zero 82 if (dqs!=0){ 83 b=1; 84 return(A); 85 } 86 } 87 for (int i=0;i<n;i++){ 88 for (int j=i;j<n;j++){ 89 for (int k=0;k<n;k++){ 90 if (derive(A[i][j],x[k],contextptr)!=0){ 91 b=3; 92 return(A); 93 } 94 } 95 } 96 } 97 b=2; 98 //(*(A[1]._VECTptr))[0]=21; 99 return(A); 100 } 101 qxa(const gen & q,const vecteur & x,GIAC_CONTEXT)102 vecteur qxa(const gen &q,const vecteur & x,GIAC_CONTEXT){ 103 //transforme une forme quadratique en une matrice symetrique A 104 //(les variables sont dans x) 105 int d; 106 //d nbre de variables 107 d=int(x.size()); 108 // int da; 109 //il faut verifier que q est quadratique 110 vecteur A; 111 int b; 112 A=quad(b,q,x,contextptr); 113 if (b==2) { 114 return(A); 115 } 116 else { 117 return vecteur(1,gensizeerr(gettext("q is not quadratic"))); 118 } 119 return 0; 120 } 121 symb_q2a(const gen & args)122 static gen symb_q2a(const gen & args){ 123 return symbolic(at_q2a,args); 124 } _q2a(const gen & args,GIAC_CONTEXT)125 gen _q2a(const gen & args,GIAC_CONTEXT){ 126 if ( args.type==_STRNG && args.subtype==-1) return args; 127 if (args.type!=_VECT) 128 return _q2a(makesequence(args,lidnt(args)),contextptr); 129 int s=int(args._VECTptr->size()); 130 if (s!=2) 131 return gendimerr(contextptr); 132 if (args._VECTptr->back().type==_VECT) 133 return qxa(args._VECTptr->front(),*args._VECTptr->back()._VECTptr,contextptr); 134 return symb_q2a(args); 135 } 136 static const char _q2a_s []="q2a"; 137 static define_unary_function_eval (__q2a,&_q2a,_q2a_s); 138 define_unary_function_ptr5( at_q2a ,alias_at_q2a,&__q2a,0,true); 139 gauss(const gen & q,const vecteur & x,vecteur & D,vecteur & U,vecteur & P,GIAC_CONTEXT)140 vecteur gauss(const gen & q, const vecteur & x, vecteur & D, vecteur & U, vecteur & P,GIAC_CONTEXT){ 141 int n=int(x.size()),b; 142 gen u1,u2,q1,l1,l2; 143 vecteur R(1); 144 vecteur PR; 145 for (int i=0;i<n-1;i++){ 146 PR.push_back(vecteur(n-1)); 147 } 148 vecteur PP; 149 for (int i=0;i<n;i++){ 150 PP.push_back(vecteur(n)); 151 } 152 vecteur I; 153 if (n) I=midn(n); 154 vecteur L; 155 156 //si q n'est pas quadratique b<>2 et on retourne q 157 vecteur A(quad(b,q,x,contextptr)); 158 if (b!=2){ 159 R[0]=q; 160 D.clear(); 161 U.clear(); 162 return R; 163 } 164 //la forme q est quadratique de matrice A 165 if (q==0) { 166 //R[0]=q; 167 vecteur vide(n); 168 D=vide; 169 U=vide; 170 P=I; 171 return vide; 172 } 173 if (n==1){ 174 gen q0=_factor(q,contextptr); 175 R[0]=q0; 176 vecteur un(1); 177 un[0]=A[0][0]; 178 D=un; 179 U=x; 180 P=I; 181 return(R); 182 } 183 int r; 184 r=n; 185 for (int i=n-1 ;i>=0;i--){ 186 if (A[i][i]!=0) { 187 r=i; 188 } 189 } 190 if (r!=n) { 191 //il y a des termes carres 192 u1=recursive_normal(rdiv(derive(q,x[r],contextptr),plus_two,contextptr),contextptr); 193 q1=recursive_normal(q-rdiv(u1*u1,A[r][r],contextptr),contextptr); 194 vecteur y; 195 //y contient les variables qui restent (on enleve x[r]) 196 for (int j=0;j<n;j++){ 197 if (j!=r){ 198 y.push_back(x[j]); 199 } 200 } 201 L=gauss(q1,y,D,U,PR,contextptr); 202 //on rajoute 1/a_r_r sur la diagonale D 203 R[0]=rdiv(1,A[r][r],contextptr); 204 D=mergevecteur(R,D); 205 //on rajoute u1 aux vecteurs constitue des formes lineaires 206 //q= 1/a_r_r*(u1)^2+... 207 R[0]=u1; 208 U=mergevecteur(R,U); 209 //on complete la matrice PR de dim n-1 en la matrice PP de dim n 210 //1iere ligne les coeff de u1 et rieme colonne doit avoir des 0 211 for (int i=0;i<n;i++){ 212 (*PP[0]._VECTptr)[i]=recursive_normal(derive(u1,x[i],contextptr),contextptr); 213 } 214 for (int i=1;i<n;i++){ 215 for (int j=0;j<r;j++){ 216 (*PP[i]._VECTptr)[j]=PR[i-1][j]; 217 } 218 for (int j=r+1;j<n;j++){ 219 (*PP[i]._VECTptr)[j]=PR[i-1][j-1]; 220 } 221 } 222 P=PP; 223 R[0]=rdiv(pow(u1,2),A[r][r],contextptr); 224 return(mergevecteur(R,L)); 225 } 226 //il n'y a pas de carres 227 int r1; 228 int r2; 229 r1=0; 230 r2=0; 231 for (int i=n-2;i>=0;i--){ 232 for (int j=i+1;j>=1;j--){ 233 if (A[i][j]!=0) { 234 r1=i; 235 r2=j; 236 } 237 } 238 } 239 l1=rdiv(derive(q,x[r1],contextptr),2,contextptr); 240 l2=rdiv(derive(q,x[r2],contextptr),2,contextptr); 241 u1=recursive_normal(l1+l2,contextptr); 242 u2=recursive_normal(l1-l2,contextptr); 243 q1=recursive_normal(q-rdiv(plus_two*l1*l2,A[r1][r2],contextptr),contextptr); 244 vecteur y; 245 for (int j=0;j<n;j++){ 246 if ((j!=r1) && (j!=r2)) { 247 y.push_back(x[j]); 248 } 249 } 250 L=gauss(q1,y,D,U,PR,contextptr); 251 //on rajoute 1/a_r1_r2 et -1/a_r1_r2 sur la diagonale D 252 R[0]=rdiv(1,plus_two*A[r1][r2],contextptr); 253 R.push_back(-R[0]); 254 D=mergevecteur(R,D); 255 //on rajoute u1 et u2 au vecteur U constitue des formes lineaires 256 //q= 1/a_r1_r2*(u1)^2 - 1/a_r1_r2*(u2)^2 + ... 257 R[0]=u1; 258 R[1]=u2; 259 U=mergevecteur(R,U); 260 //on complete la matrice PR de dim n-2 en la matrice PP de dim n 261 //1iere et 2ieme ligne les coeff de u1 et de u2 262 //r1ieme et r2ieme colonne doit avoir des 0 263 for (int i=0;i<n;i++){ 264 (*PP[0]._VECTptr)[i]=recursive_normal(derive(u1,x[i],contextptr),contextptr); 265 (*PP[1]._VECTptr)[i]=recursive_normal(derive(u2,x[i],contextptr),contextptr); 266 } 267 for (int i=2;i<n;i++){ 268 for (int j=0;j<r1;j++){ 269 (*PP[i]._VECTptr)[j]=PR[i-2][j]; 270 } 271 for (int j=r1+1;j<r2;j++){ 272 (*PP[i]._VECTptr)[j]=PR[i-2][j-1]; 273 } 274 for (int j=r2+1;j<n;j++){ 275 (*PP[i]._VECTptr)[j]=PR[i-2][j-2]; 276 } 277 } 278 P=PP; 279 R[0]=rdiv(pow(u1,2),plus_two*A[r1][r2],contextptr); 280 R[1]=rdiv(-pow(u2,2),plus_two*A[r1][r2],contextptr); 281 return(mergevecteur(R,L)); 282 } 283 gauss(const gen & q,const vecteur & x,GIAC_CONTEXT)284 vecteur gauss(const gen & q,const vecteur & x,GIAC_CONTEXT){ 285 vecteur D,U,P; 286 return gauss(q,x,D,U,P,contextptr); 287 } 288 _gauss(const gen & args,GIAC_CONTEXT)289 gen _gauss(const gen & args,GIAC_CONTEXT){ 290 if ( args.type==_STRNG && args.subtype==-1) return args; 291 if (args.type!=_VECT) 292 return _gauss(makesequence(args,lidnt(args)),contextptr); 293 int s=int(args._VECTptr->size()); 294 if (s<2) 295 return gendimerr(contextptr); 296 const gen & arg1=(*args._VECTptr)[1]; 297 if (arg1.type==_VECT){ 298 const vecteur & v=*arg1._VECTptr; 299 vecteur D,U,P; 300 gen w=gauss(args._VECTptr->front(),v.empty()?lidnt(args):v,D,U,P,contextptr); 301 w=_plus(w,contextptr); 302 if (s>2|| v.empty()) 303 return makesequence(w,D,P); 304 return w; 305 } 306 return _randNorm(args,contextptr); 307 } 308 static const char _gauss_s []="gauss"; 309 static define_unary_function_eval (__gauss,&_gauss,_gauss_s); 310 define_unary_function_ptr5( at_gauss ,alias_at_gauss,&__gauss,0,true); 311 axq(const vecteur & A,const vecteur & x,GIAC_CONTEXT)312 gen axq(const vecteur &A,const vecteur & x,GIAC_CONTEXT){ 313 //transforme une matrice carree (symetrique) en la forme quadratique q 314 //(les variables sont dans x) 315 //d nbre de variables 316 //il faut verifier que A est carree 317 //A n'est pas forcement symetrique 318 int d=int(x.size()); 319 int da=int(A.size()); 320 if (!(is_squarematrix(A)) || (da!=d) ) 321 return gensizeerr(gettext("Invalid dimension")); 322 vecteur Ax; 323 multmatvecteur(A,x,Ax); 324 return normal(dotvecteur(x,Ax),contextptr); 325 } 326 symb_a2q(const gen & args)327 static gen symb_a2q(const gen & args){ 328 return symbolic(at_a2q,args); 329 } _a2q(const gen & args,GIAC_CONTEXT)330 gen _a2q(const gen & args,GIAC_CONTEXT){ 331 if ( args.type==_STRNG && args.subtype==-1) return args; 332 if (args.type!=_VECT) 333 return symb_a2q(args); 334 int s=int(args._VECTptr->size()); 335 if (s!=2) 336 return gendimerr(contextptr); 337 if ((args._VECTptr->front().type==_VECT) && (args._VECTptr->back().type==_VECT)) 338 return axq(*args._VECTptr->front()._VECTptr,*args._VECTptr->back()._VECTptr,contextptr); 339 return symb_a2q(args); 340 } 341 static const char _a2q_s []="a2q"; 342 static define_unary_function_eval (__a2q,&_a2q,_a2q_s); 343 define_unary_function_ptr5( at_a2q ,alias_at_a2q,&__a2q,0,true); 344 qxac(const gen & q,const vecteur & x,GIAC_CONTEXT)345 vecteur qxac(const gen &q,const vecteur & x,GIAC_CONTEXT){ 346 //transforme une forme quadratique en une matrice symetrique A 347 //(les variables sont dans x) 348 int d; 349 //d nbre de variables 350 d=int(x.size()); 351 // int da; 352 //il faut verifier que q est quadratique 353 vecteur A; 354 int b; 355 A=quad(b,q,x,contextptr); 356 if (b==2) { 357 return(A); 358 } 359 else { 360 return vecteur(1,gensizeerr(gettext("q is not quadratic"))); 361 } 362 } 363 364 // rational parametrization of a conic, given cartesian equation and point over conique_ratparam(const gen & eq,const gen & M,GIAC_CONTEXT)365 gen conique_ratparam(const gen & eq,const gen & M,GIAC_CONTEXT){ 366 if (is_undef(M)) 367 return undef; 368 gen Mx,My,x(x__IDNT_e),y(y__IDNT_e),t(t__IDNT_e); 369 if (!contains(eq,x)) 370 ck_parameter_x(contextptr); 371 if (!contains(eq,y)) 372 ck_parameter_y(contextptr); 373 ck_parameter_t(contextptr); 374 reim(M,Mx,My,contextptr); 375 gen eqM=_quo(makesequence(subst(eq,makevecteur(x,y),makevecteur(Mx+x,My+t*x),false,contextptr),x),contextptr); 376 gen a,b; 377 if (!is_linear_wrt(eqM,x,a,b,contextptr)) 378 return undef; 379 return M+(-b/a)*(1+cst_i*t); 380 vecteur res=solve(eqM,x,0,contextptr); // x in terms of t 381 if (res.size()!=1) 382 return undef; 383 return M+res[0]*(1+cst_i*t); 384 } 385 386 // return a,b,c,d,e such that the parametric equation of the conic 387 // is M+(1+i*t)*(d*t+e)/(a*t^2+b*t+c) conique_ratparams(const gen & eq,const gen & M,GIAC_CONTEXT)388 vecteur conique_ratparams(const gen & eq,const gen & M,GIAC_CONTEXT){ 389 if (is_undef(M)) 390 return vecteur(1,undef); 391 gen Mx,My,x(x__IDNT_e),y(y__IDNT_e),t(t__IDNT_e); 392 if (!contains(eq,x)) 393 ck_parameter_x(contextptr); 394 if (!contains(eq,y)) 395 ck_parameter_y(contextptr); 396 ck_parameter_t(contextptr); 397 reim(M,Mx,My,contextptr); 398 gen eqM=_quo(makesequence(subst(eq,makevecteur(x,y),makevecteur(Mx+x,My+t*x),false,contextptr),x),contextptr); 399 gen num,deno,a,b,c,d,e; 400 if (!is_linear_wrt(eqM,x,deno,num,contextptr) || !is_linear_wrt(num,t,d,e,contextptr) || !is_quadratic_wrt(deno,t,a,b,c,contextptr)) 401 return vecteur(1,undef); 402 return makevecteur(at_ellipse,M,a,b,c,-d,-e); 403 } 404 conique_reduite(const gen & equation_conique,const gen & pointsurconique,const vecteur & nom_des_variables,gen & x0,gen & y0,vecteur & V0,vecteur & V1,gen & propre,gen & equation_reduite,vecteur & param_curves,gen & ratparam,bool numeric,GIAC_CONTEXT)405 bool conique_reduite(const gen & equation_conique,const gen & pointsurconique,const vecteur & nom_des_variables,gen & x0, gen & y0, vecteur & V0, vecteur &V1, gen & propre,gen & equation_reduite, vecteur & param_curves,gen & ratparam,bool numeric,GIAC_CONTEXT){ 406 ratparam=conique_ratparam(equation_conique,pointsurconique,contextptr); 407 gen q(remove_equal(equation_conique)); 408 vecteur x(nom_des_variables); 409 if (x.size()!=2) 410 return false; // setsizeerr(contextptr); 411 identificateur iT(" T"); 412 x.push_back(iT); 413 //n est le nombre de variables en geo. projective 414 int n=3; 415 //nom des nouvelles variables 416 vecteur A; 417 gen qp; 418 qp=q; 419 for (int i=0;i<n-1;i++){ 420 qp=subst(qp,x[i],rdiv(x[i],x[2],contextptr),false,contextptr); 421 } 422 qp=recursive_normal(x[2]*x[2]*qp,contextptr); 423 //qp est l'equation en projective qp est quadratique 424 A=qxac(qp,x,contextptr); 425 if (is_undef(A)) 426 return false; 427 #ifdef EMCC 428 // otherwise some implicit plots do not work 429 // but this make distance(point(0,2),y=x^2) approx... 430 if (numeric) 431 A=*evalf(A,1,contextptr)._VECTptr; 432 #endif 433 //q=ax^2+2bxy+cy^2+2dx+2ey+f 434 gen a=A[0][0]; 435 gen b=A[0][1]; 436 gen c=A[1][1]; 437 gen d=A[0][2]; 438 gen e=A[1][2]; 439 gen f=A[2][2]; 440 //propre=(a*c-b*b)*f+d*(b*e-d*c)-e*(a*e-b*d); 441 if ((a*c-b*b)*f+d*(b*e-d*c)-e*(a*e-b*d)!=0) { 442 propre=1; 443 } else { 444 propre=0; 445 } 446 gen X0; 447 gen Y0; 448 gen vp0; 449 gen vp1; 450 V0=vecteur(2); 451 V1=vecteur(2); 452 gen norme; 453 norme=normalize_sqrt(sqrt(a*a+b*b,contextptr),contextptr); 454 if (a*c-b*b==0) { 455 gen coeffy2,coeffx,coeffcst; 456 //on a une parabole ou 2dr// 457 //vecteur propre V0 (vp0 =0) V1 (vp1=a+c) 458 if (a!=0){ 459 //on calcule le nouveau repere origine=(x0,y0) base (V0,V1) 460 //(X0,Y0) sont les coord du centre ds la base (V0,V1) 461 vp0=0; 462 vp1=a+c; 463 V0[0]=rdiv(b,norme,contextptr); 464 V0[1]=rdiv(-a,norme,contextptr); 465 V1[0]=-V0[1]; V1[1]=V0[0]; 466 Y0=-rdiv(d*a+e*b,norme*vp1,contextptr); 467 coeffy2=normal(a+c,contextptr); 468 //cout<<"Y0="<<Y0<<'\n'; 469 coeffx=normal(2*(d*b-e*a)/norme,contextptr); 470 if (coeffx!=0){ 471 // parabole 472 X0=rdiv((d*a+e*b)*(d*a+e*b)-f*(a*a+b*b)*vp1,gen(2)*vp1*(d*b-e*a)*norme,contextptr); 473 equation_reduite=coeffy2*pow(x[1],2)+coeffx*x[0]; 474 } else { 475 //si d*b-e*a==0 alors X0=0 476 coeffcst=normal(f-(d*a+e*b)*(d*a+e*b)/((a+c)*(a*a+b*b)),contextptr); 477 equation_reduite=coeffy2*pow(x[1],2)+coeffcst; 478 X0=0; 479 } 480 } 481 else { 482 // a==0 => b==0 (puisque a*c-b*b==0) et c!=0 483 V0[0]=1; V0[1]=0; 484 V1[0]=0; V1[1]=1; 485 Y0=-rdiv(e,c,contextptr); 486 if (d!=0){ 487 X0=rdiv(e*e-f*c,gen(2)*d*c,contextptr); 488 coeffy2=c; 489 coeffx=normal(2*d,contextptr); 490 equation_reduite=c*pow(x[1],2)+coeffx*x[0]; 491 } else { 492 X0=0; 493 coeffy2=normal(c*c,contextptr); 494 coeffcst=normal(f*c-e*e,contextptr); 495 equation_reduite=coeffy2*pow(x[1],2)+coeffcst; 496 } 497 } 498 x0=normal(V0[0]*X0+V1[0]*Y0,contextptr); 499 y0=normal(V0[1]*X0+V1[1]*Y0,contextptr); 500 gen z0=x0+cst_i*y0; 501 gen zV0=V0[0]+cst_i*V0[1]; 502 if (coeffx==0){ 503 // coeffy2*Y^2+coeffcst=0 504 // 2 lines or empty 505 gen coeff=-coeffcst/coeffy2; 506 gen coeffs=sign(coeff,contextptr); 507 if (is_minus_one(coeffs)){ 508 #ifndef GIAC_HAS_STO_38 509 *logptr(contextptr) << gettext("Empty parabola") << '\n'; 510 #endif 511 return true; 512 } 513 #ifndef GIAC_HAS_STO_38 514 *logptr(contextptr) << gettext("2 parallel lines") << '\n'; 515 #endif 516 gen Y0=normalize_sqrt(sqrt(coeff,contextptr),contextptr); 517 // Y=Y0 : points (0,Y0), (1,Y0) 518 gen zY0=cst_i*zV0*Y0; 519 param_curves.push_back(gen(makevecteur(z0+zY0,z0+zV0+zY0),_LINE__VECT)); 520 param_curves.push_back(gen(makevecteur(z0-zY0,z0+zV0-zY0),_LINE__VECT)); 521 } 522 else { 523 // parabola coeffy2*Y^2+coeffx*X=0 -> X=-coeffy2/coeffx*Y^2 524 // X+i*Y=-coeffy2/coeffx*Y^2+i*Y 525 gen coeff=-coeffy2/coeffx; 526 #ifdef GIAC_HAS_STO_38 527 gen t(vx_var); 528 #else 529 gen t(t__IDNT_e); 530 #endif 531 ck_parameter_t(contextptr); 532 gen Z=coeff*t*t+cst_i*t; 533 Z=z0+zV0*Z; 534 if (is_undef(ratparam)) 535 ratparam=Z; 536 #if defined POCKETCAS 537 param_curves.push_back(makevecteur(Z,t,-10,10,0.1,q,ratparam)); 538 #else 539 param_curves.push_back(makevecteur(Z,t,-4,4,0.1,q,ratparam)); 540 #endif 541 } 542 } 543 else { 544 // a*c-b*b!=0 => on a une conique a centre ou 2 dr concourantes 545 // ellipse/hyperbole 546 if (b==0){ 547 vp0=a; 548 vp1=c; 549 V0[0]=1; V0[1]=0; 550 V1[0]=0; V1[1]=1; 551 } else { 552 //si b!=0 553 gen delta; 554 delta=(a-c)*(a-c)+4*b*b; 555 delta=normalize_sqrt(sqrt(delta,contextptr),contextptr); 556 vp0=ratnormal((a+c+delta)/2,contextptr); 557 vp1=ratnormal((a+c-delta)/2,contextptr); 558 gen normv1(normalize_sqrt(sqrt(b*b+(vp0-a)*(vp0-a),contextptr),contextptr)); 559 V0[0]=normal(b/normv1,contextptr); 560 V0[1]=normal((vp0-a)/normv1,contextptr); 561 V1[0]=-V0[1]; V1[1]=V0[0]; 562 } 563 //coord du centre 564 x0=(-d*c+b*e)/(a*c-b*b); 565 y0=(-a*e+d*b)/(a*c-b*b); 566 gen z0=x0+cst_i*y0; 567 gen zV0=V0[0]+cst_i*V0[1]; 568 gen coeffcst=normal(d*x0+e*y0+f,contextptr); 569 equation_reduite=vp0*pow(x[0],2)+vp1*pow(x[1],2)+ coeffcst; 570 // parametric equations 571 gen svp0(exact(sign(vp0,contextptr),contextptr)), 572 svp1(exact(sign(vp1,contextptr),contextptr)), 573 scoeffcst(exact(sign(coeffcst,contextptr),contextptr)); 574 if (svp0.type==_INT_ && svp1.type==_INT_ && scoeffcst.type==_INT_){ 575 #ifdef GIAC_HAS_STO_38 576 gen t(vx_var); 577 #else 578 gen t(t__IDNT_e); 579 #endif 580 ck_parameter_t(contextptr); 581 int sprodvp = svp0.val * svp1.val; 582 int sprodcoeff = svp0.val*scoeffcst.val; 583 if (sprodvp>0){ // ellipse 584 if (is_zero(coeffcst)){ 585 #ifndef GIAC_HAS_STO_38 586 *logptr(contextptr) << gettext("Ellipsis reduced to (") << x0 << "," << y0 << ")" << '\n'; 587 #endif 588 param_curves.push_back(z0); 589 return true; 590 } 591 if (sprodcoeff>0){ 592 #ifndef GIAC_HAS_STO_38 593 *logptr(contextptr) << gettext("Empty ellipsis") << '\n'; 594 #endif 595 return true; 596 } 597 #ifndef GIAC_HAS_STO_38 598 *logptr(contextptr) << gettext("Ellipsis of center (") << x0 << "," << y0 << ")" << '\n'; 599 #endif 600 vp0=normalize_sqrt(sqrt(-coeffcst/vp0,contextptr),contextptr); 601 vp1=normalize_sqrt(sqrt(-coeffcst/vp1,contextptr),contextptr); 602 // (x[0]/vp0)^2 + (x[1]/vp1)^2 = 1 603 // => x[0]+i*x[1]=vp0*cos(t)+i*vp1*sin(t) 604 // => x+i*y = x0+i*y0 + V0*(x[0]+i*y[0]) 605 gen tmp; 606 if (numeric){ 607 if (!is_undef(ratparam)) 608 tmp=subst(evalf(ratparam,1,contextptr),t,symbolic(at_tan,t/2),false,contextptr); 609 else { 610 tmp=evalf(vp0,1,contextptr)*symb_cos(t)+cst_i*evalf(vp1,1,contextptr)*symb_sin(t); 611 tmp=evalf(z0,1,contextptr)+evalf(zV0,1,contextptr)*tmp; 612 } 613 } 614 else { 615 tmp=vp0*symb_cos(t)+cst_i*vp1*symb_sin(t); 616 tmp=z0+zV0*tmp; 617 } 618 if (is_undef(ratparam)) 619 ratparam=z0+zV0*(vp0*(1-pow(t,2))+cst_i*vp1*plus_two*t)/(1+pow(t,2)); 620 621 bool rad = angle_radian(contextptr), deg = angle_degree(contextptr); 622 param_curves.push_back(makevecteur(tmp,t,0, rad?cst_two_pi:(deg ? 360 : 400), rad?cst_two_pi/60:(deg?6:rdiv(20,3)),q,ratparam)); //grad 623 624 } else { 625 if (is_zero(coeffcst)){ 626 // 2 secant lines at (x0,y0) 627 #ifndef GIAC_HAS_STO_38 628 *logptr(contextptr) << gettext("2 secant lines at (") << x0 << "," << y0 << ")" << '\n'; 629 #endif 630 // vp0*X^2+vp1*Y^2=0 => Y=+/-sqrt(-vp0/vp1)*X 631 gen directeur=normalize_sqrt(sqrt(-vp0/vp1,contextptr),contextptr); 632 if (is_undef(ratparam)) 633 ratparam=makevecteur(z0+zV0*(1+cst_i*directeur)*t,z0+zV0*(1-cst_i*directeur)*t); 634 param_curves.push_back(gen(makevecteur(z0,z0+zV0*(1+cst_i*directeur)),_LINE__VECT)); 635 param_curves.push_back(gen(makevecteur(z0,z0+zV0*(1-cst_i*directeur)),_LINE__VECT)); 636 return true; 637 } 638 // hyperbole 639 #ifndef GIAC_HAS_STO_38 640 *logptr(contextptr) << gettext("Hyperbola of center (") << x0 << "," << y0 << ")" << '\n'; 641 #endif 642 if (sprodcoeff<0) 643 vp0=-vp0; 644 else 645 vp1=-vp1; 646 vp0=normalize_sqrt(sqrt(coeffcst/vp0,contextptr),contextptr); 647 vp1=normalize_sqrt(sqrt(coeffcst/vp1,contextptr),contextptr); 648 gen tmp; 649 if (numeric){ 650 if (!is_undef(ratparam)) 651 tmp=subst(evalf(ratparam,1,contextptr),t,symbolic(at_tan,t/2),false,contextptr); 652 else { 653 tmp=evalf(vp0,1,contextptr)*symbolic(sprodcoeff<0?at_cosh:at_sinh,t)+cst_i*evalf(vp1,1,contextptr)*symbolic(sprodcoeff<0?at_sinh:at_cosh,t); 654 tmp=evalf(z0,1,contextptr)+evalf(zV0,1,contextptr)*tmp; 655 } 656 } 657 else { 658 tmp=vp0*symbolic(sprodcoeff<0?at_cosh:at_sinh,t)+cst_i*vp1*symbolic(sprodcoeff<0?at_sinh:at_cosh,t); 659 tmp=z0+zV0*tmp; 660 } 661 bool noratparam=is_undef(ratparam); 662 if (noratparam){ 663 ratparam=vp0*gen((sprodcoeff<0)?(t+plus_one/t)/2:(t-plus_one/t)/2)+cst_i*vp1*((sprodcoeff<0)?(t-plus_one/t)/2:(t+plus_one/t)/2); 664 ratparam=z0+zV0*ratparam; 665 } 666 #if defined POCKETCAS 667 param_curves.push_back(makevecteur(tmp,t,-10,10,0.1,q,ratparam)); 668 #else 669 param_curves.push_back(makevecteur(tmp,t,-3.14,3.14,0.0314,q,ratparam)); 670 #endif 671 if (noratparam){ 672 if (numeric){ 673 tmp=(sprodcoeff<0?-1:1)*evalf(vp0,1,contextptr)*symbolic(sprodcoeff<0?at_cosh:at_sinh,t)+(sprodcoeff<0?1:-1)*cst_i*evalf(vp1,1,contextptr)*symbolic(sprodcoeff<0?at_sinh:at_cosh,t); 674 tmp=evalf(z0,1,contextptr)+evalf(zV0,1,contextptr)*tmp; 675 } 676 else { 677 tmp=(sprodcoeff<0?-1:1)*vp0*symbolic(sprodcoeff<0?at_cosh:at_sinh,t)+(sprodcoeff<0?1:-1)*cst_i*vp1*symbolic(sprodcoeff<0?at_sinh:at_cosh,t); 678 tmp=z0+zV0*tmp; 679 } 680 #if defined POCKETCAS 681 param_curves.push_back(makevecteur(tmp,t,-10,10,0.1,q,ratparam)); 682 #else 683 param_curves.push_back(makevecteur(tmp,t,-3,3,0.1,q,ratparam)); 684 #endif 685 } 686 } 687 } 688 } 689 return true; 690 } 691 692 #ifdef RTOS_THREADX quadrique_reduite(const gen & q,const gen & M,const vecteur & vxyz,gen & x,gen & y,gen & z,vecteur & u,vecteur & v,vecteur & w,vecteur & propre,gen & equation_reduite,vecteur & param_surface,vecteur & centre,bool numeric,GIAC_CONTEXT)693 bool quadrique_reduite(const gen & q,const gen & M,const vecteur & vxyz,gen & x,gen & y,gen & z,vecteur & u,vecteur & v,vecteur & w,vecteur & propre,gen & equation_reduite,vecteur & param_surface,vecteur & centre,bool numeric,GIAC_CONTEXT){ 694 return false; 695 } 696 #else quadrique_reduite(const gen & q,const gen & M,const vecteur & vxyz,gen & x,gen & y,gen & z,vecteur & u,vecteur & v,vecteur & w,vecteur & propre,gen & equation_reduite,vecteur & param_surface,vecteur & centre,bool numeric,GIAC_CONTEXT)697 bool quadrique_reduite(const gen & q,const gen & M,const vecteur & vxyz,gen & x,gen & y,gen & z,vecteur & u,vecteur & v,vecteur & w,vecteur & propre,gen & equation_reduite,vecteur & param_surface,vecteur & centre,bool numeric,GIAC_CONTEXT){ 698 if (vxyz.size()!=3) 699 return false; // setdimerr(contextptr); 700 x=vxyz[0]; y=vxyz[1]; z=vxyz[2]; 701 identificateur idt("t"); 702 gen t(idt),upar("u",contextptr),vpar("v",contextptr); 703 ck_parameter_u(contextptr); 704 ck_parameter_v(contextptr); 705 gen Q=normal(t*t*(subst(q,vxyz,makevecteur(x/t,y/t,z/t),false,contextptr)),contextptr); // homogeneize 706 matrice A=qxa(Q,makevecteur(x,y,z,t),contextptr); 707 if (is_undef(A)) 708 return false; 709 if (numeric) 710 A=*evalf_double(A,1,contextptr)._VECTptr; 711 // unsigned r=_rank(A).val; 712 matrice B=matrice_extract(A,0,0,3,3); 713 if (is_undef(B)) return false; 714 matrice C=makevecteur(A[0][3],A[1][3],A[2][3]); 715 matrice P; 716 egv(B,P,propre,contextptr,false,false,false); 717 gen s1=propre[0],s2=propre[1],s3=propre[2]; 718 if (is_zero(s1)) s1=0; 719 if (is_zero(s2)) s2=0; 720 if (is_zero(s3)) s3=0; 721 P=mtran(P); 722 P[0]=normal(_normalize(P[0],contextptr),contextptr); 723 P[1]=normal(_normalize(P[1],contextptr),contextptr); 724 P[2]=normal(_normalize(P[2],contextptr),contextptr); 725 u=*P[0]._VECTptr; 726 v=*P[1]._VECTptr; 727 w=*P[2]._VECTptr; 728 if ( s1==0 && s2!=0 ){ 729 vecteur b(u); 730 u=v; v=w; w=b; 731 gen a(s1); 732 s1=s2; s2=s3; s3=a; 733 } 734 if ( s2==0 && s3!=0 ){ 735 vecteur b=w; 736 w=v; v=u; u=b; 737 gen a=s3; 738 s3=s2; s2=s1; s1=a; 739 } 740 gen s1g=evalf_double(sign(s1,contextptr),1,contextptr); 741 gen s2g=evalf_double(sign(s2,contextptr),1,contextptr); 742 gen s3g=evalf_double(sign(s3,contextptr),1,contextptr); 743 if (s1g.type!=_DOUBLE_ || s2g.type!=_DOUBLE_ || s3g.type!=_DOUBLE_){ 744 *logptr(contextptr) << (gettext("Can't check sign ")+s1g.print(contextptr)+gettext(" or ")+s2g.print(contextptr)+gettext(" or ")+s3g.print(contextptr)) << '\n'; 745 return false; 746 } 747 int s1s=int(s1g._DOUBLE_val), s2s=int(s2g._DOUBLE_val), s3s=int(s3g._DOUBLE_val); 748 if (s3!=0){ // hence s1!=0 && s2!=0 749 if (s1s*s2s<0){ 750 if (s1s*s3s>0){ // exchange s2 and s3 751 swap(v,w); 752 swap(s2,s3); 753 swap(s2s,s3s); 754 } 755 else { // s1s and s2s not same sign, s1s and s3s not same sign 756 // therefore s2s and s3s have the same sign 757 vecteur b(u); 758 u=v; v=w; w=b; 759 gen a(s1); 760 s1=s2; s2=s3; s3=a; 761 int as(s1s); 762 s1s=s2s; s2s=s3s; s3s=as; 763 } 764 } 765 // now s1 and s2 have the same sign 766 } 767 P=mtran(makevecteur(u,v,w)); 768 vecteur CP=multvecteurmat(C,P); 769 /* gen CPxyz=dotvecteur(CP,vxyz); 770 gen c=normal(derive(CPxyz,vxyz),contextptr); 771 if (c.type!=_VECT || c._VECTptr->size()!=3) 772 return false; 773 */ 774 gen c=normal(CP,contextptr),d; 775 gen c1(c._VECTptr->front()),c2((*c._VECTptr)[1]),c3(c._VECTptr->back()); 776 gen ustep=_USTEP; 777 ustep.subtype=_INT_PLOT; 778 gen vstep=_VSTEP; 779 vstep.subtype=_INT_PLOT; 780 if (!is_zero(s1)){ 781 if (!is_zero(s2)){ 782 if (!is_zero(s3)){ 783 gen tmp=normal(-c1/s1*u-c2/s2*v-c3/s3*w,contextptr); 784 if (tmp.type!=_VECT || tmp._VECTptr->size()!=3) 785 return false; 786 centre=*tmp._VECTptr; 787 d=normal(subst(q,vxyz,centre,false,contextptr),contextptr); 788 equation_reduite=s1*pow(x,2)+s2*pow(y,2)+s3*pow(z,2)+d; 789 gen dg=evalf_double(sign(d,contextptr),1,contextptr); 790 if (dg.type!=_DOUBLE_) 791 return false; // cksignerr(d); 792 int ds=int(dg._DOUBLE_val); 793 if (ds==0){ 794 // if s3s*s1s>0 solution=1 point, else cone 795 if (s3s*s1s>0) 796 param_surface.push_back(centre); 797 else { // s1*x^2+s2*y^2=-s3*z^2 -> x^2/a^2+y^2/b^2=z^2 798 gen a(sqrt(-s3/s1,contextptr)),b(sqrt(-s3/s2,contextptr)); 799 gen eq=makevecteur(a*upar*symb_cos(vpar),b*upar*symb_sin(vpar),upar); 800 *logptr(contextptr) << gettext("Cone of center ") << centre << '\n'; 801 eq=centre+multmatvecteur(P,*eq._VECTptr); 802 gen ueq=symbolic(at_equal,makesequence(upar,symb_interval(-5,5))); 803 gen veq=symbolic(at_equal,makesequence(vpar,symb_interval(0,cst_two_pi))); 804 ustep=symb_equal(ustep,1./2); 805 vstep=symb_equal(vstep,cst_two_pi/20); 806 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 807 } 808 } // end if ds==0 809 else { 810 if (s1s*s3s>0){ 811 if (s1s*ds>0) 812 *logptr(contextptr) << gettext("Empty ellipsoid") << '\n'; 813 else { 814 gen a=sqrt(-d/s1,contextptr),b=sqrt(-d/s2,contextptr),c=sqrt(-d/s3,contextptr); 815 // x^2/a^2+y^2/b^2+z^2/c^2=1 816 gen eq=makevecteur(a*symb_sin(upar)*symb_cos(vpar),b*symb_sin(upar)*symb_sin(vpar),c*symb_cos(upar)); 817 *logptr(contextptr) << gettext("Ellipsoid of center ") << centre << '\n'; 818 eq=centre+multmatvecteur(P,*eq._VECTptr); 819 gen ueq=symbolic(at_equal,makesequence(upar,symb_interval(0,cst_pi))); 820 gen veq=symbolic(at_equal,makesequence(vpar,symb_interval(0,cst_two_pi))); 821 ustep=symb_equal(ustep,cst_pi/20); 822 vstep=symb_equal(vstep,cst_two_pi/20); 823 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 824 } 825 } // end s1 and s3 of same sign 826 else { // s1 and s2 have same sign, opposite to s3 827 if (s1s*ds>0){ 828 gen a=sqrt(d/s1,contextptr),b=sqrt(d/s2,contextptr),c=sqrt(-d/s3,contextptr); 829 // x^2/a^2+y^2/b^2+1=z^2/c^2, hyperboloide, 2 nappes 830 gen eq=makevecteur(a*symb_sinh(upar)*symb_cos(vpar),b*symb_sinh(upar)*symb_sin(vpar),c*symb_cosh(upar)); 831 eq=centre+multmatvecteur(P,*eq._VECTptr); 832 *logptr(contextptr) << gettext("2-fold hyperboloid of center ") << centre << '\n'; 833 gen ueq=symbolic(at_equal,makesequence(upar,symb_interval(0,3))); 834 gen veq=symbolic(at_equal,makesequence(vpar,symb_interval(0,cst_two_pi))); 835 ustep=symb_equal(ustep,3./20); 836 vstep=symb_equal(vstep,cst_two_pi/20); 837 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 838 eq=makevecteur(a*symb_sinh(upar)*symb_cos(vpar),b*symb_sinh(upar)*symb_sin(vpar),-c*symb_cosh(upar)); 839 eq=centre+multmatvecteur(P,*eq._VECTptr); 840 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 841 } 842 else { // s1, s2 opposite sign to s3,d 843 gen a=sqrt(-d/s1,contextptr),b=sqrt(-d/s2,contextptr),c=sqrt(d/s3,contextptr); 844 // x^2/a^2+y^2/b^2=z^2/c^2+1, hyperboloide, 2 nappes 845 gen eq=makevecteur(a*symb_cosh(upar)*symb_cos(vpar),b*symb_cosh(upar)*symb_sin(vpar),c*symb_sinh(upar)); 846 *logptr(contextptr) << gettext("2-fold hyperboloid of center ") << centre << '\n'; 847 eq=centre+multmatvecteur(P,*eq._VECTptr); 848 gen ueq=symbolic(at_equal,makesequence(upar,symb_interval(-3,3))); 849 gen veq=symbolic(at_equal,makesequence(vpar,symb_interval(0,cst_two_pi))); 850 ustep=symb_equal(ustep,3./20); 851 vstep=symb_equal(vstep,cst_two_pi/20); 852 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 853 } 854 } 855 } 856 return true; 857 } // end if (s3!=0) 858 // s3==0, s1!=0, s2!=0 859 if (is_zero(c3)){ 860 gen tmp=normal(multmatvecteur(P,makevecteur(-c1/s1,-c2/s2,0)),contextptr); 861 if (tmp.type!=_VECT || tmp._VECTptr->size()!=3) return false; 862 centre=*tmp._VECTptr; 863 gen d(normal(subst(q,vxyz,centre,false,contextptr),contextptr)); 864 gen dg=evalf_double(sign(d,contextptr),1,contextptr); 865 if (dg.type!=_DOUBLE_) 866 return false; // cksignerr(d); 867 int ds=int(dg._DOUBLE_val); 868 equation_reduite=s1*pow(x,2)+s2*pow(y,2)+d; 869 if (s1s*s2s>0){ 870 *logptr(contextptr) << gettext("Elliptic cylinder around ") << centre << '\n'; 871 872 if (is_zero(d)) // line (cylinder of radius 0) 873 param_surface.push_back(makevecteur(centre,centre+w)); 874 else { 875 // elliptic cylinder (maybe empty) 876 if (s1s*ds<0){ // s1*x^2+s2*y^2=-d 877 gen a(sqrt(-d/s1,contextptr)),b(sqrt(-d/s2,contextptr)); 878 gen eq=makevecteur(a*symb_cos(vpar),b*symb_sin(vpar),upar); 879 eq=centre+multmatvecteur(P,*eq._VECTptr); 880 gen ueq=symbolic(at_equal,makesequence(upar,symb_interval(-5,5))); 881 gen veq=symbolic(at_equal,makesequence(vpar,symb_interval(0,cst_two_pi))); 882 ustep=symb_equal(ustep,1./2); 883 vstep=symb_equal(vstep,cst_two_pi/20); 884 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 885 } 886 } 887 return true; 888 } // end s1 and s2 of same sign 889 else { // s1 and s2 have opposite signs, s1*x^2+s2*y^2+d=0 890 if (is_zero(d)){ // 2 plans 891 *logptr(contextptr) << gettext("2 plans intersecting at ") << centre << '\n'; 892 gen n=u+sqrt(-s2/s1,contextptr)*v; 893 param_surface.push_back(symbolic(at_hyperplan,gen(makevecteur(n,centre),_SEQ__VECT))); 894 n=u-sqrt(-s2/s1,contextptr)*v; 895 param_surface.push_back(symbolic(at_hyperplan,gen(makevecteur(n,centre),_SEQ__VECT))); 896 return true; 897 } 898 else { // hyperbolic cylinder 899 *logptr(contextptr) << gettext("Hyperbolic cylinder around ") << centre << '\n'; 900 gen ueq=symbolic(at_equal,makesequence(upar,symb_interval(-5,5))); 901 gen veq=symbolic(at_equal,makesequence(vpar,symb_interval(-3,3))); 902 ustep=symb_equal(ustep,1./2); 903 vstep=symb_equal(vstep,0.3); 904 if (s1s*ds<0){ // x^2/(-d/s1) - y^2/(d/s2)=1 905 gen a(sqrt(-d/s1,contextptr)),b(sqrt(d/s2,contextptr)); 906 gen eq=makevecteur(a*symb_cosh(vpar),b*symb_sinh(vpar),upar); 907 eq=centre+multmatvecteur(P,*eq._VECTptr); 908 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 909 eq=makevecteur(-a*symb_cosh(vpar),b*symb_sinh(vpar),upar); 910 eq=centre+multmatvecteur(P,*eq._VECTptr); 911 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 912 return true; 913 } 914 else { // x^2/(d/s1) - y^2/(-d/s2)=-1 915 gen a(sqrt(d/s1,contextptr)),b(sqrt(-d/s2,contextptr)); 916 gen eq=makevecteur(a*symb_sinh(vpar),b*symb_cosh(vpar),upar); 917 eq=centre+multmatvecteur(P,*eq._VECTptr); 918 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 919 eq=makevecteur(a*symb_sinh(vpar),-b*symb_cosh(vpar),upar); 920 eq=centre+multmatvecteur(P,*eq._VECTptr); 921 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 922 return true; 923 } 924 } 925 } 926 } // end c3==0 927 else { 928 // s3==0, s1!=0, s2!=0, c3!=0 929 gen tmp=normal(-c1/s1*u-c2/s2*v,contextptr); 930 if (tmp.type!=_VECT) return false; 931 gen dred=subst(q,vxyz,*tmp._VECTptr,false,contextptr); 932 tmp=normal(tmp-dred/(2*c3)*w,contextptr); 933 if (tmp.type!=_VECT || tmp._VECTptr->size()!=3) return false; 934 centre=*tmp._VECTptr; 935 equation_reduite=s1*pow(x,2)+s2*pow(y,2)+2*c3*z; 936 // parametrization of s1*x^2+s2*y^2+2*c3*z=0 937 if (s1s*s2s>0){ 938 *logptr(contextptr) << gettext("Elliptic paraboloid of center ") << centre << '\n'; 939 // if (s1s*s2s>0) x^2+y^2/(s1/s2)=-2*c3*z/s1 940 // x=u*cos(t), y=u*sqrt(s1/s2)*sin(t), z=-u^2*s1/2/c3 941 gen ueq=symbolic(at_equal,makesequence(upar,symb_interval(0,5))); 942 ustep=symb_equal(ustep,1./2); 943 gen veq=symbolic(at_equal,makesequence(vpar,symb_interval(0,cst_two_pi))); 944 vstep=symb_equal(vstep,cst_two_pi/20); 945 gen a(sqrt(s1/s2,contextptr)),b(-s1/2/c3); 946 gen eq=makevecteur(upar*symb_cos(vpar),a*upar*symb_sin(vpar),b*pow(upar,2)); 947 eq=centre+multmatvecteur(P,*eq._VECTptr); 948 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 949 return true; 950 } 951 else { 952 // if (s1s*s2s<0) x^2-y^2/(-s1/s2)=-2*c3*z/s1 953 *logptr(contextptr) << gettext("Hyperbolic paraboloid of center ") << centre << '\n'; 954 gen ueq=symbolic(at_equal,makesequence(upar,symb_interval(-3,3))); 955 ustep=symb_equal(ustep,0.3); 956 gen veq=symbolic(at_equal,makesequence(vpar,symb_interval(-3,3))); 957 vstep=symb_equal(vstep,0.3); 958 gen a(-s1/s2),b(s1/2/c3); 959 gen eq=makevecteur(upar,sqrt(a,contextptr)*vpar,b*(pow(vpar,2)-pow(upar,2))); 960 eq=centre+multmatvecteur(P,*eq._VECTptr); 961 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 962 return true; 963 } 964 } 965 } // end if (!is_zero(s2)) 966 // here s3==0, s2==0, s1!=0 967 gen tmp=normal(-c1/s1*u,contextptr); 968 // gen tmp=normal(multmatvecteur(P,makevecteur(-c1/s1,0,0)),contextptr); 969 if (tmp.type!=_VECT || tmp._VECTptr->size()!=3) return false; 970 // gen tmp=normal(multmatvecteur(P,makevecteur(-c1/s1,0,0)),contextptr); 971 if (!is_zero(c2) || !is_zero(c3)){ 972 gen c4=normalize_sqrt(sqrt(c2*c2+c3*c3,contextptr),contextptr); 973 gen v1=normal(multmatvecteur(P,makevecteur(0,c3/c4,-c2/c4)),contextptr); 974 gen w1=normal(multmatvecteur(P,makevecteur(0,c2/c4,c3/c4)),contextptr); 975 gen dred=subst(q,vxyz,*tmp._VECTptr,false,contextptr); 976 P=mtran(makevecteur(u,v1,w1)); 977 v=*v1._VECTptr; w=*w1._VECTptr; 978 tmp=tmp+normal(-dred/(2*c4)*w1,contextptr); 979 // tmp=tmp+normal(multmatvecteur(P,makevecteur(0,0,-d/(2*c4))),contextptr); 980 if (tmp.type!=_VECT || tmp._VECTptr->size()!=3) return false; 981 centre=*tmp._VECTptr; 982 // ??? dred=normal(subst(q,vxyz,centre),contextptr); 983 equation_reduite=s1*pow(x,2)+2*c4*z; // ???+dred; 984 gen ueq=symbolic(at_equal,makesequence(upar,symb_interval(-5,5))); 985 gen veq=symbolic(at_equal,makesequence(vpar,symb_interval(-5,5))); 986 ustep=symb_equal(ustep,1./2); 987 vstep=symb_equal(vstep,1./2); 988 *logptr(contextptr) << gettext("Paraboloid cylinder") << '\n'; 989 gen eq=makevecteur(upar,vpar,-s1*pow(upar,2)/2/c4); 990 eq=centre+multmatvecteur(P,*eq._VECTptr); 991 param_surface.push_back(makevecteur(eq,ueq,veq,ustep,vstep)); 992 return true; 993 } 994 else { // c2==0 and c3==0 995 *logptr(contextptr) << gettext("2 parallel plans") << '\n'; 996 centre=*tmp._VECTptr; 997 gen dred=normal(subst(q,vxyz,centre,false,contextptr),contextptr); 998 equation_reduite=s1*pow(x,2)+dred; 999 if (is_zero(dred)){ // a single plan multiplicity 2 1000 param_surface.push_back(symbolic(at_hyperplan,gen(makevecteur(u,centre),_SEQ__VECT))); 1001 } 1002 gen dg=evalf_double(sign(dred,contextptr),1,contextptr); 1003 if (dg.type!=_DOUBLE_) 1004 return false; // cksignerr(d); 1005 int ds=int(dg._DOUBLE_val); 1006 if (s1s*ds<0){ // 2 plans x = +/- sqrt(-dred/s1) 1007 gen a(sqrt(-dred/s1,contextptr)); 1008 param_surface.push_back(symbolic(at_hyperplan,gen(makevecteur(u,centre+a*u),_SEQ__VECT))); 1009 param_surface.push_back(symbolic(at_hyperplan,gen(makevecteur(u,centre-a*u),_SEQ__VECT))); 1010 } 1011 return true; 1012 } 1013 } // end if !is_zero(s1) 1014 return false; 1015 } 1016 #endif // RTOS_THREADX 1017 conique_quadrique_reduite(const gen & args,GIAC_CONTEXT,bool conique)1018 gen conique_quadrique_reduite(const gen & args,GIAC_CONTEXT,bool conique){ 1019 vecteur v(gen2vecteur(args)); 1020 int s=int(v.size()); 1021 if (!s || s>4) 1022 return gendimerr(contextptr); 1023 if (s==4) 1024 v=makevecteur(v[0],makevecteur(v[1],v[2],v[3])); 1025 if (s==3) 1026 v=makevecteur(v[0],makevecteur(v[1],v[2])); 1027 if (s==1){ 1028 v.push_back(conique?makevecteur(x__IDNT_e,y__IDNT_e):makevecteur(x__IDNT_e,y__IDNT_e,z__IDNT_e)); 1029 } 1030 if (v[0].type==_SYMB && v[1].type==_VECT){ 1031 gen x0,y0,z0,eq_reduite,propre,ratparam; 1032 vecteur V0,V1,V2,param_curves,centre,proprev; 1033 if (v[1]._VECTptr->size()==3){ 1034 quadrique_reduite(v[0],undef,*v[1]._VECTptr,x0,y0,z0,V0,V1,V2,proprev,eq_reduite,param_curves,centre,false,contextptr); 1035 return makevecteur(centre,mtran(makevecteur(V0,V1,V2)),proprev,eq_reduite,param_curves); 1036 } 1037 else { 1038 if (!conique_reduite(v[0],undef,*v[1]._VECTptr,x0,y0,V0,V1,propre,eq_reduite,param_curves,ratparam,false,contextptr)) 1039 return gensizeerr(contextptr); 1040 return makevecteur(makevecteur(x0,y0),mtran(makevecteur(V0,V1)),propre,eq_reduite,param_curves); 1041 } 1042 } 1043 return gentypeerr(contextptr); 1044 } _conique_reduite(const gen & args,GIAC_CONTEXT)1045 gen _conique_reduite(const gen & args,GIAC_CONTEXT){ 1046 if ( args.type==_STRNG && args.subtype==-1) return args; 1047 return conique_quadrique_reduite(args,contextptr,true); 1048 } 1049 static const char _conique_reduite_s []="reduced_conic"; 1050 static define_unary_function_eval (__conique_reduite,&_conique_reduite,_conique_reduite_s); 1051 define_unary_function_ptr5( at_conique_reduite ,alias_at_conique_reduite,&__conique_reduite,0,true); 1052 _quadrique_reduite(const gen & args,GIAC_CONTEXT)1053 gen _quadrique_reduite(const gen & args,GIAC_CONTEXT){ 1054 if ( args.type==_STRNG && args.subtype==-1) return args; 1055 return conique_quadrique_reduite(args,contextptr,false); 1056 } 1057 static const char _quadrique_reduite_s []="reduced_quadric"; 1058 static define_unary_function_eval (__quadrique_reduite,&_quadrique_reduite,_quadrique_reduite_s); 1059 define_unary_function_ptr5( at_quadrique_reduite ,alias_at_quadrique_reduite,&__quadrique_reduite,0,true); 1060 1061 1062 #ifndef NO_NAMESPACE_GIAC 1063 } // namespace giac 1064 #endif // ndef NO_NAMESPACE_GIAC 1065 1066