1 /* -*- mode:C++ ; compile-command: "g++-3.4 -I. -I.. -I../include -g -c plot3d.cc -DIN_GIAC -DHAVE_CONFIG_H " -*- */ 2 // NB: Using gnuplot optimally requires patching and recompiling gnuplot 3 // If you use the -DGNUPLOT_IO compile flag, you 4 // MUST compile gnuplot with interactive mode enabled, file src/plot.c 5 // line 448 6 /* 7 diff plot.c plot.c~ 8 448c448 9 < interactive = TRUE; // isatty(fileno(stdin)); 10 --- 11 > interactive = isatty(fileno(stdin)); 12 */ 13 14 /* 15 * Copyright (C) 2000/2014 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 16 * implicitplot3d code adapted from 17 * http://astronomy.swin.edu.au/~pbourke/modelling/polygonise 18 * by Paul Bourke and Cory Gene Bloyd 19 * 20 * This program is free software; you can redistribute it and/or modify 21 * it under the terms of the GNU General Public License as published by 22 * the Free Software Foundation; either version 3 of the License, or 23 * (at your option) any later version. 24 * 25 * This program is distributed in the hope that it will be useful, 26 * but WITHOUT ANY WARRANTY; without even the implied warranty of 27 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 28 * GNU General Public License for more details. 29 * 30 * You should have received a copy of the GNU General Public License 31 * along with this program. If not, see <http://www.gnu.org/licenses/>. 32 */ 33 34 #include "giacPCH.h" 35 36 using namespace std; 37 #if !defined NSPIRE && !defined FXCG && !defined KHICAS 38 #if defined VISUALC13 && !defined BESTA_OS 39 #undef clock 40 #undef clock_t 41 #endif 42 #include <iomanip> 43 #endif 44 #include "vector.h" 45 #include <algorithm> 46 #include <cmath> 47 48 // C headers 49 #include <stdio.h> 50 51 // Giac headers 52 #include "gen.h" 53 #include "usual.h" 54 #include "plot.h" 55 #include "plot3d.h" 56 #include "prog.h" 57 #include "rpn.h" 58 #include "identificateur.h" 59 #include "subst.h" 60 #include "symbolic.h" 61 #include "derive.h" 62 #include "solve.h" 63 #include "intg.h" 64 #include "path.h" 65 #include "sym2poly.h" 66 #include "input_parser.h" 67 #include "input_lexer.h" 68 #include "ti89.h" 69 #include "isom.h" 70 #include "ifactor.h" 71 #include "gauss.h" 72 #include "giacintl.h" 73 74 #ifndef NO_NAMESPACE_GIAC 75 namespace giac { 76 #endif // ndef NO_NAMESPACE_GIAC 77 do_point3d(const gen & g)78 gen do_point3d(const gen & g){ 79 gen tmp(g); 80 if (tmp.type==_VECT) 81 tmp.subtype=_POINT__VECT; 82 return tmp; 83 } 84 rand_3d()85 vecteur rand_3d(){ 86 int i=std_rand(),j=std_rand(),k=std_rand(); 87 i=i/(RAND_MAX/10)-5; 88 j=j/(RAND_MAX/10)-5; 89 k=k/(RAND_MAX/10)-5; 90 return makevecteur(i,j,k); 91 } 92 hyperplan_normal(const gen & g)93 vecteur hyperplan_normal(const gen & g){ 94 vecteur n,P; 95 if (!hyperplan_normal_point(g,n,P)) 96 return vecteur(3,gensizeerr(gettext("hyperplan_normal"))); 97 return n; 98 } hyperplan_normal_point(const gen & g,vecteur & n,vecteur & P)99 bool hyperplan_normal_point(const gen & g,vecteur & n,vecteur & P){ 100 gen h=remove_at_pnt(g); 101 if (h.is_symb_of_sommet(at_hyperplan)) 102 h=h._SYMBptr->feuille; 103 if (h.type!=_VECT || h._VECTptr->size()!=2 || h._VECTptr->front().type!=_VECT || h._VECTptr->back().type!=_VECT) 104 return false; // setsizeerr(contextptr); 105 n=*h._VECTptr->front()._VECTptr; 106 P=*h._VECTptr->back()._VECTptr; 107 return true; 108 } 109 remove_pnt_vect(const gen & g)110 gen remove_pnt_vect(const gen & g){ 111 gen res=remove_at_pnt(g); 112 if (res.type==_VECT && res.subtype==_VECTOR__VECT && res._VECTptr->size()==2) 113 res=res._VECTptr->back()-res._VECTptr->front(); 114 return res; 115 } 116 _plan(const gen & args,GIAC_CONTEXT)117 gen _plan(const gen & args,GIAC_CONTEXT){ 118 if ( args.type==_STRNG && args.subtype==-1) return args; 119 if (args.type==_INT_ || (args.type==_VECT && args._VECTptr->empty()) ) 120 return mkrand2d3d(3,3,_plan,contextptr); 121 if (args.type==_SYMB) 122 return droite_by_equation(vecteur(1,args),true,contextptr); 123 if (args.type!=_VECT || args._VECTptr->size()<2) 124 return gensizeerr(contextptr); 125 vecteur v = *args._VECTptr; 126 vecteur attributs(1,default_color(contextptr)); 127 int s=read_attributs(v,attributs,contextptr); 128 v=vecteur(v.begin(),v.begin()+s); 129 if (!v.empty() && is_equal(v.front())) 130 return droite_by_equation(*args._VECTptr,true,contextptr); 131 if (s) 132 v[0]=remove_at_pnt(v[0]); 133 if (s>1) 134 v[1]=remove_pnt_vect(v[1]); 135 if (s==2){ 136 if (v[0].type==_VECT && v[0]._VECTptr->size()==2 && v[1].type==_VECT && v[1]._VECTptr->size()==2){ 137 // plane in space defined by 2 lines: must be parallel or secant 138 gen A=v[0]._VECTptr->front(),B=v[0]._VECTptr->back(), 139 C=v[1]._VECTptr->front(),D=v[1]._VECTptr->back(); 140 if (!check3dpoint(A) || !check3dpoint(B) || !check3dpoint(C) || !check3dpoint(D)) 141 return gensizeerr(contextptr); 142 vecteur v1(subvecteur(*B._VECTptr,*A._VECTptr)); 143 vecteur v2(subvecteur(*D._VECTptr,*C._VECTptr)); 144 gen M,N,coeff; 145 vecteur n; 146 if (est_parallele_vecteur(v1,v2,coeff,contextptr)) 147 return _plan(gen(makevecteur(A,B,C),_SEQ__VECT),contextptr); 148 if (perpendiculaire_commune(v[0],v[1],M,N,n,contextptr) && is_zero(M-N)) 149 return pnt_attrib(symbolic(at_hyperplan,gen(makevecteur(n,M),_SEQ__VECT)),attributs,contextptr); 150 return gensizeerr(contextptr); 151 } 152 if (v[1].type==_VECT && v[1]._VECTptr->size()==2){ 153 s++; 154 v.push_back(v[1]._VECTptr->back()); 155 v[1]=v[1]._VECTptr->front(); 156 } 157 else 158 return pnt_attrib(symbolic(at_hyperplan,gen(v,args.subtype)),attributs,contextptr); 159 } 160 if (s==3){ 161 v[2]=remove_pnt_vect(v[2]); 162 if (v[0].type==_VECT && v[0]._VECTptr->size()==3 && v[1].type==_VECT && v[1]._VECTptr->size()==3 && v[2].type==_VECT && v[2]._VECTptr->size()==3){ 163 // given by 3 points, compute normal vector 164 gen v1=v[1]-v[0]; 165 gen v2=v[2]-v[0]; 166 gen n=cross(*v1._VECTptr,*v2._VECTptr,contextptr); 167 return pnt_attrib(symbolic(at_hyperplan,gen(makevecteur(n,v[0]),_SEQ__VECT)),attributs,contextptr); 168 } 169 } 170 return gensizeerr(contextptr); 171 } 172 static const char _plan_s []="plane"; 173 static define_unary_function_eval (__plan,&_plan,_plan_s); 174 define_unary_function_ptr5( at_plan ,alias_at_plan,&__plan,0,true); 175 176 // args=center,radius _sphere(const gen & args,GIAC_CONTEXT)177 gen _sphere(const gen & args,GIAC_CONTEXT){ 178 if ( args.type==_STRNG && args.subtype==-1) return args; 179 // if ( (args.type==_SYMB) && (args._SYMBptr->sommet==at_equal) ) 180 // return sphere_by_equation(vecteur(1,args),contextptr); 181 if (args.is_symb_of_sommet(at_equal)) 182 return _plotimplicit(makesequence(args,x__IDNT_e,y__IDNT_e,z__IDNT_e),contextptr); 183 if (args.type!=_VECT || args._VECTptr->size()<2) 184 return gensizeerr(contextptr); 185 gen errcode=checkanglemode(contextptr); 186 if (is_undef(errcode)) return errcode; 187 // if ((args._VECTptr->size()>=2) && (args._VECTptr->front().type==_SYMB) && (args._VECTptr->front()._SYMBptr->sommet==at_equal)) 188 // return sphere_by_equation(*args._VECTptr,contextptr); 189 vecteur v = *args._VECTptr; 190 vecteur attributs(1,default_color(contextptr)); 191 int s=read_attributs(v,attributs,contextptr); 192 v=vecteur(v.begin(),v.begin()+s); 193 v[0]=remove_at_pnt(v[0]); 194 if (v[0].type!=_VECT) 195 v[0]=gen(makevecteur(v[0],0,0),_POINT__VECT); 196 v[1]=remove_at_pnt(v[1]); 197 if (v[1].type==_VECT){ 198 gen tmp=v[1]; 199 if (v[1].subtype==_POINT__VECT){ 200 tmp=(v[1]-v[0])/2; 201 if (tmp.type!=_VECT) 202 return gensizeerr(contextptr); 203 v[0]=(v[0]+v[1])/2; 204 } 205 v[1]=l2norm(*tmp._VECTptr,contextptr); 206 } 207 else { 208 if (is_strictly_positive(-v[1],contextptr)) 209 return gensizeerr(contextptr); 210 } 211 return pnt_attrib(symbolic(at_hypersphere,gen(v,args.subtype)),attributs,contextptr); 212 } 213 static const char _sphere_s []="sphere"; 214 static define_unary_function_eval (__sphere,&_sphere,_sphere_s); 215 define_unary_function_ptr5( at_sphere ,alias_at_sphere,&__sphere,0,true); 216 option_adjust(int & nstep,int & jstep,int & kstep)217 static void option_adjust(int & nstep,int & jstep,int & kstep){ 218 if (nstep){ 219 jstep=int(std::sqrt(double(nstep))); 220 kstep=int(std::sqrt(double(nstep))); 221 } 222 if (kstep<1) 223 kstep=10; 224 if (jstep<1) 225 jstep=10; 226 } 227 cone(const gen & args,bool cone_complet,GIAC_CONTEXT)228 gen cone(const gen & args,bool cone_complet,GIAC_CONTEXT){ 229 if (args.type!=_VECT) 230 return gensizeerr(contextptr); 231 vecteur attributs(1,default_color(contextptr)); 232 int s=read_attributs(*args._VECTptr,attributs,contextptr); 233 double xmin=gnuplot_xmin,xmax=gnuplot_xmax,ymin=gnuplot_ymin,ymax=gnuplot_ymax,zmin=gnuplot_zmin,zmax=gnuplot_zmax; 234 int jstep=0,kstep=0,nstep=0; 235 vecteur vtmp; 236 read_option(*args._VECTptr,xmin,xmax,ymin,ymax,zmin,zmax,vtmp,nstep,jstep,kstep,contextptr); 237 option_adjust(nstep,jstep,kstep); 238 if (s<3) 239 return gensizeerr(contextptr); 240 gen errcode=checkanglemode(contextptr); 241 if (is_undef(errcode)) return errcode; 242 ck_parameter_x(contextptr); 243 ck_parameter_y(contextptr); 244 ck_parameter_z(contextptr); 245 vecteur v = *args._VECTptr; 246 gen P=remove_at_pnt(v[0]),theta=v[2]; 247 v[1]=remove_pnt_vect(v[1]); 248 if (v[1].type!=_VECT || P.type!=_VECT) 249 return gensizeerr(contextptr); 250 vecteur xyz(makevecteur(x__IDNT_e,y__IDNT_e,z__IDNT_e)); 251 vecteur xyzP=subvecteur(xyz,*P._VECTptr); 252 vecteur n=*v[1]._VECTptr,n1,n2; 253 if (!normal3d(v[1],n1,n2)) 254 return gensizeerr(contextptr); 255 n=divvecteur(n,abs_norm(n,contextptr)); 256 n1=divvecteur(n1,abs_norm(n1,contextptr)); 257 n2=divvecteur(n2,abs_norm(n2,contextptr)); 258 gen eq=normal(pow(dotvecteur(xyzP,n),2)*pow(sin(theta,contextptr),2)-(pow(dotvecteur(xyzP,n1),2)+pow(dotvecteur(xyzP,n2),2))*pow(cos(theta,contextptr),2),contextptr); 259 gen M=P+u__IDNT_e*(cos(theta,contextptr)*n+sin(theta,contextptr)*(cos(v__IDNT_e,contextptr)*n1+sin(v__IDNT_e,contextptr)*n2)); 260 double uscale=gnuplot_tmax-gnuplot_tmin; 261 bool cercles=false; 262 if (s>3){ 263 gen tmp=evalf_double(v[3]/cos(theta,contextptr),eval_level(contextptr),contextptr); 264 if (tmp.type==_DOUBLE_){ 265 uscale=tmp._DOUBLE_val; 266 cercles=true; 267 } 268 } 269 vecteur uv(makevecteur(u__IDNT_e,v__IDNT_e)); 270 gen res= plotparam3d(M,uv,xmin,xmax,ymin,ymax,zmin,zmax,cone_complet?-uscale:0,uscale,0,2*M_PI,false,false,attributs,uscale/jstep,M_PI/kstep,eq,xyz,contextptr); 271 if (!cercles) 272 return res; 273 theta=evalf_double(theta,1,contextptr); 274 if (theta.type!=_DOUBLE_) 275 return res; 276 double thetad=theta._DOUBLE_val; 277 // add disque center P+uscale*cos(theta)*n, perp to n, r=uscale*sin(theta) 278 vecteur vres(1,res); 279 M=P+uscale*cos(theta,contextptr)*n+u__IDNT_e*(cos(v__IDNT_e,contextptr)*n1+sin(v__IDNT_e,contextptr)*n2); 280 res=plotparam3d(M,uv,xmin,xmax,ymin,ymax,zmin,zmax,0,uscale*std::sin(thetad),0,2*M_PI,false,false,attributs,uscale*std::sin(thetad),M_PI/kstep,undef,xyz,contextptr); 281 vres.push_back(res); 282 if (cone_complet){ 283 M=P-uscale*cos(theta,contextptr)*n+u__IDNT_e*(cos(v__IDNT_e,contextptr)*n1+sin(v__IDNT_e,contextptr)*n2); 284 res=plotparam3d(M,uv,xmin,xmax,ymin,ymax,zmin,zmax,0,uscale*std::sin(thetad),0,2*M_PI,false,false,attributs,uscale*std::sin(thetad),M_PI/kstep,undef,xyz,contextptr); 285 vres.push_back(res); 286 } 287 return vres; // gen(vres,_SEQ__VECT); 288 } 289 // args=point, direction, angle _cone(const gen & args,GIAC_CONTEXT)290 gen _cone(const gen & args,GIAC_CONTEXT){ 291 if ( args.type==_STRNG && args.subtype==-1) return args; 292 return cone(args,true,contextptr); 293 } 294 static const char _cone_s []="cone"; 295 static define_unary_function_eval (__cone,&_cone,_cone_s); 296 define_unary_function_ptr5( at_cone ,alias_at_cone,&__cone,0,true); 297 298 // args=point, direction, angle _demi_cone(const gen & args,GIAC_CONTEXT)299 gen _demi_cone(const gen & args,GIAC_CONTEXT){ 300 if ( args.type==_STRNG && args.subtype==-1) return args; 301 return cone(args,false,contextptr); 302 } 303 static const char _demi_cone_s []="half_cone"; 304 static define_unary_function_eval (__demi_cone,&_demi_cone,_demi_cone_s); 305 define_unary_function_ptr5( at_demi_cone ,alias_at_demi_cone,&__demi_cone,0,true); 306 307 // args=point, direction, radius _cylindre(const gen & args,GIAC_CONTEXT)308 gen _cylindre(const gen & args,GIAC_CONTEXT){ 309 if ( args.type==_STRNG && args.subtype==-1) return args; 310 if (args.type!=_VECT) 311 return gensizeerr(contextptr); 312 vecteur attributs(1,default_color(contextptr)); 313 int s=read_attributs(*args._VECTptr,attributs,contextptr); 314 double xmin=gnuplot_xmin,xmax=gnuplot_xmax,ymin=gnuplot_ymin,ymax=gnuplot_ymax,zmin=gnuplot_zmin,zmax=gnuplot_zmax; 315 int jstep=0,kstep=0,nstep=0; 316 vecteur vtmp; 317 read_option(*args._VECTptr,xmin,xmax,ymin,ymax,zmin,zmax,vtmp,nstep,jstep,kstep,contextptr); 318 option_adjust(nstep,jstep,kstep); 319 if (s<3) 320 return gensizeerr(contextptr); 321 double uscale=gnuplot_tmax-gnuplot_tmin; 322 bool cercles=false; 323 vecteur v = *args._VECTptr; 324 if (s>3){ 325 gen tmp=evalf_double(v[3],eval_level(contextptr),contextptr); 326 if (tmp.type==_DOUBLE_){ 327 uscale=tmp._DOUBLE_val; 328 cercles=true; 329 } 330 } 331 gen errcode=checkanglemode(contextptr); 332 if (is_undef(errcode)) return errcode; 333 ck_parameter_x(contextptr); 334 ck_parameter_y(contextptr); 335 ck_parameter_z(contextptr); 336 gen P=remove_at_pnt(v[0]),r=v[2]; 337 v[1]=remove_pnt_vect(v[1]); 338 if (v[1].type!=_VECT || P.type!=_VECT) 339 return gensizeerr(contextptr); 340 vecteur xyz(makevecteur(x__IDNT_e,y__IDNT_e,z__IDNT_e)); 341 vecteur xyzP=subvecteur(xyz,*P._VECTptr); 342 vecteur n=*v[1]._VECTptr,n1,n2; 343 if (!normal3d(v[1],n1,n2)) 344 return gensizeerr(contextptr); 345 n=divvecteur(n,abs_norm(n,contextptr)); 346 n1=divvecteur(n1,abs_norm(n1,contextptr)); 347 n2=divvecteur(n2,abs_norm(n2,contextptr)); 348 gen M=P+u__IDNT_e*n+r*(cos(v__IDNT_e,contextptr)*n1+sin(v__IDNT_e,contextptr)*n2); 349 gen eq=normal(pow(r,2)-(pow(dotvecteur(xyzP,n1),2)+pow(dotvecteur(xyzP,n2),2)),contextptr); 350 vecteur uv(makevecteur(u__IDNT_e,v__IDNT_e)); 351 gen res=plotparam3d(M,uv,xmin,xmax,ymin,ymax,zmin,zmax,0,uscale,0,2*M_PI,false,false,attributs,uscale/jstep,M_PI/kstep,eq,xyz,contextptr); 352 if (!cercles) 353 return res; 354 // add disque center P and P+uscale*n, perp to n, radius r 355 r=evalf_double(r,1,contextptr); 356 if (r.type!=_DOUBLE_) 357 return res; 358 double rd=r._DOUBLE_val; 359 vecteur vres(1,res); 360 M=P+u__IDNT_e*(cos(v__IDNT_e,contextptr)*n1+sin(v__IDNT_e,contextptr)*n2); 361 res=plotparam3d(M,uv,xmin,xmax,ymin,ymax,zmin,zmax,0,rd,0,2*M_PI,false,false,attributs,rd,M_PI/kstep,undef,xyz,contextptr); 362 vres.push_back(res); 363 M=M+n*gen(uscale); 364 res=plotparam3d(M,uv,xmin,xmax,ymin,ymax,zmin,zmax,0,rd,0,2*M_PI,false,false,attributs,rd,M_PI/kstep,undef,xyz,contextptr); 365 vres.push_back(res); 366 return vres; // gen(vres,_SEQ__VECT); 367 } 368 static const char _cylindre_s []="cylinder"; 369 static define_unary_function_eval (__cylindre,&_cylindre,_cylindre_s); 370 define_unary_function_ptr5( at_cylindre ,alias_at_cylindre,&__cylindre,0,true); 371 372 // find the 2 points of d1 and d2 and a common normal vector to d1 d2 perpendiculaire_commune(const gen & d1,const gen & d2,gen & M,gen & N,vecteur & n,GIAC_CONTEXT)373 bool perpendiculaire_commune(const gen & d1,const gen & d2,gen & M, gen & N,vecteur & n,GIAC_CONTEXT){ 374 gen D1=remove_at_pnt(d1); 375 gen D2=remove_at_pnt(d2); 376 if (D1.type!=_VECT || D1._VECTptr->size()!=2 || D2.type!=_VECT || D2._VECTptr->size()!=2) 377 return false; 378 gen & A=D1._VECTptr->front(); 379 gen & B=D1._VECTptr->back(); 380 gen & C=D2._VECTptr->front(); 381 gen & D=D2._VECTptr->back(); 382 if (!check3dpoint(A)){ 383 return false; 384 } 385 if (!check3dpoint(B)) 386 return false; 387 if (!check3dpoint(C)) 388 return false; 389 if (!check3dpoint(D)) 390 return false; 391 vecteur v1(subvecteur(*B._VECTptr,*A._VECTptr)); 392 vecteur v2(subvecteur(*D._VECTptr,*C._VECTptr)); 393 n=*normal(cross(v1,v2,contextptr),contextptr)._VECTptr; 394 if (is_zero(n)) 395 return false; 396 // M=A+u*v1, N=C-v*v2, find u and v such that 397 // NM.v1=0 and NM.v2=0 398 // where NM= -(C-A)+ u*v1+v*v2 399 // Hence solve the system of matrix 400 // v1.v1*u + v2.v1*v = (C-A).v1 401 // v1.v2*u + v2.v2*v = (C-A).v2 402 vecteur AC(subvecteur(*C._VECTptr,*A._VECTptr)); 403 gen v11(dotvecteur(v1,v1)),v22(dotvecteur(v2,v2)),v12(dotvecteur(v1,v2)); 404 gen AC1(dotvecteur(v1,AC)),AC2(dotvecteur(v2,AC)); 405 gen det(v11*v22-v12*v12); 406 gen u= (v22*AC1-v12*AC2)/det,v=(v11*AC2-v12*AC1)/det; 407 M=A+u*v1; 408 N=C-v*v2; 409 M.subtype=_POINT__VECT; 410 N.subtype=_POINT__VECT; 411 return true; 412 } _perpendiculaire_commune(const gen & args,GIAC_CONTEXT)413 gen _perpendiculaire_commune(const gen & args,GIAC_CONTEXT){ 414 if ( args.type==_STRNG && args.subtype==-1) return args; 415 if ( (args.type!=_VECT) || (args._VECTptr->size()<2)) 416 return gensizeerr(contextptr); 417 vecteur attributs(1,default_color(contextptr)); 418 read_attributs(*args._VECTptr,attributs,contextptr); 419 gen M,N; 420 vecteur n; 421 if (!perpendiculaire_commune(args._VECTptr->front(),args._VECTptr->back(),M,N,n,contextptr)) 422 return gensizeerr(gettext("Parallel lines")); 423 return pnt_attrib(gen(makevecteur(M,N),_LINE__VECT),attributs,contextptr); 424 } 425 static const char _perpendiculaire_commune_s []="common_perpendicular"; 426 static define_unary_function_eval (__perpendiculaire_commune,&_perpendiculaire_commune,_perpendiculaire_commune_s); 427 define_unary_function_ptr5( at_perpendiculaire_commune ,alias_at_perpendiculaire_commune,&__perpendiculaire_commune,0,true); 428 429 gen _polyedre(const gen & args,GIAC_CONTEXT); 430 431 // Given a list of 3-d points, make a convex polyedre polyedre(const gen & g,GIAC_CONTEXT)432 static vecteur polyedre(const gen & g,GIAC_CONTEXT){ 433 if (g.type!=_VECT || g._VECTptr->size()<3) 434 return vecteur(1,gensizeerr(contextptr)); 435 vecteur v =*g._VECTptr; 436 // Construct faces: easy algorithm 437 // Make all possibles plans with 3 points, find equation 438 // Add to the face any point that is in the plane (eq=0) 439 // All other points must have the same sign 440 // Otherwise break and try next triple of points 441 iterateur i=v.begin(),ie=v.end(); 442 for (;i!=ie;++i){ 443 *i = remove_at_pnt(*i); 444 if (i->type!=_VECT || i->_VECTptr->size()!=3) 445 return vecteur(1,gendimerr(contextptr)); 446 } 447 vecteur faces; 448 for (i=v.begin();i!=ie;++i){ 449 const_iterateur j=i+1; 450 for (;j!=ie;++j){ 451 const_iterateur k=j+1; 452 vecteur v1(subvecteur(*j->_VECTptr,*i->_VECTptr)); 453 for (;k!=ie;++k){ 454 vecteur v2(subvecteur(*k->_VECTptr,*i->_VECTptr)); 455 vecteur n(*normal(cross(v1,v2,contextptr),contextptr)._VECTptr); 456 if (is_zero(n)) 457 continue; 458 const_iterateur l=v.begin(); 459 gen s; // sign 460 gen eq,eqs; 461 vecteur currentface; 462 for (;l!=ie;++l){ 463 if (l==i || l==j || l==k){ 464 currentface.push_back(*l); 465 continue; 466 } 467 eq=normal(dotvecteur(n,subvecteur(*l->_VECTptr,*i->_VECTptr)),contextptr); 468 if (is_zero(eq)){ 469 currentface.push_back(*l); 470 continue; 471 } 472 eqs=evalf_double(sign(eq,contextptr),1,contextptr); 473 if (eqs.type!=_DOUBLE_) 474 return vecteur(1,gensizeerr(gettext("Unable to check sign ")+eq.print(contextptr))); 475 if (is_zero(s)) 476 s=eqs; 477 if (eqs!=s) 478 break; 479 } // end for l 480 if (l==ie) 481 faces.push_back(currentface); 482 } 483 } 484 } 485 return faces; 486 } polyedre_face(vecteur & v,const vecteur & attributs,GIAC_CONTEXT)487 static gen polyedre_face(vecteur & v,const vecteur & attributs,GIAC_CONTEXT){ 488 iterateur it=v.begin(),itend=v.end(); 489 for (;it!=itend;++it){ 490 if (it->type!=_VECT) 491 return gensizeerr(gettext("Each element must be a face (vector of 3/4 points)")); 492 vecteur w(*it->_VECTptr); 493 int s=int(w.size()); 494 if (s<3) 495 return gensizeerr(gettext("at least 3 points by face")); 496 iterateur jt=w.begin(),jtend=w.end(); 497 for (;jt!=jtend;++jt){ 498 *jt=remove_at_pnt(*jt); 499 if (jt->type==_VECT) 500 jt->subtype=_POINT__VECT; 501 } 502 *it=w; 503 it->subtype=_GROUP__VECT; 504 } 505 return pnt_attrib(gen(v,_POLYEDRE__VECT),attributs,contextptr); 506 } _polyedre(const gen & args,GIAC_CONTEXT)507 gen _polyedre(const gen & args,GIAC_CONTEXT){ 508 if ( args.type==_STRNG && args.subtype==-1) return args; 509 if (args.type!=_VECT) 510 return gensizeerr(contextptr); 511 vecteur v(*args._VECTptr); 512 vecteur attributs(1,default_color(contextptr)); 513 int s=read_attributs(v,attributs,contextptr); 514 v=vecteur(v.begin(),v.begin()+s); 515 if (!s) 516 return gendimerr(contextptr); 517 if (s==2){ 518 // Base, sommet 519 gen base=remove_at_pnt(v.front()); 520 if (base.type!=_VECT) 521 return gensizeerr(contextptr); 522 gen & sommet=v.back(); 523 vecteur w=*base._VECTptr; 524 vecteur nv; 525 if (w.front()!=w.back()) 526 w.push_back(w.front()); 527 int s=int(w.size()); 528 if (s<3) 529 return gendimerr(contextptr); 530 for (int i=0;i<s;++i){ 531 nv.push_back(makevecteur(w[i],w[(i+1)%s],sommet)); 532 } 533 nv.push_back(base); 534 return polyedre_face(nv,attributs,contextptr); 535 } 536 gen g=remove_at_pnt(v.front()); 537 if (g.type==_VECT && g._VECTptr->size()==3 && remove_at_pnt(g._VECTptr->front()).type!=_VECT ) 538 v=polyedre(v,contextptr); 539 return polyedre_face(v,attributs,contextptr); 540 } 541 static const char _polyedre_s []="polyhedron"; 542 static define_unary_function_eval (__polyedre,&_polyedre,_polyedre_s); 543 define_unary_function_ptr5( at_polyedre ,alias_at_polyedre,&__polyedre,0,true); 544 _prisme(const gen & args,GIAC_CONTEXT)545 gen _prisme(const gen & args,GIAC_CONTEXT){ 546 if ( args.type==_STRNG && args.subtype==-1) return args; 547 if (args.type!=_VECT) 548 return gensizeerr(contextptr); 549 vecteur & v=*args._VECTptr; 550 vecteur attributs(1,default_color(contextptr)); 551 int sv=read_attributs(v,attributs,contextptr); 552 if (sv!=2) 553 return gendimerr(contextptr); 554 gen base=remove_at_pnt(v[0]), sommet=remove_at_pnt(v[1]); 555 if (base.type!=_VECT || base._VECTptr->size()<2) 556 return gensizeerr(contextptr); 557 vecteur w = *base._VECTptr; 558 gen x=sommet-w[0]; 559 int s=int(w.size()); 560 vecteur faces; 561 for (int i=0;i<s;++i){ 562 faces.push_back(makevecteur(w[i],w[(i+1)%s],w[(i+1)%s]+x,w[i]+x)); 563 } 564 faces.push_back(base); 565 for (int i=0;i<s;++i) 566 w[i]=w[i]+x; 567 faces.push_back(w); 568 return polyedre_face(faces,attributs,contextptr); 569 } 570 static const char _prisme_s []="prism"; 571 static define_unary_function_eval (__prisme,&_prisme,_prisme_s); 572 define_unary_function_ptr5( at_prisme ,alias_at_prisme,&__prisme,0,true); 573 parallelepipede4(const gen & A0,const gen & B0,const gen & C0,const gen & D0,const vecteur & attributs,GIAC_CONTEXT)574 static gen parallelepipede4(const gen & A0,const gen & B0,const gen & C0,const gen & D0,const vecteur & attributs,GIAC_CONTEXT){ 575 gen A(A0),B(B0),C(C0),D(D0); 576 A.subtype=_POINT__VECT; 577 B.subtype=_POINT__VECT; 578 C.subtype=_POINT__VECT; 579 D.subtype=_POINT__VECT; 580 gen AB=B-A,AC=C-A,AD=D-A; 581 gen E=A+AB+AC,F=A+AC+AD,G=A+AB+AD,H=A+AB+AC+AD; 582 E.subtype=_POINT__VECT; 583 F.subtype=_POINT__VECT; 584 G.subtype=_POINT__VECT; 585 H.subtype=_POINT__VECT; 586 vecteur res; 587 // Face 1 A B // C E=A+AB+AC 588 res.push_back(makevecteur(A,B,E,C)); 589 // Face 6 D G // F H 590 res.push_back(makevecteur(D,G,H,F)); 591 // Face 2 A C // D F=A+AC+AD 592 res.push_back(makevecteur(A,D,F,C)); 593 // Face 3 A B // D G=A+AB+AD 594 res.push_back(makevecteur(A,B,G,D)); 595 // Face 4 B E // G H 596 res.push_back(makevecteur(B,E,H,G)); 597 // Face 5 C E // F H 598 res.push_back(makevecteur(C,F,H,E)); 599 return polyedre_face(res,attributs,contextptr); 600 } _parallelepipede(const gen & args,GIAC_CONTEXT)601 gen _parallelepipede(const gen & args,GIAC_CONTEXT){ 602 if ( args.type==_STRNG && args.subtype==-1) return args; 603 if (args.type!=_VECT) 604 return gensizeerr(contextptr); 605 vecteur v(*args._VECTptr); 606 vecteur attributs(1,default_color(contextptr)); 607 int s=read_attributs(v,attributs,contextptr); 608 if (s==3){ 609 v.insert(v.begin(),makevecteur(0,0,0)); 610 ++s; 611 } 612 if (s!=4) 613 return gendimerr(contextptr); 614 gen A=remove_at_pnt(v[0]); 615 gen B=remove_at_pnt(v[1]); 616 if (B.type==_VECT && B.subtype==_VECTOR__VECT && B._VECTptr->size()==2) 617 B=A+B._VECTptr->back()-B._VECTptr->front(); 618 gen C=remove_at_pnt(v[2]); 619 if (C.type==_VECT && C.subtype==_VECTOR__VECT && C._VECTptr->size()==2) 620 C=A+C._VECTptr->back()-C._VECTptr->front(); 621 gen D=remove_at_pnt(v[3]); 622 if (D.type==_VECT && D.subtype==_VECTOR__VECT && D._VECTptr->size()==2) 623 D=A+D._VECTptr->back()-D._VECTptr->front(); 624 return parallelepipede4(A,B,C,D,attributs,contextptr); 625 } 626 static const char _parallelepipede_s []="parallelepiped"; 627 static define_unary_function_eval (__parallelepipede,&_parallelepipede,_parallelepipede_s); 628 define_unary_function_ptr5( at_parallelepipede ,alias_at_parallelepipede,&__parallelepipede,0,true); 629 pyramide4(const gen & A0,const gen & B0,const gen & C0,const gen & D0,const vecteur & attributs,GIAC_CONTEXT)630 static gen pyramide4(const gen & A0,const gen & B0,const gen & C0,const gen & D0,const vecteur & attributs,GIAC_CONTEXT){ 631 vecteur res; 632 gen A(A0),B(B0),C(C0),D(D0); 633 A.subtype=_POINT__VECT; 634 B.subtype=_POINT__VECT; 635 C.subtype=_POINT__VECT; 636 D.subtype=_POINT__VECT; 637 // Face 1 A B C 638 res.push_back(makevecteur(A,B,C)); 639 // Face 2 A C D 640 res.push_back(makevecteur(A,C,D)); 641 // Face 3 A B D 642 res.push_back(makevecteur(A,B,D)); 643 // Face 4 B C D 644 res.push_back(makevecteur(B,C,D)); 645 return polyedre_face(res,attributs,contextptr); 646 } 647 _pyramide(const gen & args,GIAC_CONTEXT)648 gen _pyramide(const gen & args,GIAC_CONTEXT){ 649 if ( args.type==_STRNG && args.subtype==-1) return args; 650 if (args.type!=_VECT) 651 return gensizeerr(contextptr); 652 vecteur v(*args._VECTptr); 653 vecteur attributs(1,default_color(contextptr)); 654 int s=read_attributs(v,attributs,contextptr); 655 if (s<2) 656 return gendimerr(contextptr); 657 v=vecteur(v.begin(),v.begin()+s); 658 gen A=remove_at_pnt(v[0]); 659 if (s==2){ 660 gen r=abs(v[1],contextptr); 661 v[1]=A+r*gen(makevecteur(1,0,0)); 662 v.push_back(A+r*gen(makevecteur(plus_one_half,plus_sqrt3_2,0))); 663 ++s; 664 } 665 gen B=remove_at_pnt(v[1]); 666 gen C=remove_at_pnt(v[2]); 667 if (s==3){ // tetraedre 668 gen AB=B-A,AC=C-A; 669 if (AB.type!=_VECT || AB._VECTptr->size()!=3 || AC.type!=_VECT || AC._VECTptr->size()!=3) 670 return gensizeerr(contextptr); 671 vecteur v1(*AB._VECTptr),v2(*AC._VECTptr); 672 vecteur n(cross(v1,v2,contextptr)); 673 v2=cross(n,v1,contextptr); 674 // Normalize 675 gen a(dotvecteur(v1,v1)); 676 v2=multvecteur(sqrt(3*a/dotvecteur(v2,v2),contextptr),v2); 677 C = A + divvecteur(v1,2) + divvecteur(v2,2) ; 678 n= multvecteur(sqrt(2*a/3/dotvecteur(n,n),contextptr),n); 679 gen D = A + divvecteur(v1,2) + divvecteur(v2,6) + n; 680 return pyramide4(A,B,C,D,attributs,contextptr); 681 } 682 gen D=remove_at_pnt(v[3]); 683 return pyramide4(A,B,C,D,attributs,contextptr); 684 } 685 static const char _pyramide_s []="pyramid"; 686 static define_unary_function_eval (__pyramide,&_pyramide,_pyramide_s); 687 define_unary_function_ptr5( at_pyramide ,alias_at_pyramide,&__pyramide,0,true); 688 689 static const char _tetraedre_s []="tetrahedron"; 690 static define_unary_function_eval (__tetraedre,&_pyramide,_tetraedre_s); 691 define_unary_function_ptr5( at_tetraedre ,alias_at_tetraedre,&__tetraedre,0,true); 692 693 // Find A,B,C,D such that AB=AC=AD and all are orthogonal cube_octaedre(const gen & args,gen & A,gen & B,gen & C,gen & D,vecteur & attributs,GIAC_CONTEXT)694 static bool cube_octaedre(const gen & args,gen & A,gen & B,gen & C,gen & D,vecteur & attributs,GIAC_CONTEXT){ 695 if (args.type!=_VECT) 696 return false; // gensizeerr(contextptr); 697 vecteur &v(*args._VECTptr); 698 int s=read_attributs(v,attributs,contextptr); 699 if (s<2) 700 return false; // gendimerr(contextptr); 701 A=v[0]; 702 B=v[1]; 703 if (s==2){ 704 gen r=abs(B,contextptr); 705 B=A+r*gen(makevecteur(r,0,0)); 706 C=A+r*gen(makevecteur(0,r,0)); 707 } 708 else 709 C=v[2]; 710 gen AB(B-A),AC(C-A); 711 if (AB.type!=_VECT || AB._VECTptr->size()!=3 || AC.type!=_VECT || AC._VECTptr->size()!=3) 712 return false; // gensizeerr(contextptr); 713 // AB cross AC -> normal direction to ABC plan 714 gen AB2(normal(scalar_product(AB,AB,contextptr),contextptr)); 715 if (is_undef(AB2)) 716 return false; 717 vecteur AD=*normal(cross(*AB._VECTptr,*AC._VECTptr,contextptr),contextptr)._VECTptr; 718 D=A+AD*sqrt(normal(AB2/dotvecteur(AD,AD),contextptr),contextptr); 719 // binormal direction gives the 2nd direction AB, AE, AD 720 vecteur AE=*normal(cross(AD,*AB._VECTptr,contextptr),contextptr)._VECTptr; 721 C=A+AE*sqrt(normal(AB2/dotvecteur(AE,AE),contextptr),contextptr); 722 A.subtype=_POINT__VECT; 723 B.subtype=_POINT__VECT; 724 C.subtype=_POINT__VECT; 725 D.subtype=_POINT__VECT; 726 return true; 727 } 728 _tetraedre_centre(const gen & args,GIAC_CONTEXT)729 gen _tetraedre_centre(const gen & args,GIAC_CONTEXT){ 730 if ( args.type==_STRNG && args.subtype==-1) return args; 731 gen O,b,c,d; 732 vecteur attributs(1,default_color(contextptr)); 733 if (!cube_octaedre(args,O,b,c,d,attributs,contextptr)) 734 return gensizeerr(contextptr); 735 gen v1(normal(b-O,contextptr)),v2(normal(c-O,contextptr)),v3(normal(d-O,contextptr)); 736 gen A = b; // O + v1 737 gen B = normal(O - v1/3 - sqrt(2,contextptr)*v2/3 - sqrt(6,contextptr)*v3/3,contextptr); 738 gen C = normal(O - v1/3 - sqrt(2,contextptr)*v2/3 + sqrt(6,contextptr)*v3/3,contextptr); 739 gen D = normal(O - v1/3 + 2*sqrt(2,contextptr)*v2/3,contextptr); 740 return pyramide4(A,B,C,D,attributs,contextptr); 741 } 742 static const char _tetraedre_centre_s []="centered_tetrahedron"; 743 static define_unary_function_eval (__tetraedre_centre,&_tetraedre_centre,_tetraedre_centre_s); 744 define_unary_function_ptr5( at_tetraedre_centre ,alias_at_tetraedre_centre,&__tetraedre_centre,0,true); 745 746 // args= 3 points A, B, C _cube(const gen & args,GIAC_CONTEXT)747 gen _cube(const gen & args,GIAC_CONTEXT){ 748 if ( args.type==_STRNG && args.subtype==-1) return args; 749 gen A,B,C,D; 750 vecteur attributs(1,default_color(contextptr)); 751 if (!cube_octaedre(args,A,B,C,D,attributs,contextptr)) 752 return gensizeerr(contextptr); 753 return parallelepipede4(A,B,C,D,attributs,contextptr); 754 } 755 static const char _cube_s []="cube"; 756 static define_unary_function_eval (__cube,&_cube,_cube_s); 757 define_unary_function_ptr5( at_cube ,alias_at_cube,&__cube,0,true); 758 759 // args= center A, vertex B, point C such that ABC is a symmetry plan _cube_centre(const gen & args,GIAC_CONTEXT)760 gen _cube_centre(const gen & args,GIAC_CONTEXT){ 761 if ( args.type==_STRNG && args.subtype==-1) return args; 762 if (args.type!=_VECT || args._VECTptr->size()<3) 763 return gensizeerr(contextptr); 764 gen A,B,C,D; 765 vecteur attributs(1,default_color(contextptr)); 766 if (!cube_octaedre(args,A,B,C,D,attributs,contextptr)) 767 return gensizeerr(contextptr); 768 gen x = (B-A)/3; 769 // Take C1=A+cos(theta)*(B-A)+sin(theta)*(C-A), cos(theta)=1/3 770 gen C1 = normal(A+x+2*sqrt(2,contextptr)/3*(C-A),contextptr); 771 gen D1 = normal(A+x-sqrt(2,contextptr)/3*(C-A)+sqrt(6,contextptr)/3*(D-A),contextptr); 772 if (!cube_octaedre(makevecteur(B,D1,C1),A,B,C,D,attributs,contextptr)) 773 return gensizeerr(contextptr); 774 return parallelepipede4(A,B,C,D,attributs,contextptr); 775 } 776 static const char _cube_centre_s []="centered_cube"; 777 static define_unary_function_eval (__cube_centre,&_cube_centre,_cube_centre_s); 778 define_unary_function_ptr5( at_cube_centre ,alias_at_cube_centre,&__cube_centre,0,true); 779 780 // args= center A, point B, C such that ABC = symmetry plan with 4 vertices _octaedre(const gen & args,GIAC_CONTEXT)781 gen _octaedre(const gen & args,GIAC_CONTEXT){ 782 if ( args.type==_STRNG && args.subtype==-1) return args; 783 gen A,B,C,D; 784 vecteur attributs(1,default_color(contextptr)); 785 if (!cube_octaedre(args,A,B,C,D,attributs,contextptr)) 786 return gensizeerr(contextptr); 787 // B, C, D are 3 vertices 788 gen E,F,G; 789 E = A - (B-A); 790 F = A - (C-A); 791 G = A - (D-A); 792 vecteur faces; 793 faces.push_back(makevecteur(B,C,D)); 794 faces.push_back(makevecteur(B,C,G)); 795 faces.push_back(makevecteur(B,F,D)); 796 faces.push_back(makevecteur(B,F,G)); 797 faces.push_back(makevecteur(E,C,D)); 798 faces.push_back(makevecteur(E,C,G)); 799 faces.push_back(makevecteur(E,F,D)); 800 faces.push_back(makevecteur(E,F,G)); 801 return polyedre_face(faces,attributs,contextptr); 802 } 803 static const char _octaedre_s []="octahedron"; 804 static define_unary_function_eval (__octaedre,&_octaedre,_octaedre_s); 805 define_unary_function_ptr5( at_octaedre ,alias_at_octaedre,&__octaedre,0,true); 806 res_push(vecteur & res,gen * s,int i,int j,int k)807 static void res_push(vecteur & res,gen * s, int i,int j,int k){ 808 res.push_back(makevecteur(s[i],s[j],s[k])); 809 } 810 // args= centre, sommet1, sommet2 (direction of) _icosaedre(const gen & args,GIAC_CONTEXT)811 gen _icosaedre(const gen & args,GIAC_CONTEXT){ 812 if ( args.type==_STRNG && args.subtype==-1) return args; 813 if (args.type!=_VECT) 814 return gensizeerr(contextptr); 815 gen errcode=checkanglemode(contextptr); 816 if (is_undef(errcode)) return errcode; 817 vecteur & v = *args._VECTptr; 818 vecteur attributs(1,default_color(contextptr)); 819 int sv=read_attributs(v,attributs,contextptr); 820 if (sv!=3) 821 return gendimerr(contextptr); 822 gen s[12]; 823 gen centre=v[0],s1=v[1],s2=v[2]; 824 gen v1g(s1-centre),v2g(s2-centre); 825 // Icosaedre=s1+symetric of s1 with respect to center 826 // + 2* 5 points as a pentagon on 2 plans perpendicular to centre->s1 827 // If the distance of the 2 // plans is 1 to the center 828 // Then the 5 vertices are at distance 2 to the intersection axe/plan 829 // (abscisse=+/-1, sqrt(y^2+z^2)=2) 830 // and |s1-centre|=sqrt(5) 831 if (v1g.type!=_VECT|| v2g.type!=_VECT) 832 return gensizeerr(contextptr); 833 vecteur v1(*v1g._VECTptr),v2(*v2g._VECTptr); 834 vecteur n(cross(v1,v2,contextptr)); 835 v2=divvecteur(cross(n,v1,contextptr),sqrt(dotvecteur(n,n),contextptr)); 836 // norm=distance(centre,sommet) 837 n=multvecteur(sqrt(dotvecteur(v1,v1)/dotvecteur(n,n),contextptr),n); 838 // centre +/- (v1/sqrt(5) + 2/sqrt(5)*(cos(2*k*pi/5)*v2 +sin(2*k*pi/5)*n)) 839 s[0]=s1; 840 s[11]=s1-multvecteur(2,v1); 841 for (int i=0;i<5;++i){ 842 context ctmp; 843 gen tmp = gen(1)/sqrt(5,contextptr)*(gen(v1) + 2*(cos(2*i*cst_pi/5,&ctmp)*gen(v2)+sin(2*i*cst_pi/5,&ctmp)*n)); 844 s[1+i] = centre + tmp; 845 s[10-i] = centre - tmp; 846 } 847 // Make 5 faces s[0] with s[1+i], s[2+i] for i in 1..4 and with 5,1 848 // 5 with s[11] with s[11-i], s[10-i] for i in 1..4 and with 849 // 1,7,8 + 1,7,2 + 2,6,7 + 2,6,3 + 3,10,6 + 3,10,4 + 4,9,10 + 4,9,5 + 5,8,9+ 5,8,1 850 vecteur res; 851 for (int i=1;i<5;++i){ 852 res_push(res,s,0,i,1+i); res_push(res,s,11,11-i,10-i); 853 } 854 res_push(res,s,0,5,1); res_push(res,s,11,6,10); 855 res_push(res,s,1,7,8); res_push(res,s,1,7,2); 856 res_push(res,s,2,6,7); res_push(res,s,2,6,3); 857 res_push(res,s,3,10,6); res_push(res,s,3,10,4); 858 res_push(res,s,4,9,10); res_push(res,s,4,9,5); 859 res_push(res,s,5,8,9); res_push(res,s,5,8,1); 860 return polyedre_face(res,attributs,contextptr); 861 } 862 static const char _icosaedre_s []="icosahedron"; 863 static define_unary_function_eval (__icosaedre,&_icosaedre,_icosaedre_s); 864 define_unary_function_ptr5( at_icosaedre ,alias_at_icosaedre,&__icosaedre,0,true); 865 res_push(vecteur & res,gen * s,int i,int j,int k,int l,int m)866 static void res_push(vecteur & res,gen * s, int i,int j,int k,int l,int m){ 867 res.push_back(makevecteur(s[i],s[j],s[k],s[l],s[m])); 868 } 869 // args= centre, sommet1, 3rd point defining a plan containing the axis 870 // Example dodecaedre([0,0,0],[0,2,sqrt(5)/2+3/2],[0,0,1]) _dodecaedre(const gen & args,GIAC_CONTEXT)871 gen _dodecaedre(const gen & args,GIAC_CONTEXT){ 872 if ( args.type==_STRNG && args.subtype==-1) return args; 873 if (args.type!=_VECT) 874 return gensizeerr(contextptr); 875 vecteur attributs(1,default_color(contextptr)); 876 int sv=read_attributs(*args._VECTptr,attributs,contextptr); 877 if (sv!=3) 878 return gendimerr(contextptr); 879 gen errcode=checkanglemode(contextptr); 880 if (is_undef(errcode)) return errcode; 881 vecteur v = *evalf(args,1,contextptr)._VECTptr; 882 gen centre=v[0],s1=v[1],s2=v[2]; 883 gen v1g(s1-centre),v2g(s2-centre); 884 if (v1g.type!=_VECT|| v2g.type!=_VECT) 885 return gensizeerr(contextptr); 886 vecteur v1(*v1g._VECTptr),v2(*v2g._VECTptr); 887 gen phi=evalf((sqrt(5,contextptr)+1)/2,1,0); 888 gen r2(dotvecteur(v1,v1)); // r = |v1| = sqrt(6+3phi)*unit of length 889 vecteur w2(cross(v1,v2,contextptr)); // y direction 890 vecteur v3(cross(w2,v1,contextptr)); // v1,v3 contains the axis v1.v3=0 891 // Now normalize w2 to 1 unit of length 892 w2=multvecteur(sqrt(r2/dotvecteur(w2,w2)/(6+3*phi),contextptr),w2); 893 v3=multvecteur(sqrt(r2/dotvecteur(v3,v3),contextptr),v3); // |v3|=|v1|=r 894 gen w1=(2*gen(v1)-(phi+1)*gen(v3))/(6+3*phi); // |w1|=sqrt(6+3phi)*norm(v1)/(6+3phi) 895 gen w3=((phi+1)*v1+2*gen(v3))/(6+3*phi); // = |w2|=1 unit of length 896 // Dodecaedre at center 0. Edge length=sqrt(10-2*sqrt(5)) 897 // Golden ratio phi=(sqrt(5)+1)/2 (phi^2=phi+1) 898 // Vertices at +/-(2*cos(2*pi/5),2*sin(2*pi/5),phi+1) 899 // and +/-(2*phi*cos(2*pi/5),2*phi*sin(2*pi/5),phi-1) 900 // Sphere radius ^2 = 4 + (phi+1)^2 = 6+3*phi=(15+3*sqrt(5))/2 901 gen s[20]; 902 context ctmp; 903 for (int i=0;i<5;++i){ 904 s[i]=centre+evalf(2*cos(2*i*cst_pi/5,&ctmp)*w1+2*sin(2*i*cst_pi/5,&ctmp)*w2+(phi+1)*w3,1,contextptr); 905 s[15+i]=centre-s[i]; 906 s[5+i]=centre+evalf(2*phi*(cos(2*i*cst_pi/5,&ctmp)*w1+sin(2*i*cst_pi/5,&ctmp)*w2)+(phi-1)*w3,1,contextptr); 907 s[10+i]=centre-s[5+i]; 908 } 909 vecteur res; 910 res_push(res,s,0,1,2,3,4); res_push(res,s,15,16,17,18,19); 911 for (int i=0;i<5;++i){ 912 res_push(res,s,(i+1)%5,5+(i+1)%5,10+(3+i)%5,5+i,i); 913 res_push(res,s,15+(i+1)%5,10+(i+1)%5,5+(3+i)%5,10+i,15+i); 914 } 915 return polyedre_face(res,attributs,contextptr); 916 } 917 static const char _dodecaedre_s []="dodecahedron"; 918 static define_unary_function_eval (__dodecaedre,&_dodecaedre,_dodecaedre_s); 919 define_unary_function_ptr5( at_dodecaedre ,alias_at_dodecaedre,&__dodecaedre,0,true); 920 _aretes(const gen & args,GIAC_CONTEXT)921 gen _aretes(const gen & args,GIAC_CONTEXT){ 922 if ( args.type==_STRNG && args.subtype==-1) return args; 923 bool tmp=show_point(contextptr); 924 show_point(false,contextptr); 925 gen g=remove_at_pnt(args); 926 vecteur v(gen2vecteur(g)); 927 vecteur res; 928 const_iterateur it=v.begin(),itend=v.end(); 929 for (;it!=itend;++it){ 930 if (!ckmatrix(*it)) 931 return gensizeerr(contextptr); 932 const_iterateur jt=it->_VECTptr->begin(),jtend=it->_VECTptr->end(); 933 for (;jt+1!=jtend;++jt){ 934 res.push_back(_segment(makesequence(*jt,*(jt+1)),contextptr)); 935 } 936 res.push_back(_segment(makesequence(*jt,it->_VECTptr->front()),contextptr)); 937 } 938 show_point(tmp,contextptr); 939 return res; 940 } 941 static const char _aretes_s []="line_segments"; 942 static define_unary_function_eval (__aretes,&_aretes,_aretes_s); 943 define_unary_function_ptr5( at_aretes ,alias_at_aretes,&__aretes,0,true); 944 rotation3d(const gen & elem,const gen & b,GIAC_CONTEXT)945 static gen rotation3d(const gen & elem,const gen & b,GIAC_CONTEXT){ 946 if (elem.type==_VECT && elem._VECTptr->size()==2){ 947 gen A=elem._VECTptr->front(); 948 gen M=elem._VECTptr->back(); 949 gen res=A+M*(b-A); 950 res.subtype=_POINT__VECT; 951 return res; 952 } 953 return gensizeerr(contextptr); 954 } 955 similitude3d(const vecteur & centrev,const gen & angle,const gen & rapport,const gen & b,int symrot,GIAC_CONTEXT)956 gen similitude3d(const vecteur & centrev,const gen & angle,const gen & rapport,const gen & b,int symrot,GIAC_CONTEXT){ 957 if (centrev.size()!=2 || centrev.front().type!=_VECT || centrev.back().type!=_VECT) 958 return gensizeerr(contextptr); 959 vecteur A(*centrev.front()._VECTptr),B(*centrev.back()._VECTptr); 960 vecteur AB(subvecteur(B,A)); 961 if (AB.size()!=3) 962 return gendimerr(contextptr); 963 // Find rotation matrix from isom.cc/h 964 gen M=rapport*mkisom(makevecteur(AB,angle),symrot,contextptr); 965 gen elem(makevecteur(A,M)); 966 if (b.type==_VECT) 967 return symb_pnt(apply3d(elem,b,contextptr,rotation3d),default_color(contextptr),contextptr); 968 if (b.is_symb_of_sommet(at_hypersphere)){ 969 gen c,r; 970 centre_rayon(b,c,r,false,contextptr); 971 c=A+M*(c-A); 972 return _sphere(makesequence(c,r),contextptr); 973 } 974 if (b.is_symb_of_sommet(at_hyperplan)){ 975 vecteur n,P; 976 if (!hyperplan_normal_point(b,n,P)) 977 return gensizeerr(contextptr); 978 gen Pr=A+M*(P-A); 979 gen nr=M*n; 980 return _plan(makesequence(nr,Pr),contextptr); 981 } 982 return curve_surface_apply(elem,b,rotation3d,contextptr); 983 } 984 985 plotparam3d(const gen & f,const vecteur & vars,double function_xmin,double function_xmax,double function_ymin,double function_ymax,double function_zmin,double function_zmax,double function_umin,double function_umax,double function_vmin,double function_vmax,bool densityplot,bool f_autoscale,const vecteur & attributs,double ustep,double vstep,const gen & eq,const vecteur & eqvars,GIAC_CONTEXT)986 gen plotparam3d(const gen & f,const vecteur & vars,double function_xmin,double function_xmax, double function_ymin, double function_ymax,double function_zmin,double function_zmax,double function_umin,double function_umax,double function_vmin,double function_vmax,bool densityplot,bool f_autoscale,const vecteur & attributs,double ustep,double vstep,const gen & eq,const vecteur & eqvars,GIAC_CONTEXT){ 987 int color=default_color(contextptr); 988 gen attribut=attributs.empty()?color:attributs[0]; 989 if (attribut.type==_INT_) 990 color=attribut.val; 991 identificateur u("u"),v("v"); 992 vecteur vsub(1,u); 993 vsub.push_back(v); 994 gen ff=subst(f,vars,vsub,false,contextptr); 995 gen r=symbolic(at_plotparam,makesequence(f,vars)); 996 if (is_zero(derive(ff,v,contextptr))){ 997 vecteur res; 998 double x=function_umin; 999 double dx=ustep; 1000 int n=int((function_umax-function_umin)/ustep+.5); 1001 vecteur values; 1002 gen prec,cur,tmp; 1003 double Dx=function_xmax-function_xmin,Dy=function_ymax-function_ymin,Dz=function_zmax-function_zmin,fmin=function_umin; 1004 for (int i=0;i<=n;++i,x+=dx){ 1005 cur=evalf_double(subst(f,vars[0],x,false,contextptr),1,contextptr); 1006 if (cur.type!=_VECT || cur._VECTptr->size()!=3) 1007 continue; 1008 cur.subtype=_POINT__VECT; 1009 vecteur & curt=*cur._VECTptr; 1010 if (curt[0].type!=_DOUBLE_ || curt[1].type!=_DOUBLE_ || curt[2].type!=_DOUBLE_) 1011 continue; 1012 bool joindre=true; 1013 if (prec.type==_VECT && prec._VECTptr->size()==3 && !values.empty()){ 1014 vecteur & precv=*prec._VECTptr; 1015 double oldx=precv[0]._DOUBLE_val,oldy=precv[1]._DOUBLE_val,oldz=precv[2]._DOUBLE_val; 1016 double curx=curt[0]._DOUBLE_val,cury=curt[1]._DOUBLE_val,curz=curt[2]._DOUBLE_val; 1017 if (absdouble(curx-oldx)>Dx/10 || absdouble(cury-oldy)>Dy/10 || absdouble(curz-oldz)>Dz/10 ){ 1018 tmp=evalf_double(subst(f,vars[0],x+dx/2,false,contextptr),1,contextptr); 1019 if (tmp.type!=_VECT || tmp._VECTptr->size()!=3) 1020 continue; 1021 tmp.subtype=_POINT__VECT; 1022 vecteur & tmpv=*tmp._VECTptr; 1023 double entrex=tmpv[0]._DOUBLE_val,entrey=tmpv[1]._DOUBLE_val,entrez=tmpv[2]._DOUBLE_val; 1024 if ( (entrex-oldx)*(curx-entrex)<0 || (entrey-oldy)*(cury-entrey)<0 || (entrez-oldz)*(curz-entrez)<0 ) 1025 joindre=false; 1026 } 1027 } 1028 if (joindre) 1029 values.push_back(cur); 1030 else { 1031 res.push_back(symb_pnt(symb_curve(gen(makevecteur(f,vars[0],fmin,x,gen(values,_GROUP__VECT)),_PNT__VECT),undef),color,contextptr)); 1032 fmin=x; 1033 values.clear(); 1034 prec=0; 1035 } 1036 prec=cur; 1037 } 1038 r=symb_pnt(symb_curve(gen(makevecteur(f,vars[0],fmin,function_umax,gen(values,_GROUP__VECT)),_PNT__VECT),undef),color,contextptr); 1039 if (!res.empty()){ 1040 res.push_back(r); 1041 r=res; // gen(res,_SEQ__VECT); 1042 } 1043 } // if is_zero(derive(ff,v)) 1044 else { 1045 vecteur vals(2); 1046 double x=function_umin,y=function_vmin; 1047 double dx=ustep; 1048 double dy=vstep; 1049 int nu=int((function_umax-function_umin)/ustep+.5),nv=int((function_vmax-function_vmin)/vstep+.5); 1050 // Compute a grid of values 1051 vecteur values; 1052 for (int i=0;i<=nu;++i,x+=dx){ 1053 y=function_vmin; 1054 vals[0]=x; 1055 vecteur tmp; 1056 for (int j=0;j<=nv;++j,y+=dy){ 1057 vals[1]=y; 1058 gen tmppnt=evalf_double(subst(f,vars,vals,false,contextptr),1,contextptr); 1059 if (tmppnt.type==_VECT) 1060 tmppnt.subtype=_POINT__VECT; 1061 tmp.push_back(tmppnt); 1062 } 1063 values.push_back(gen(tmp,_GROUP__VECT)); 1064 } 1065 r=symb_pnt(hypersurface(gen(makevecteur(f,vars,makevecteur(function_umin,function_vmin),makevecteur(function_umax,function_vmax),gen(values,_GROUP__VECT)),_PNT__VECT),eq,eqvars),color,contextptr); 1066 } 1067 #ifdef WITH_GNUPLOT 1068 int out_handle; 1069 bool clrplot=false; 1070 FILE * gnuplot_out_readstream,* stream = open_gnuplot(clrplot,gnuplot_out_readstream,out_handle); 1071 #ifdef IPAQ 1072 fprintf(stream,"set samples 10\n"); 1073 // fprintf(stream,"show samples\n"); 1074 #endif 1075 r.subtype=gnuplot_fileno; 1076 reset_gnuplot_hidden3d(stream); 1077 if (debug_infolevel) 1078 printf("set urange [%g:%g]\n",function_umin,function_umax); 1079 fprintf(stream,"set urange [%g:%g]\n",function_umin,function_umax); 1080 if (debug_infolevel) 1081 printf("set vrange [%g:%g]\n",function_vmin,function_vmax); 1082 fprintf(stream,"set vrange [%g:%g]\n",function_vmin,function_vmax); 1083 if (!f_autoscale){ 1084 if (debug_infolevel) 1085 printf("set xrange [%g:%g]\n",function_xmin,function_xmax); 1086 fprintf(stream,"set xrange [%g:%g]\n",function_xmin,function_xmax); 1087 if (debug_infolevel) 1088 printf("set yrange [%g:%g]\n",function_ymin,function_ymax); 1089 fprintf(stream,"set yrange [%g:%g]\n",function_ymin,function_ymax); 1090 if (debug_infolevel) 1091 printf("set zrange [%g:%g]\n",function_zmin,function_zmax); 1092 fprintf(stream,"set zrange [%g:%g]\n",function_zmin,function_zmax); 1093 } 1094 if (clrplot || gnuplot_do_splot){ 1095 if (debug_infolevel) 1096 printf("%s","splot "); 1097 fprintf(stream,"%s","splot "); 1098 } 1099 else { 1100 if (debug_infolevel) 1101 printf("%s","replot "); 1102 fprintf(stream,"%s","replot "); 1103 } 1104 gnuplot_do_splot=false; 1105 if (ff.type!=_VECT){ 1106 if (debug_infolevel) 1107 printf("%s notitle\n",gnuplot_traduit(ff).c_str()); 1108 fprintf(stream,"%s notitle\n",gnuplot_traduit(ff).c_str()); 1109 } 1110 else { 1111 string tmp(gnuplot_traduit(ff)); 1112 // cout << tmp.substr(1,tmp.size()-2) << endl; 1113 if (debug_infolevel) 1114 printf("%s notitle\n",tmp.substr(1,tmp.size()-2).c_str()); 1115 fprintf(stream,"%s notitle\n",tmp.substr(1,tmp.size()-2).c_str()); 1116 } 1117 win9x_gnuplot(stream); 1118 gnuplot_wait(out_handle,gnuplot_out_readstream,gnuplot_wait_times); 1119 ++gnuplot_fileno; 1120 // return gnuplot_fileno-1; 1121 return r; 1122 #endif // GNUPLOT 1123 return r; 1124 } 1125 normal3d(const gen & nn,vecteur & v1,vecteur & v2)1126 bool normal3d(const gen & nn,vecteur & v1,vecteur & v2){ 1127 if (nn.type!=_VECT || nn._VECTptr->size()!=3) 1128 return false; 1129 vecteur & n = *nn._VECTptr; 1130 if (is_zero(n[0])) 1131 v1=makevecteur(1,0,0); 1132 else 1133 v1=makevecteur(n[1],-n[0],0); 1134 v2=cross(n,v1,context0); 1135 return true; 1136 } 1137 hypersphere_equation(const gen & g,const vecteur & xyz)1138 gen hypersphere_equation(const gen & g,const vecteur & xyz){ 1139 gen centre,rayon; 1140 if (!centre_rayon(g,centre,rayon,false,0) ||centre.type!=_VECT) 1141 return gensizeerr(gettext("hypersphere_equation")); 1142 vecteur & v=*centre._VECTptr; 1143 if (v.size()!=3) 1144 return gendimerr(gettext("hypersphere_equation")); 1145 vecteur xyzc(subvecteur(xyz,v)); 1146 gen eq=ratnormal(dotvecteur(xyzc,xyzc)-pow(rayon,2),context0); 1147 return eq; 1148 } hypersphere_parameq(const gen & g,const vecteur & st)1149 vecteur hypersphere_parameq(const gen & g,const vecteur & st){ 1150 gen centre,rayon; 1151 if (!centre_rayon(g,centre,rayon,false,0) ||centre.type!=_VECT) 1152 return vecteur(1,gensizeerr(gettext("hypersphere_parameq"))); 1153 vecteur & v=*centre._VECTptr; 1154 if (v.size()!=3) 1155 return vecteur(1,gendimerr(gettext("hypersphere_parameq"))); 1156 vecteur res(4); 1157 res[0]=centre+makevecteur(rayon*symb_cos(st[0])*symb_cos(st[1]),rayon*symb_cos(st[0])*symb_sin(st[1]),rayon*symb_sin(st[0])); 1158 res[1]=st; 1159 res[2]=makevecteur(-cst_pi_over_2,0); 1160 res[3]=makevecteur(cst_pi_over_2,cst_two_pi); 1161 return res; 1162 } 1163 hypersurface_equation(const gen & g,const vecteur & xyz,GIAC_CONTEXT)1164 gen hypersurface_equation(const gen & g,const vecteur & xyz,GIAC_CONTEXT){ 1165 if (!g.is_symb_of_sommet(at_hypersurface)) 1166 return gensizeerr(contextptr); 1167 gen & f=g._SYMBptr->feuille; 1168 if (f.type!=_VECT) return gensizeerr(contextptr); 1169 vecteur & fv=*f._VECTptr; 1170 if (fv.size()==3 && fv[1].type!=_VECT && fv[2].type==_VECT){ 1171 gen eq=fv[1]; 1172 if (is_undef(eq)){ 1173 gen point=fv[0]; 1174 if (point.type==_VECT && point._VECTptr->size()>=2){ 1175 gen f=(*point._VECTptr)[0]; 1176 gen vars=(*point._VECTptr)[1]; 1177 if (vars.type==_VECT && vars._VECTptr->size()==2 && f.type==_VECT && f._VECTptr->size()==3 && xyz.size()==3){ 1178 vecteur lv(*vars._VECTptr); 1179 lvar(f,lv); 1180 if (lv==vars){ // resultant -> eq 1181 gen tmp1=_resultant(makesequence(f[0]-xyz[0],f[1]-xyz[1],vars[0]),contextptr); 1182 if (is_undef(tmp1)) 1183 return tmp1; 1184 gen tmp2=_resultant(makesequence(f[0]-xyz[0],f[2]-xyz[2],vars[0]),contextptr); 1185 if (is_undef(tmp2)) 1186 return tmp2; 1187 gen res=_resultant(makesequence(tmp1,tmp2,vars[1]),contextptr); 1188 return res; 1189 } 1190 } 1191 } 1192 } 1193 return subst(fv[1],*fv[2]._VECTptr,xyz,false,contextptr); 1194 } 1195 return gensizeerr(gettext("Hypersurface w/o equation")); 1196 } 1197 1198 // a must be a vector of length 2, b a symbolic interdroitehyperplan(const gen & a,const gen & b,GIAC_CONTEXT)1199 vecteur interdroitehyperplan(const gen & a,const gen &b,GIAC_CONTEXT){ 1200 if (a.type!=_VECT || b.type!=_SYMB || a._VECTptr->size()!=2) 1201 return vecteur(1,gensizeerr(contextptr)); 1202 // D inter H 1203 gen A=a._VECTptr->front(),B=a._VECTptr->back(); // D=(AB) 1204 gen & f=b._SYMBptr->feuille; 1205 gen AB=B-A; 1206 if (f.type!=_VECT || f._VECTptr->size()!=2 ) 1207 return vecteur(1,gensizeerr(contextptr)); 1208 gen C=f._VECTptr->back(),Hn=f._VECTptr->front(); // H= C normal is n 1209 gen AC=C-A; 1210 if (Hn.type!=_VECT || AB.type!=_VECT || AC.type!=_VECT) 1211 return vecteur(1,gensizeerr(contextptr)); 1212 vecteur v(*AB._VECTptr),n(*Hn._VECTptr); 1213 gen vn(normal(dotvecteur(v,n),contextptr)); 1214 if (is_zero(vn)){ // D is parallel to H 1215 return vecteur(0); // FIXME should be D if D is in H 1216 } 1217 // H inter D = A + t*v, with t such that A+t*v-C is normal to n 1218 // Hence t=(C-A).n/(v.n) 1219 gen t=dotvecteur(*AC._VECTptr,n)/vn; 1220 gen M(_point(A+t*v,contextptr)); 1221 return remove_not_in_segment(A,B,a.subtype,vecteur(1,M),contextptr); 1222 } 1223 1224 // a hyperplan, b hypersphere interplansphere(const gen & a,const gen & b,GIAC_CONTEXT)1225 vecteur interplansphere(const gen & a,const gen & b,GIAC_CONTEXT){ 1226 gen cg,r; 1227 if (!centre_rayon(b,cg,r,false,contextptr)) return vecteur(1,gensizeerr(contextptr)); 1228 if (cg.type!=_VECT || cg._VECTptr->size()!=3) 1229 return vecteur(1,gensizeerr(contextptr)); 1230 vecteur c(*cg._VECTptr); 1231 vecteur n,p,res; 1232 vecteur attributs(1,default_color(contextptr)); 1233 if (!hyperplan_normal_point(a,n,p)) 1234 return vecteur(1,gensizeerr(contextptr)); 1235 gen n2=dotvecteur(n,n),r2=ratnormal(r*r,contextptr); 1236 // find x such that c+x*n must belong to a, 1237 // hence (c-p+x*n).n=0 -> x=(p-c).n/n.n 1238 gen x=dotvecteur(subvecteur(p,c),n)/n2; 1239 // if ||x*n||=r, one point c+x*n 1240 // if ||x*n||>r, empty 1241 // if ||x*n||<r, circle of radius sqrt(xn2-r2), 1242 // centered at c+x*n inside the plan b 1243 gen xn2=ratnormal(x*x*n2,contextptr); 1244 gen center=c+x*n; 1245 identificateur id(" plansphere"); 1246 gen T__IDNT_e(id); 1247 if (xn2==r2) 1248 res.push_back(_point(center,contextptr)); 1249 else { 1250 if (is_strictly_greater(r2,xn2,contextptr)){ 1251 vecteur v1,v2; 1252 if (!normal3d(n,v1,v2)) 1253 return vecteur(1,gensizeerr(contextptr)); 1254 gen v12(dotvecteur(v1,v1)); 1255 gen v22(dotvecteur(v2,v2)); 1256 v12=sqrt((r2-xn2)/v12,contextptr); 1257 v22=sqrt((r2-xn2)/v22,contextptr); 1258 res.push_back(plotparam3d(center+cos(T__IDNT_e,contextptr)*v12*v1+sin(T__IDNT_e,contextptr)*v22*v2, 1259 makevecteur(T__IDNT_e,u__IDNT_e), 1260 gnuplot_xmin,gnuplot_xmax,gnuplot_ymin,gnuplot_ymax,gnuplot_zmin,gnuplot_zmax, 1261 0,2*M_PI,0,0,false,false,attributs,M_PI/30,0, 1262 undef /* FIXME: equation */,makevecteur(T__IDNT_e,u__IDNT_e),contextptr)); 1263 } 1264 } 1265 return res; 1266 } 1267 inter_solve(const gen & args,GIAC_CONTEXT)1268 static gen inter_solve(const gen & args,GIAC_CONTEXT){ 1269 bool b=all_trig_sol(contextptr); 1270 all_trig_sol(false,contextptr); 1271 gen res=_solve(args,contextptr); 1272 all_trig_sol(b,contextptr); 1273 return res; 1274 } 1275 1276 // a=hypersurface, b=hypersurface inter2hypersurface(const gen & a,const gen & b,GIAC_CONTEXT)1277 vecteur inter2hypersurface(const gen & a,const gen &b,GIAC_CONTEXT){ 1278 gen & af=a._SYMBptr->feuille; 1279 gen & bf=b._SYMBptr->feuille; 1280 if (af.type!=_VECT || bf.type!=_VECT || af._VECTptr->empty() || bf._VECTptr->empty() ) 1281 return vecteur(1,gensizeerr(contextptr)); 1282 vecteur av=*af._VECTptr,bv=*bf._VECTptr; 1283 bool aparam=av[0].type==_VECT,bparam=bv[0].type==_VECT; 1284 if (!aparam && bparam) 1285 return inter2hypersurface(b,a,contextptr); 1286 vecteur res; 1287 identificateur ids(" s"),idt(" t"),idu(" u"),idv(" v"); 1288 gen s(ids),t(idt),u(idu),v(idv); 1289 if (aparam ){ 1290 av=*av[0]._VECTptr; 1291 gen A=subst(av[0],av[1],makevecteur(s,t),false,contextptr); 1292 if (bparam && bv.size()<3){ 1293 bv=*bv[0]._VECTptr; 1294 // av[0]=point on hypersurface a, av[1]=parameters 1295 // bv[]= same on b 1296 // Rename parameters so that they do not have the same name for a and b 1297 gen B=subst(bv[0],bv[1],makevecteur(u,v),false,contextptr); 1298 if (A.type!=_VECT || A._VECTptr->size()!=3 || B.type!=_VECT || B._VECTptr->size()!=3 ) 1299 return vecteur(1,gensizeerr(contextptr)); 1300 // we have now to solve A[0]=B[0], A[1]=B[1], A[2]=B[2] 1301 // that should give us t,u,v as a function of s 1302 // then we will draw the parametric curve A with respect to s (t replaced) 1303 vecteur veq(makevecteur(A[0]-B[0],A[1]-B[1],A[2]-B[2])); 1304 gen sol=inter_solve(gen(makevecteur(veq,makevecteur(t,u,v)),_SEQ__VECT),contextptr); 1305 if (sol.type!=_VECT) 1306 return vecteur(1,gensizeerr(contextptr)); 1307 // for each element of sol, get the first component t, subst in A 1308 int nsol=int(sol._VECTptr->size()); 1309 for (int i=0;i<nsol;i++){ 1310 if (sol[i].type==_VECT && sol[i]._VECTptr->size()==3){ 1311 gen As=ratnormal(subst(A,t,sol[i]._VECTptr->front(),false,contextptr),contextptr); 1312 // now make the parametric curves [A,s] 1313 res.push_back(paramplotparam(gen(makevecteur(As,s),_SEQ__VECT),false,contextptr)); 1314 } 1315 } 1316 } // end both parametric 1317 else { 1318 if (bv.size()<3) 1319 return vecteur(1,gensizeerr(contextptr)); 1320 // a parametric, b by equation bv[1] parameters in bv[2] 1321 // find curve equation with respect to s and t 1322 gen curveeq=subst(bv[1],bv[2],A,false,contextptr); 1323 bool swapped=is_zero(derive(curveeq,t,contextptr)); 1324 if (swapped) 1325 std::swap(s,t); 1326 gen sol=inter_solve(gen(makevecteur(symbolic(at_equal,makesequence(curveeq,0)),t),_SEQ__VECT),contextptr); 1327 if (sol.type!=_VECT) 1328 return vecteur(1,gensizeerr(contextptr)); 1329 if (!swapped && sol._VECTptr->empty()){ 1330 std::swap(s,t); 1331 sol=inter_solve(gen(makevecteur(symbolic(at_equal,makesequence(curveeq,0)),t),_SEQ__VECT),contextptr); 1332 if (sol.type!=_VECT) 1333 return vecteur(1,gensizeerr(contextptr)); 1334 } 1335 // for each element of sol, get the first component t, subst in A 1336 int nsol=int(sol._VECTptr->size()); 1337 for (int i=0;i<nsol;i++){ 1338 gen As=ratnormal(subst(A,t,sol[i],false,contextptr),contextptr); 1339 gen smin=gnuplot_tmin,smax=gnuplot_tmax; 1340 if (av.size()>=4 && av[2].type==_VECT && av[3].type==_VECT){ 1341 smin=av[2][swapped?1:0]; 1342 smax=av[3][swapped?1:0]; 1343 } 1344 // now make the parametric curves [A,s] 1345 res.push_back(paramplotparam(gen(makevecteur(As,symb_equal(s,symb_interval(smin,smax))),_SEQ__VECT),false,contextptr)); 1346 } 1347 } // end b by equation 1348 } // end a parametric 1349 // both by equation 1350 return res; 1351 } 1352 1353 // a=hypersurface, b=curve interhypersurfacecurve(const gen & a,const gen & b,GIAC_CONTEXT)1354 vecteur interhypersurfacecurve(const gen & a,const gen &b,GIAC_CONTEXT){ 1355 gen & af=a._SYMBptr->feuille; 1356 gen & bf=b._SYMBptr->feuille; 1357 if (af.type!=_VECT || bf.type!=_VECT || af._VECTptr->empty() || bf._VECTptr->empty() ) 1358 return vecteur(1,gensizeerr(contextptr)); 1359 vecteur av=*af._VECTptr; 1360 // av[0]=point on hypersurface, av[1]=parameters 1361 gen & bf0=bf._VECTptr->front(); 1362 if (bf0.type!=_VECT || bf0._VECTptr->size()<2 || bf0._VECTptr->front().type!=_VECT) 1363 return vecteur(1,gensizeerr(contextptr)); 1364 vecteur & bv=*bf0._VECTptr; // bv[0]=point on curve, bv[1]=parameter 1365 if (av.size()==3 && av[1].type!=_VECT && av[2].type==_VECT){ 1366 // Hypersurface with an equation 1367 // av[1]=equation, av[2]=variables 1368 vecteur & vars=*av[2]._VECTptr; 1369 gen eq(subst(av[1],vars,*bv[0]._VECTptr,false,contextptr)); 1370 vecteur sol; 1371 #ifndef NO_STDEXCEPT 1372 try { 1373 #endif 1374 sol=solve(eq,bv[1],0,contextptr); 1375 #ifndef NO_STDEXCEPT 1376 } 1377 catch (std::runtime_error & ){ 1378 last_evaled_argptr(contextptr)=NULL; 1379 } 1380 #endif 1381 vecteur res; 1382 iterateur it=sol.begin(),itend=sol.end(); 1383 for (;it!=itend;++it){ 1384 res.push_back(_point(subst(bv[0],bv[1],*it,false,contextptr),contextptr)); 1385 } 1386 return res; 1387 } 1388 // Hypersurface without equation (parametrized only) 1389 if (av.size()<2 || av[1].type!=_VECT) 1390 return vecteur(1,gensizeerr(contextptr)); 1391 gen eq(av[0]-bv[0]); 1392 vecteur vars(*av[1]._VECTptr); 1393 vars.push_back(bv[1]); 1394 vecteur sol; 1395 #ifndef NO_STDEXCEPT 1396 try { 1397 #endif 1398 sol=solve(eq,vars,0,contextptr); 1399 #ifndef NO_STDEXCEPT 1400 } 1401 catch (std::runtime_error & ){ 1402 return vecteur(1,gensizeerr(contextptr)); 1403 } 1404 #endif 1405 vecteur res; 1406 iterateur it=sol.begin(),itend=sol.end(); 1407 for (;it!=itend;++it){ 1408 res.push_back(_point(subst(bv[0],bv[1],it->_VECTptr->back(),false,contextptr),contextptr)); 1409 } 1410 return res; 1411 } 1412 hyperplan2hypersurface(const gen & g)1413 gen hyperplan2hypersurface(const gen & g){ 1414 if (!g.is_symb_of_sommet(at_hyperplan)) 1415 return gensizeerr(gettext("hyperplan2hypersurface")); 1416 vecteur n,P; 1417 if (!hyperplan_normal_point(g,n,P)) 1418 return gensizeerr(gettext("hyperplan2hypersurface")); 1419 if (n.size()!=3) 1420 return gendimerr(gettext("hyperplan2hypersurface")); 1421 vecteur xyz(makevecteur(x__IDNT_e,y__IDNT_e,z__IDNT_e)); 1422 gen eq=dotvecteur(subvecteur(xyz,P),n); 1423 vecteur v1,v2; 1424 if (!normal3d(n,v1,v2)) 1425 return gensizeerr(gettext("hyperplan2hypersurface")); 1426 vecteur parameq(makevecteur(addvecteur(P,addvecteur(multvecteur(u__IDNT,v1),multvecteur(v__IDNT,v2))),makevecteur(u__IDNT,v__IDNT))); 1427 return hypersurface(parameq,eq,xyz); 1428 } 1429 hypersphere2hypersurface(const gen & g)1430 gen hypersphere2hypersurface(const gen & g){ 1431 if (!g.is_symb_of_sommet(at_hypersphere)) 1432 return gensizeerr(gettext("hypersphere2hypersurface")); 1433 vecteur xyz(makevecteur(x__IDNT_e,y__IDNT_e,z__IDNT_e)); 1434 vecteur uv(makevecteur(u__IDNT_e,v__IDNT_e)); 1435 return hypersurface(hypersphere_parameq(g,uv),hypersphere_equation(g,xyz),xyz); 1436 } 1437 1438 // Currently works only if v is made of lines 1439 // For each line, get the intersections of the polygone ABCD with the line 1440 // If there are 2 intersections return a segment, else return a point or void remove_face(const vecteur & face,const vecteur & v,GIAC_CONTEXT)1441 vecteur remove_face(const vecteur & face,const vecteur & v,GIAC_CONTEXT){ 1442 vecteur ABCD=face; 1443 if (ABCD.size()<3) 1444 return vecteur(1,gendimerr(contextptr)); 1445 if (ABCD.back()!=ABCD.front()) 1446 ABCD.push_back(ABCD.front()); 1447 vecteur res; 1448 const_iterateur it=v.begin(),itend=v.end(); 1449 for (;it!=itend;++it){ 1450 gen tmp=remove_at_pnt(*it); 1451 if (tmp.type!=_VECT || tmp._VECTptr->size()!=2) 1452 res.push_back(*it); // FIXME, arc of circles etc. 1453 vecteur v(interpolygone(ABCD,*it,contextptr)); 1454 if (is_undef(v)) 1455 return v; 1456 int s=int(v.size()); 1457 if (!s) 1458 continue; 1459 if (s==1) 1460 res.push_back(v.front()); 1461 if (s==2) 1462 res.push_back(symb_pnt(gen(makevecteur(remove_at_pnt(v.front()),remove_at_pnt(v.back())),_GROUP__VECT),contextptr)); 1463 } 1464 return res; 1465 } 1466 segments2polygone(const vecteur & v,GIAC_CONTEXT)1467 static vecteur segments2polygone(const vecteur & v,GIAC_CONTEXT){ 1468 vecteur polylines,other; 1469 int s=int(v.size()); 1470 for (int i=0;i<s;++i){ 1471 gen g=remove_at_pnt(v[i]); 1472 if (g.type!=_VECT || g._VECTptr->size()!=2){ 1473 other.push_back(g); 1474 continue; 1475 } 1476 gen a=g._VECTptr->front(),b=g._VECTptr->back(); 1477 int ps=int(polylines.size()),j=0; 1478 for (int j=0;j<ps;++j){ 1479 if (polylines[j].type==_VECT && !polylines[j]._VECTptr->empty()){ 1480 vecteur w=*polylines[j]._VECTptr; 1481 if (a==w.back()){ 1482 w.push_back(b); 1483 polylines[j]=gen(w,_GROUP__VECT); 1484 break; 1485 } 1486 if (b==w.back()){ 1487 w.push_back(a); 1488 polylines[j]=gen(w,_GROUP__VECT); 1489 break; 1490 } 1491 if (a==w.front()){ 1492 w.insert(w.begin(),b); 1493 polylines[j]=gen(w,_GROUP__VECT); 1494 break; 1495 } 1496 if (b==w.front()){ 1497 w.insert(w.begin(),a); 1498 polylines[j]=gen(w,_GROUP__VECT); 1499 break; 1500 } 1501 } 1502 } 1503 if (j==ps) 1504 polylines.push_back(g); 1505 } 1506 s=int(polylines.size()); 1507 for (int i=0;i<s;++i){ 1508 other.push_back(symb_pnt(polylines[i],contextptr)); 1509 } 1510 return other; 1511 } 1512 1513 // Intersection polyedre with something 1514 // Find intersection of each face with something 1515 // For each face, find intersection of hyperplan with something interpolyedre(const vecteur & p,const gen & bb,GIAC_CONTEXT)1516 vecteur interpolyedre(const vecteur & p,const gen & bb,GIAC_CONTEXT){ 1517 vecteur res; 1518 const_iterateur it=p.begin(),itend=p.end(); 1519 for (;it!=itend;++it){ 1520 if (it->type!=_VECT) 1521 continue; 1522 vecteur v=*it->_VECTptr; 1523 int s=int(v.size()); 1524 if (s<3) 1525 continue; 1526 gen AB(v[1]-v[0]),AC(v[2]-v[0]); 1527 if (AB.type!=_VECT || AB._VECTptr->size()!=3 || AC.type!=_VECT || AC._VECTptr->size()!=3) 1528 continue; 1529 vecteur n=cross(*AB._VECTptr,*AC._VECTptr,contextptr); 1530 gen tmp=symbolic(at_hyperplan,makesequence(n,v[0])); 1531 bool b=show_point(contextptr); 1532 show_point(false,contextptr); 1533 vecteur w(inter(tmp,bb,contextptr)); 1534 show_point(b,contextptr); 1535 vecteur restmp=remove_face(*it->_VECTptr,w,contextptr); 1536 if (is_undef(restmp)) 1537 return restmp; 1538 res=mergevecteur(res,restmp); 1539 } 1540 return segments2polygone(res,contextptr); 1541 } 1542 interhyperplan(const gen & p1,const gen & p2,GIAC_CONTEXT)1543 vecteur interhyperplan(const gen & p1,const gen & p2,GIAC_CONTEXT){ 1544 vecteur P1,n1,P2,n2; 1545 if (!hyperplan_normal_point(p1,n1,P1) || !hyperplan_normal_point(p2,n2,P2)) 1546 return vecteur(1,gensizeerr(contextptr)); 1547 vecteur n=cross(n1,n2,contextptr); // direction of intersection 1548 vecteur n3=cross(n,n1,contextptr); // perpendicular to n1 hence in P1 and not // inter 1549 // Find a point on intersection: P1+t n3 -P2 perpendicular to n2 1550 // hence (P2-P1).n2=t n3.n2 1551 gen P=do_point3d(P1-scalar_product(P1-P2,n2,contextptr)/dotvecteur(n3,n2)*n3); 1552 gen Q=do_point3d(P+n); 1553 return makevecteur(symb_pnt(gen(makevecteur(P,Q),_LINE__VECT),contextptr)); 1554 } 1555 1556 1557 // equation f -> geometric object g equation2geo3d(const gen & f0,const gen & M,const gen & x,const gen & y,const gen & z,gen & g,double umin,double umax,double ustep,double vmin,double vmax,double vstep,bool numeric,const context * contextptr)1558 static bool equation2geo3d(const gen & f0,const gen & M,const gen & x,const gen & y,const gen & z,gen & g,double umin,double umax,double ustep,double vmin,double vmax,double vstep,bool numeric,const context * contextptr){ 1559 gen f=_fxnd(remove_equal(f0),contextptr)._VECTptr->front(); 1560 gen fx(derive(f,x,contextptr)),fy(derive(f,y,contextptr)),fz(derive(f,z,contextptr)); 1561 bool fx0=is_zero(fx),fy0=is_zero(fy),fz0=is_zero(fz); 1562 if (fx0 && fy0 && fz0) 1563 return false; 1564 gen fxx(derive(fx,x,contextptr)),fxy(derive(fx,y,contextptr)),fyy(derive(fy,y,contextptr)),fxz(derive(fx,z,contextptr)),fyz(derive(fy,z,contextptr)),fzz(derive(fz,z,contextptr)); 1565 if (is_undef(fx)||is_undef(fy) || is_undef(fz) || is_undef(fxx) || is_undef(fxy) || is_undef(fxz) || is_undef(fyy) || is_undef(fyz) || is_undef(fzz)) 1566 return false; 1567 if ( is_zero(derive(fxx,x,contextptr)) && is_zero(derive(fxy,x,contextptr)) && is_zero(derive(fyy,x,contextptr)) && is_zero(derive(fxz,x,contextptr)) && is_zero(derive(fyz,x,contextptr)) && is_zero(derive(fzz,x,contextptr)) && 1568 is_zero(derive(fxx,y,contextptr)) && is_zero(derive(fxy,y,contextptr)) && is_zero(derive(fyy,y,contextptr)) && is_zero(derive(fxz,y,contextptr)) && is_zero(derive(fyz,y,contextptr)) && is_zero(derive(fzz,y,contextptr)) && 1569 is_zero(derive(fxx,z,contextptr)) && is_zero(derive(fxy,z,contextptr)) && is_zero(derive(fyy,z,contextptr)) && is_zero(derive(fxz,z,contextptr)) && is_zero(derive(fyz,z,contextptr)) && is_zero(derive(fzz,z,contextptr)) 1570 ){ 1571 vecteur vxyz(makevecteur(x,y,z)),v0(3,0); 1572 gen c=ratnormal(subst(f,vxyz,v0,false,contextptr),contextptr); 1573 fxx=ratnormal(fxx,contextptr); fyy=ratnormal(fyy,contextptr); fxy=ratnormal(fxy,contextptr); 1574 fxz=ratnormal(fxz,contextptr); fyz=ratnormal(fyz,contextptr); fzz=ratnormal(fzz,contextptr); 1575 if (is_zero(fxy) && is_zero(fxz) && is_zero(fyz)){ 1576 if (is_zero(fxx) && is_zero(fyy) && is_zero(fzz)){ 1577 gen d=gcd(gcd(fx,fy),fz); 1578 fx=normal(fx/d,contextptr); fy=normal(fy/d,contextptr); fz=normal(fz/d,contextptr); c=normal(c/d,contextptr); 1579 vecteur n(makevecteur(fx,fy,fz)); 1580 // plan 1581 if (!fx0){ 1582 gen tmp=makevecteur(ratnormal(-c/fx,contextptr),0,0); 1583 g=symbolic(at_hyperplan,makesequence(n,tmp)); 1584 } 1585 else { 1586 if (!fy0){ 1587 gen tmp=makevecteur(0,ratnormal(-c/fy,contextptr),0); 1588 g=symbolic(at_hyperplan,makesequence(n,tmp)); 1589 } 1590 else { 1591 gen tmp=makevecteur(0,0,ratnormal(-c/fz,contextptr)); 1592 g=symbolic(at_hyperplan,makesequence(n,tmp)); 1593 } 1594 } 1595 return true; 1596 } 1597 // may check for a sphere (fxx=fyy=fzz) 1598 } 1599 // conique 1600 gen x0,y0,z0,equation_reduite; 1601 vecteur V0,V1,V2,param_curves,propre,centre; 1602 quadrique_reduite(f,M,vxyz,x0,y0,z0,V0,V1,V2,propre,equation_reduite,param_curves,centre,numeric,contextptr); 1603 vecteur res; 1604 int n=int(param_curves.size()); 1605 for (int i=0;i<n;++i){ 1606 gen & obj=param_curves[i]; 1607 if (obj.type==_VECT){ 1608 vecteur & objv=*obj._VECTptr; 1609 int s=int(objv.size()); 1610 if (s==2) 1611 res.push_back(obj); 1612 if (s==5){ 1613 gen tmp=paramplotparam(gen(objv,_SEQ__VECT),false,contextptr); 1614 tmp=remove_at_pnt(tmp); 1615 if (tmp.is_symb_of_sommet(at_hypersurface) && tmp._SYMBptr->feuille.type==_VECT){ 1616 vecteur tmpv=*tmp._SYMBptr->feuille._VECTptr; 1617 if (tmpv.size()==3){ 1618 tmpv[1]=f; 1619 tmpv[2]=vxyz; 1620 tmp=symbolic(at_hypersurface,gen(tmpv,_SEQ__VECT)); 1621 } 1622 } 1623 res.push_back(tmp); 1624 } 1625 } 1626 else 1627 res.push_back(obj); 1628 } 1629 g= (res.size()==1)? res.front() : res; // gen(res,_SEQ__VECT); 1630 return true; 1631 } 1632 return false; 1633 } 1634 1635 #if !defined(RTOS_THREADX) && !defined(EMCC) 1636 // 3-d implicit surface using the marching cube algorithm 1637 /* Adapted from http://astronomy.swin.edu.au/~pbourke/modelling/ 1638 by Paul Bourke 1639 Given a grid cell and an isolevel, calculate the triangular 1640 facets required to represent the isosurface through the cell. 1641 Return the number of triangular facets, the array "triangles" 1642 will be loaded up with the vertices at most 5 triangular facets. 1643 0 will be returned if the grid cell is either totally above 1644 of totally below the isolevel. 1645 */ 1646 int const edgeTable[256]={ 1647 0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 1648 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, 1649 0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 1650 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, 1651 0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c, 1652 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 1653 0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac, 1654 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, 1655 0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c, 1656 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60, 1657 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc, 1658 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, 1659 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c, 1660 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, 1661 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc , 1662 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 1663 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 1664 0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0, 1665 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 1666 0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650, 1667 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 1668 0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 1669 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, 1670 0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460, 1671 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 1672 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0, 1673 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 1674 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230, 1675 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 1676 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190, 1677 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 1678 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 }; 1679 int const triTable[256][16] = 1680 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1681 {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1682 {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1683 {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1684 {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1685 {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1686 {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1687 {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1}, 1688 {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1689 {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1690 {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1691 {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1}, 1692 {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1693 {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1}, 1694 {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1}, 1695 {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1696 {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1697 {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1698 {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1699 {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1}, 1700 {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1701 {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1}, 1702 {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, 1703 {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1}, 1704 {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1705 {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1}, 1706 {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, 1707 {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1}, 1708 {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1}, 1709 {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1}, 1710 {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1}, 1711 {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1}, 1712 {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1713 {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1714 {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1715 {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1}, 1716 {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1717 {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, 1718 {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1}, 1719 {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1}, 1720 {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1721 {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, 1722 {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, 1723 {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1}, 1724 {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1}, 1725 {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1}, 1726 {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1}, 1727 {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1}, 1728 {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1729 {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1}, 1730 {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1}, 1731 {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1732 {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1}, 1733 {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1}, 1734 {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1}, 1735 {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1}, 1736 {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1}, 1737 {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1}, 1738 {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1}, 1739 {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1}, 1740 {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1}, 1741 {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1}, 1742 {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1}, 1743 {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1744 {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1745 {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1746 {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1747 {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, 1748 {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1749 {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1}, 1750 {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1}, 1751 {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1}, 1752 {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1753 {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, 1754 {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, 1755 {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1}, 1756 {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1}, 1757 {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1}, 1758 {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1}, 1759 {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1}, 1760 {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1761 {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1}, 1762 {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, 1763 {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1}, 1764 {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1}, 1765 {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1}, 1766 {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1}, 1767 {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1}, 1768 {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, 1769 {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1}, 1770 {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1}, 1771 {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1}, 1772 {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1}, 1773 {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1}, 1774 {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1}, 1775 {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1}, 1776 {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1777 {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1}, 1778 {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1}, 1779 {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1}, 1780 {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1}, 1781 {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1}, 1782 {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1783 {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1}, 1784 {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1}, 1785 {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1}, 1786 {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1}, 1787 {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1}, 1788 {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1}, 1789 {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1}, 1790 {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1}, 1791 {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1792 {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1}, 1793 {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1}, 1794 {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1}, 1795 {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1}, 1796 {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1}, 1797 {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1}, 1798 {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1}, 1799 {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1800 {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1}, 1801 {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1}, 1802 {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1}, 1803 {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1}, 1804 {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1}, 1805 {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1806 {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1}, 1807 {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1808 {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1809 {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1810 {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1811 {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, 1812 {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1813 {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, 1814 {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, 1815 {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1}, 1816 {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1817 {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1}, 1818 {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1}, 1819 {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1}, 1820 {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1}, 1821 {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1}, 1822 {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1}, 1823 {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1}, 1824 {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1825 {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1}, 1826 {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1}, 1827 {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1}, 1828 {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1}, 1829 {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1}, 1830 {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1}, 1831 {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1}, 1832 {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1}, 1833 {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1834 {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1}, 1835 {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1}, 1836 {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1}, 1837 {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1}, 1838 {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1}, 1839 {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1840 {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1841 {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, 1842 {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, 1843 {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1}, 1844 {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, 1845 {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1}, 1846 {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1}, 1847 {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1}, 1848 {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1}, 1849 {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1}, 1850 {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1}, 1851 {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1}, 1852 {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1}, 1853 {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1}, 1854 {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1}, 1855 {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1}, 1856 {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1}, 1857 {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1}, 1858 {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1}, 1859 {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1}, 1860 {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1}, 1861 {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1}, 1862 {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1}, 1863 {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1}, 1864 {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1}, 1865 {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1}, 1866 {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1}, 1867 {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1868 {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1}, 1869 {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1}, 1870 {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1871 {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1872 {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1873 {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1}, 1874 {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1}, 1875 {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1}, 1876 {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1}, 1877 {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1}, 1878 {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1}, 1879 {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1}, 1880 {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1}, 1881 {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1}, 1882 {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1}, 1883 {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1}, 1884 {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1885 {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1}, 1886 {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1}, 1887 {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1888 {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1}, 1889 {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1}, 1890 {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1}, 1891 {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1}, 1892 {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1}, 1893 {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1}, 1894 {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1}, 1895 {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1896 {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1}, 1897 {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1}, 1898 {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1}, 1899 {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1}, 1900 {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1}, 1901 {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1902 {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1}, 1903 {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1904 {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1}, 1905 {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1}, 1906 {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1}, 1907 {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1}, 1908 {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1}, 1909 {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1}, 1910 {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1}, 1911 {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1}, 1912 {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1}, 1913 {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1}, 1914 {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1}, 1915 {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1916 {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1}, 1917 {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1}, 1918 {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1919 {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1920 {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1921 {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1}, 1922 {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1}, 1923 {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1924 {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1}, 1925 {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1}, 1926 {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1927 {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1928 {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1}, 1929 {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1930 {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1}, 1931 {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1932 {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1933 {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1934 {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 1935 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}; 1936 1937 struct GRIDCELL { 1938 vecteur p[8]; 1939 double val[8]; 1940 } ; 1941 VertexInterp(double isolevel,const vecteur & p1,const vecteur & p2,double valp1,double valp2)1942 static vecteur VertexInterp(double isolevel,const vecteur & p1,const vecteur & p2,double valp1,double valp2){ 1943 double mu; 1944 1945 if (absdouble(isolevel-valp1) < 0.00001) 1946 return(p1); 1947 if (absdouble(isolevel-valp2) < 0.00001) 1948 return(p2); 1949 if (absdouble(valp1-valp2) < 0.00001) 1950 return(p1); 1951 mu = (isolevel - valp1) / (valp2 - valp1); 1952 1953 return addvecteur(p1,multvecteur(mu,subvecteur(p2,p1))); 1954 } 1955 in_plotimplicit(const gen & f_orig,const gen & x,const gen & y,const gen & z,double xmin,double xmax,double ymin,double ymax,double zmin,double zmax,int nxstep,int nystep,int nzstep,double eps,const vecteur & attributs,const context * contextptr)1956 static gen in_plotimplicit(const gen& f_orig,const gen&x,const gen & y,const gen & z,double xmin,double xmax,double ymin,double ymax,double zmin,double zmax,int nxstep,int nystep,int nzstep,double eps,const vecteur & attributs,const context * contextptr){ 1957 if (f_orig.is_symb_of_sommet(at_inv) || (is_zero(derive(f_orig,x,contextptr)) && is_zero(derive(f_orig,y,contextptr))) ) 1958 return vecteur(0); // gen(vecteur(0),_SEQ__VECT); 1959 if (f_orig.is_symb_of_sommet(at_prod) && f_orig._SYMBptr->feuille.type==_VECT){ 1960 vecteur res; 1961 vecteur & fv = *f_orig._SYMBptr->feuille._VECTptr; 1962 int s=int(fv.size()); 1963 for (int i=0;i<s;++i){ 1964 gen tmp=in_plotimplicit(fv[i],x,y,z,xmin,xmax,ymin,ymax,zmin,zmax, 1965 nxstep,nystep,nzstep,eps,attributs,contextptr); 1966 res=mergevecteur(res,gen2vecteur(tmp)); 1967 } 1968 return res; // gen(res,_SEQ__VECT); 1969 } 1970 if (f_orig.is_symb_of_sommet(at_pow)){ 1971 gen farg=f_orig._SYMBptr->feuille; 1972 if (farg.type==_VECT && farg._VECTptr->size()==2){ 1973 gen arg=farg._VECTptr->front(); 1974 gen expo=farg._VECTptr->back(); 1975 if (ck_is_positive(expo,contextptr)) 1976 return in_plotimplicit(farg,x,y,z,xmin,xmax,ymin,ymax,zmin,zmax,nxstep,nystep,nzstep,eps,attributs,contextptr); 1977 else 1978 return vecteur(0); // gen(vecteur(0),_SEQ__VECT); 1979 } 1980 } 1981 gen attribut=attributs.empty()?default_color(contextptr):attributs[0]; 1982 gen lieu_geo; 1983 if (equation2geo3d(f_orig,undef,x,y,z,lieu_geo,gnuplot_tmin,gnuplot_tmax,gnuplot_tstep,gnuplot_tmin,gnuplot_tmax,gnuplot_tstep,true,contextptr)) 1984 return put_attributs(lieu_geo,attributs,contextptr); 1985 //return undef; 1986 if (nxstep*double(nystep)*nzstep>8000){ 1987 nxstep=10; 1988 nystep=10; 1989 nzstep=10; 1990 } 1991 double xstep=(xmax-xmin)/nxstep; 1992 double ystep=(ymax-ymin)/nystep; 1993 double zstep=(zmax-zmin)/nzstep; 1994 identificateur xloc(" xloc"),yloc(" yloc"),zloc(" zloc"); 1995 vecteur xyz(makevecteur(x,y,z)),locvar(makevecteur(xloc,yloc,zloc)); 1996 gen f=quotesubst(f_orig,xyz,locvar,contextptr).evalf(1,contextptr); 1997 vecteur localvar(makevecteur(xloc,yloc,zloc)); 1998 context * newcontextptr=(context *) contextptr; 1999 int protect=giac_bind(makevecteur(xmin,ymin,zmin),localvar,newcontextptr); 2000 vector< vector< vector<double> > > fxyz(nxstep+1,vector< vector<double> >(nystep+1,vector<double>(nzstep+1))); 2001 gen gtmp; 2002 // initialize each cell value of f 2003 local_sto_double(xmin,xloc,newcontextptr); 2004 // xloc.localvalue->back()._DOUBLE_val = xmin; 2005 for (int i=0;i<=nxstep;++i){ 2006 local_sto_double(ymin,yloc,newcontextptr); 2007 // yloc.localvalue->back()._DOUBLE_val = ymin ; 2008 for (int j=0;j<=nystep;++j){ 2009 local_sto_double(zmin,zloc,newcontextptr); 2010 // zloc.localvalue->back()._DOUBLE_val = zmin ; 2011 for (int k=0;k<=nzstep;++k){ 2012 gtmp=f.evalf_double(eval_level(contextptr),newcontextptr); 2013 if (gtmp.type==_DOUBLE_) 2014 fxyz[i][j][k]=gtmp._DOUBLE_val==0?1e-200:gtmp._DOUBLE_val; 2015 else 2016 fxyz[i][j][k]=0; 2017 local_sto_double_increment(zstep,zloc,newcontextptr); 2018 // zloc.localvalue->back()._DOUBLE_val += zstep; 2019 } 2020 local_sto_double_increment(ystep,yloc,newcontextptr); 2021 // yloc.localvalue->back()._DOUBLE_val += ystep; 2022 } 2023 local_sto_double_increment(xstep,xloc,newcontextptr); 2024 // xloc.localvalue->back()._DOUBLE_val += xstep; 2025 } 2026 leave(protect,localvar,newcontextptr); 2027 2028 GRIDCELL grid; 2029 vecteur triangle(4),triangulation; 2030 double xcur=xmin,ycur,zcur; 2031 for (int i=0;i<nxstep;++i,xcur+=xstep){ 2032 ycur=ymin; 2033 for (int j=0;j<nystep;++j,ycur+=ystep){ 2034 zcur=zmin; 2035 for (int k=0;k<nzstep;++k,zcur+=zstep){ 2036 // Set current gridcell 2037 grid.val[0]=fxyz[i][j][k]; 2038 grid.p[0]=makevecteur(xcur,ycur,zcur); 2039 grid.val[1]=fxyz[i+1][j][k]; 2040 grid.p[1]=makevecteur(xcur+xstep,ycur,zcur); 2041 grid.val[2]=fxyz[i+1][j][k+1]; 2042 grid.p[2]=makevecteur(xcur+xstep,ycur,zcur+zstep); 2043 grid.val[3]=fxyz[i][j][k+1]; 2044 grid.p[3]=makevecteur(xcur,ycur,zcur+zstep); 2045 grid.val[4]=fxyz[i][j+1][k]; 2046 grid.p[4]=makevecteur(xcur,ycur+ystep,zcur); 2047 grid.val[5]=fxyz[i+1][j+1][k]; 2048 grid.p[5]=makevecteur(xcur+xstep,ycur+ystep,zcur); 2049 grid.val[6]=fxyz[i+1][j+1][k+1]; 2050 grid.p[6]=makevecteur(xcur+xstep,ycur+ystep,zcur+zstep); 2051 grid.val[7]=fxyz[i][j+1][k+1]; 2052 grid.p[7]=makevecteur(xcur,ycur+ystep,zcur+zstep); 2053 2054 /* 2055 Determine the index into the edge table which 2056 tells us which vertices are inside of the surface 2057 */ 2058 int cubeindex=0; 2059 if (grid.val[0] < 0) cubeindex |= 1; 2060 if (grid.val[1] < 0) cubeindex |= 2; 2061 if (grid.val[2] < 0) cubeindex |= 4; 2062 if (grid.val[3] < 0) cubeindex |= 8; 2063 if (grid.val[4] < 0) cubeindex |= 16; 2064 if (grid.val[5] < 0) cubeindex |= 32; 2065 if (grid.val[6] < 0) cubeindex |= 64; 2066 if (grid.val[7] < 0) cubeindex |= 128; 2067 2068 /* Cube is entirely in/out of the surface */ 2069 if (edgeTable[cubeindex] == 0) 2070 continue; 2071 2072 vecteur vertlist[12]; 2073 /* Find the vertices where the surface intersects the cube */ 2074 if (edgeTable[cubeindex] & 1) 2075 vertlist[0] = 2076 VertexInterp(0,grid.p[0],grid.p[1],grid.val[0],grid.val[1]); 2077 if (edgeTable[cubeindex] & 2) 2078 vertlist[1] = 2079 VertexInterp(0,grid.p[1],grid.p[2],grid.val[1],grid.val[2]); 2080 if (edgeTable[cubeindex] & 4) 2081 vertlist[2] = 2082 VertexInterp(0,grid.p[2],grid.p[3],grid.val[2],grid.val[3]); 2083 if (edgeTable[cubeindex] & 8) 2084 vertlist[3] = 2085 VertexInterp(0,grid.p[3],grid.p[0],grid.val[3],grid.val[0]); 2086 if (edgeTable[cubeindex] & 16) 2087 vertlist[4] = 2088 VertexInterp(0,grid.p[4],grid.p[5],grid.val[4],grid.val[5]); 2089 if (edgeTable[cubeindex] & 32) 2090 vertlist[5] = 2091 VertexInterp(0,grid.p[5],grid.p[6],grid.val[5],grid.val[6]); 2092 if (edgeTable[cubeindex] & 64) 2093 vertlist[6] = 2094 VertexInterp(0,grid.p[6],grid.p[7],grid.val[6],grid.val[7]); 2095 if (edgeTable[cubeindex] & 128) 2096 vertlist[7] = 2097 VertexInterp(0,grid.p[7],grid.p[4],grid.val[7],grid.val[4]); 2098 if (edgeTable[cubeindex] & 256) 2099 vertlist[8] = 2100 VertexInterp(0,grid.p[0],grid.p[4],grid.val[0],grid.val[4]); 2101 if (edgeTable[cubeindex] & 512) 2102 vertlist[9] = 2103 VertexInterp(0,grid.p[1],grid.p[5],grid.val[1],grid.val[5]); 2104 if (edgeTable[cubeindex] & 1024) 2105 vertlist[10] = 2106 VertexInterp(0,grid.p[2],grid.p[6],grid.val[2],grid.val[6]); 2107 if (edgeTable[cubeindex] & 2048) 2108 vertlist[11] = 2109 VertexInterp(0,grid.p[3],grid.p[7],grid.val[3],grid.val[7]); 2110 2111 /* Create the triangles */ 2112 2113 for (int i=0;triTable[cubeindex][i]!=-1;i+=3) { 2114 triangle[0]=gen(vertlist[triTable[cubeindex][i]],_POINT__VECT); 2115 triangle[1]=gen(vertlist[triTable[cubeindex][i+1]],_POINT__VECT); 2116 triangle[2]=gen(vertlist[triTable[cubeindex][i+2]],_POINT__VECT); 2117 triangle[3]=triangle[0]; 2118 triangulation.push_back(gen(triangle,_GROUP__VECT)); 2119 } 2120 2121 } // end for k 2122 } // end for j 2123 } // end for i 2124 2125 // create hypersurface 2126 gen tmp=gen(makevecteur(undef,makevecteur(x,y,z),makevecteur(xmin,ymin,zmin),makevecteur(xmax,ymax,zmax),gen(triangulation,_POLYEDRE__VECT)),_PNT__VECT); 2127 gen r=symb_pnt(hypersurface(tmp,f,makevecteur(x,y,z)),attribut,contextptr); 2128 return r; 2129 } 2130 #else //RTOS_THREADX in_plotimplicit(const gen & f_orig,const gen & x,const gen & y,const gen & z,double xmin,double xmax,double ymin,double ymax,double zmin,double zmax,int nxstep,int nystep,int nzstep,double eps,const vecteur & attributs,const context * contextptr)2131 static gen in_plotimplicit(const gen& f_orig,const gen&x,const gen & y,const gen & z,double xmin,double xmax,double ymin,double ymax,double zmin,double zmax,int nxstep,int nystep,int nzstep,double eps,const vecteur & attributs,const context * contextptr){ 2132 return undef; 2133 } 2134 #endif 2135 plotimplicit(const gen & f_orig,const gen & x,const gen & y,const gen & z,double xmin,double xmax,double ymin,double ymax,double zmin,double zmax,int nxstep,int nystep,int nzstep,double eps,const vecteur & attributs,bool unfactored,const context * contextptr)2136 gen plotimplicit(const gen& f_orig,const gen&x,const gen & y,const gen & z,double xmin,double xmax,double ymin,double ymax,double zmin,double zmax,int nxstep,int nystep,int nzstep,double eps,const vecteur & attributs,bool unfactored,const context * contextptr){ 2137 if ( x.type!=_IDNT || y.type!=_IDNT || z.type!=_IDNT ) 2138 return gensizeerr(gettext("Variables must be free")); 2139 if (!nystep || !nzstep){ 2140 nxstep=int(std::sqrt(double(absdouble(nxstep)))); 2141 nystep=nxstep; 2142 nzstep=nxstep; 2143 } 2144 gen ff( (unfactored || has_num_coeff(f_orig))?f_orig:factor(f_orig,false,contextptr)); 2145 return in_plotimplicit(ff,x,y,z,xmin,xmax,ymin,ymax,zmin,zmax,nxstep,nystep,nzstep,eps,attributs,contextptr); 2146 } 2147 2148 // FIXME Move to plot.cc is3d(const gen & g)2149 bool is3d(const gen & g){ 2150 if (g.type==_VECT) 2151 return !g._VECTptr->empty() && is3d(g._VECTptr->back()); 2152 if (g.is_symb_of_sommet(at_animation)) 2153 return is3d(g._SYMBptr->feuille); 2154 if (!g.is_symb_of_sommet(at_pnt)) 2155 return false; 2156 gen f =g._SYMBptr->feuille; 2157 if (f.type!=_VECT || f._VECTptr->empty()) 2158 return false; 2159 f=f._VECTptr->front(); 2160 if (f.type==_VECT){ 2161 if (f.subtype==_POLYEDRE__VECT || f.subtype==_POINT__VECT) 2162 return true; 2163 if (f._VECTptr->size()==3 && f.subtype!=_GROUP__VECT && f.subtype!=_LINE__VECT && f.subtype!=_HALFLINE__VECT){ 2164 vecteur & v =*f._VECTptr; 2165 return v[0].type!=_CPLX && v[1].type!=_CPLX && v[2].type!=_CPLX; 2166 } 2167 if (!f._VECTptr->empty()) 2168 return check3dpoint(f._VECTptr->back()); 2169 } 2170 if (f.type!=_SYMB) 2171 return false; 2172 if (f._SYMBptr->sommet==at_hyperplan || f._SYMBptr->sommet==at_hypersphere || f._SYMBptr->sommet == at_hypersurface) 2173 return true; 2174 if (f._SYMBptr->sommet==at_curve && f._SYMBptr->feuille.type==_VECT && !f._SYMBptr->feuille._VECTptr->empty()){ 2175 // is it a 3-d curve? 2176 f = f._SYMBptr->feuille._VECTptr->front(); 2177 // f = vect[ pnt,var,xmin,xmax ] 2178 if (f.type==_VECT && !f._VECTptr->empty()) 2179 return check3dpoint(f._VECTptr->front()); 2180 } 2181 return false; 2182 } 2183 _quadrique(const gen & args,GIAC_CONTEXT)2184 gen _quadrique(const gen & args,GIAC_CONTEXT){ 2185 if ( args.type==_STRNG && args.subtype==-1) return args; 2186 if (args.type!=_VECT) 2187 return _quadrique(gen(vecteur(1,args),_SEQ__VECT),contextptr); 2188 vecteur attributs(1,default_color(contextptr)); 2189 vecteur & v =*args._VECTptr; 2190 int s=read_attributs(v,attributs,contextptr); 2191 if (!s) return gendimerr(contextptr); 2192 bool numeric=true; 2193 if (v[s-1]==at_exact){ 2194 numeric=false; 2195 --s; 2196 if (!s) return gendimerr(contextptr); 2197 } 2198 gen g; 2199 if (s==1){ 2200 if (equation2geo3d(v[0],undef,x__IDNT_e,y__IDNT_e,z__IDNT_e,g,gnuplot_tmin,gnuplot_tmax,gnuplot_tstep,gnuplot_tmin,gnuplot_tmax,gnuplot_tstep,numeric,contextptr)) 2201 return put_attributs(g,attributs,contextptr); 2202 } 2203 if (s<=5) 2204 return _plotimplicit(args,contextptr); 2205 if (s==9){ // defined by 9 points 2206 vecteur w(9),m(9),ligne(10); 2207 gen wx,wy,wz; 2208 // find m such that m*[a,b,c,d,e,f,g,h,i,j]=0 where 2209 // a*x^2+b*x*y+c*y^2+d*x+e*y+f+g*z^2+h*z*x+i*z*y+j*z=0 2210 for (int i=0;i<9;++i){ 2211 w[i]=remove_at_pnt(v[i]); 2212 if (w[i].type!=_VECT || w[i]._VECTptr->size()!=3) 2213 return gensizeerr(contextptr); 2214 vecteur & wv = *w[i]._VECTptr; 2215 wx=wv[0]; 2216 wy=wv[1]; 2217 wz=wv[2]; 2218 ligne[0]=wx*wx; 2219 ligne[1]=wx*wy; 2220 ligne[2]=wy*wy; 2221 ligne[3]=wx; 2222 ligne[4]=wy; 2223 ligne[5]=1; 2224 ligne[6]=wz*wz; 2225 ligne[7]=wz*wx; 2226 ligne[8]=wz*wy; 2227 ligne[9]=wz; 2228 m[i]=ligne; 2229 } 2230 // find ker(m) 2231 gen me=m; // exact(m,contextptr); 2232 vecteur base; 2233 if (me.type==_VECT) 2234 base=mker(*me._VECTptr,contextptr); 2235 if (is_undef(base) || base.empty() || base.front().type!=_VECT || base.front()._VECTptr->size()!=6) 2236 return gensizeerr(gettext("Bug in quadrique reducing ")+gen(m).print(contextptr)); 2237 vecteur & res = *base.front()._VECTptr; 2238 identificateur x(" x"),y(" y"),z( "z"); 2239 gen eq=res[0]*x*x+res[1]*x*y+res[2]*y*y+res[3]*x+res[4]*y+res[5]+res[6]*z*z+res[7]*z*x+res[8]*z*y+res[9]*z; 2240 if (equation2geo3d(eq,w[0],x,y,z,g,gnuplot_tmin,gnuplot_tmax,gnuplot_tstep,gnuplot_tmin,gnuplot_tmax,gnuplot_tstep,numeric,contextptr)) 2241 return put_attributs(g,attributs,contextptr); 2242 else 2243 return gensizeerr(gettext("Bug in quadrique, equation ")+eq.print(contextptr)); 2244 } 2245 return gendimerr(contextptr); 2246 } 2247 static const char _quadrique_s []="quadric"; 2248 static define_unary_function_eval (__quadrique,&_quadrique,_quadrique_s); 2249 define_unary_function_ptr5( at_quadrique ,alias_at_quadrique,&__quadrique,0,true); 2250 est_cospherique(const gen & a,const gen & b,const gen & c,const gen & d,const gen & f,GIAC_CONTEXT)2251 bool est_cospherique(const gen & a,const gen & b,const gen & c,const gen & d,const gen & f,GIAC_CONTEXT){ 2252 gen ab=b-a,ac=c-a,ad=d-a,af=f-a; 2253 if (is_zero(ab) || is_zero(ac) || is_zero(ad) || is_zero(af)) 2254 return true; 2255 return est_coplanaire(a+ab/abs_norm2(ab,contextptr),a+ac/abs_norm2(ac,contextptr),a+ad/abs_norm2(ad,contextptr),a+af/abs_norm2(af,contextptr),contextptr); 2256 } _est_cospherique(const gen & args,GIAC_CONTEXT)2257 gen _est_cospherique(const gen & args,GIAC_CONTEXT){ 2258 if ( args.type==_STRNG && args.subtype==-1) return args; 2259 if (args.type!=_VECT) 2260 return gensizeerr(contextptr); 2261 vecteur & v=*args._VECTptr ; 2262 int s=int(v.size()); 2263 gen a(v[0]),b(undef),c(undef),d(undef); 2264 for (int i=1;i<s;++i){ 2265 if (is_undef(b)){ 2266 if (!is_zero(v[0]-v[i])) 2267 b=v[i]; 2268 } 2269 else { 2270 if (est_aligne(a,b,v[i],contextptr)) 2271 return 0; 2272 if (is_undef(c)) 2273 c=v[i]; 2274 else { 2275 if (est_cocyclique(a,b,c,v[i],contextptr)) 2276 continue; 2277 if (is_undef(d)) 2278 d=v[i]; 2279 else { 2280 if (!est_cospherique(a,b,c,d,v[i],contextptr)) 2281 return 0; 2282 } 2283 } 2284 } 2285 } 2286 return 1; 2287 } 2288 static const char _est_cospherique_s []="is_cospherical"; 2289 static define_unary_function_eval (__est_cospherique,&_est_cospherique,_est_cospherique_s); 2290 define_unary_function_ptr5( at_est_cospherique ,alias_at_est_cospherique,&__est_cospherique,0,true); 2291 in_convert3d(const gen & h0,GIAC_CONTEXT)2292 gen in_convert3d(const gen & h0,GIAC_CONTEXT){ 2293 if (h0.type!=_VECT) 2294 return h0; 2295 // convert complex to 3d vectors 2296 vecteur w=*h0._VECTptr; 2297 gen r,I; 2298 for (unsigned i=0;i<w.size();++i){ 2299 reim(w[i],r,I,contextptr); 2300 w[i]=gen(makevecteur(r,I,0),_POINT__VECT); 2301 } 2302 return gen(w,h0.subtype); 2303 } 2304 convert3d(const gen & g,GIAC_CONTEXT)2305 gen convert3d(const gen & g,GIAC_CONTEXT){ 2306 if (g.type==_VECT) 2307 return apply(g,convert3d,contextptr); 2308 if (!g.is_symb_of_sommet(at_pnt)) 2309 return g; 2310 gen h=g._SYMBptr->feuille; 2311 if (h.type!=_VECT || h._VECTptr->size()<2) 2312 return g; 2313 vecteur v=*h._VECTptr; 2314 gen h0=v.front(); 2315 if (h0.is_symb_of_sommet(at_curve)){ 2316 gen c=h0._SYMBptr->feuille; 2317 if (c.type!=_VECT || c._VECTptr->size()<2) 2318 return g; 2319 vecteur vc=*c._VECTptr; 2320 gen c0=vc[0]; 2321 gen c1=vc[1]; 2322 if (c0.type==_VECT && !c0._VECTptr->empty()){ 2323 vecteur v0=*c0._VECTptr; 2324 gen r,i; 2325 reim(v0[0],r,i,contextptr); 2326 v0[0]=gen(makevecteur(r,i,0),_POINT__VECT); 2327 vc[0]=gen(v0,c0.subtype); 2328 } 2329 vc[1]=in_convert3d(c1,contextptr); 2330 h0=symbolic(at_curve,gen(vc,h0.subtype)); 2331 v.front()=h0; 2332 h=gen(v,h.subtype); 2333 return symbolic(at_pnt,h); 2334 } 2335 if (h0.is_symb_of_sommet(at_cercle)){ 2336 gen centre,rayon; 2337 if (!centre_rayon(h0,centre,rayon,false,contextptr)) 2338 return gensizeerr(contextptr); // don't care about radius sign 2339 double nstep=50.0; 2340 gen tmin=0,tmax=2*M_PI; 2341 if (h0.type==_VECT && h0._VECTptr->size()>4){ 2342 tmin=(*h0._VECTptr)[2]; 2343 tmax=(*h0._VECTptr)[3]; 2344 } 2345 gen tstep=(tmax-tmin)/nstep; 2346 vecteur curve; 2347 for (int i=0;i<=nstep;++i){ 2348 gen t=tmin+i*tstep; 2349 if (i==nstep){ 2350 if (h0.type==_VECT && h0._VECTptr->size()>4) 2351 t=tmax; 2352 else 2353 t=tmin; 2354 } 2355 gen z=centre+rayon*exp(cst_i*t,contextptr),zx,zy; 2356 reim(z,zx,zy,contextptr); 2357 curve.push_back(gen(makevecteur(zx,zy,0),_POINT__VECT)); 2358 } 2359 v[0]=gen(curve,_GROUP__VECT); 2360 return symbolic(at_pnt,gen(v,h.subtype)); 2361 } 2362 if (h0.type==_VECT){ 2363 v.front()=in_convert3d(h0,contextptr); 2364 h=gen(v,h.subtype); 2365 return symbolic(at_pnt,h); 2366 } 2367 gen r,i; 2368 reim(h0,r,i,contextptr); 2369 h0=gen(makevecteur(r,i,0),_POINT__VECT); 2370 v[0]=h0; 2371 h=gen(v,h.subtype); 2372 return symbolic(at_pnt,h); 2373 } 2374 static const char _convert3d_s []="convert3d"; 2375 static define_unary_function_eval (__convert3d,&convert3d,_convert3d_s); 2376 define_unary_function_ptr5( at_convert3d ,alias_at_convert3d,&__convert3d,0,true); 2377 2378 2379 #ifndef NO_NAMESPACE_GIAC 2380 } // namespace giac 2381 #endif // ndef NO_NAMESPACE_GIAC 2382