00001
00002
00003
00004
00005 # pragma once
00006
00007
00009 # include <realroot/ring_monomial_tensor.hpp>
00011
00012 # include <realroot/polynomial.hpp>
00013 # include <realroot/homography.hpp>
00014 # include <realroot/sign_variation.hpp>
00015 # include <realroot/Interval.hpp>
00016
00017 # include <stack>
00018
00019 # include <realroot/solver_ucf.hpp>
00020 # include <realroot/solver_uv_bspline.hpp>
00021
00022
00023 # define TMPL template<class POL>
00024 # define DPOL polynomial<double,with<MonomialTensor> >
00025 # define BOX box_rep<POL>
00026
00027 # define INFO 1
00028 # define TODOUBLE as<double>
00029
00030
00031 namespace mmx {
00032
00033 namespace realroot {
00034
00035 template<class POL, class FT> int topological_degree_2d(BOX* b);
00036
00038 TMPL
00039 inline bool no_variation(POL & p)
00040 {
00041 return !has_sign_variation(p.begin(),p.end()) ;
00042 };
00043
00044
00045
00046
00047 TMPL
00048 inline DPOL * jacobian(Seq<POL*>& S0)
00049 {
00050 int n=S0.size();
00051 DPOL * jac;
00052 POL df;
00053 jac= new DPOL [n*n];
00054
00055 for( int i=0;i<n;i++)
00056 {
00057 for( int k=0;k<n;k++)
00058 {
00059 df= diff( *(S0[i]),k);
00060 let::assign( jac[i*n+k] , df);
00061 }
00062
00063 }
00064 return jac;
00065 }
00066
00067
00068 template<class T> inline
00069 void precondition (Seq<T> val, DPOL * mat, int n)
00070 {
00071 int m(val.size());
00072 T *ev;
00073 ev= new T [m*n];
00074
00075 for( int i=0;i<m;i++)
00076 {
00077 for( int k=0;k<n;k++)
00078 {
00079 tensor::eval( ev[i*m+k] , mat[i*m+k].rep(),
00080 val , m );
00081 }
00082 }
00083 return ev;
00084 }
00085
00086
00087
00088
00089 template<class T> inline
00090 T * eval_poly_matrix(Seq<T> val, DPOL * mat, int n)
00091 {
00092 int m(val.size());
00093 T *ev;
00094 ev= new T [m*n];
00095
00096 for( int i=0;i<m;i++)
00097 {
00098 for( int k=0;k<n;k++)
00099 {
00100 tensor::eval( ev[i*m+k] , mat[i*m+k].rep(),
00101 val , m );
00102 }
00103 }
00104 return ev;
00105 }
00106
00107 template<class T> inline
00108 int signof_det(T *mat,const int& n)
00109 {
00110 double ev;
00111 ev= det(mat,n);
00112
00113 if (ev< .1e-15) return (0);
00114 return (ev>0?1:-1);
00115 }
00116
00117 template<class T>
00118 T* submatrix(T *matrix, int order, int i, int j)
00119 {
00120 int n=order-1;
00121 T * subm;
00122 int p, q;
00123 int a = 0, b;
00124 subm = new T [n*n];
00125
00126 for(p = 0; p < order; p++) {
00127 if(p==i) continue;
00128
00129 b = 0;
00130
00131 for(q = 0; q < order; q++) {
00132 if(q==j) continue;
00133 subm[a*order + b++] = matrix[p*order+q];
00134 }
00135 a++;
00136 }
00137 return subm;
00138 }
00139
00140
00141 template<class T>
00142 T det(T *mat, int order)
00143 {
00144 if(order == 1)
00145 return mat[0];
00146 int i;
00147 T ev(0);
00148
00149 for(i = 0; i < order; i++)
00150 {
00151
00152 ev += T(i%2==0?1:-1)* mat[i*order] * det(submatrix(mat, order, i, 0), order - 1);
00153
00154
00155 }
00156 return ev;
00157 }
00158
00159
00160
00161 template<class FT>
00162 inline int signof(DPOL * p, Seq<FT> t)
00163 {
00164
00165 Interval<double> ev;
00166
00167
00168 tensor::eval( ev , p->rep() ,
00169 t , t.size() );
00170
00171 if (ev< .1e-10) return (0);
00172 return (ev>0?1:-1);
00173 }
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 TMPL
00193 bool inline miranda_test(Seq<POL>& S,const int i,const int j)
00194 {
00195 POL u,l;
00196
00197 tensor::face(l, S[i], j, 0);
00198 tensor::face(u, S[i], j, 1);
00199
00200
00201 return ( no_variation(l) &&
00202 no_variation(u) &&
00203 (l[0]>0)!=(u[0]>0) );
00204 };
00205
00207
00208 template<class POL>
00209 bool include1(BOX* b, DPOL * J)
00210 {
00211 Interval<double> ev(0);
00212 Interval<double> * m;
00213 unsigned i,j,c;
00214 POL u,l;
00215 Seq<POL> SS= b->system();
00216 unsigned dim= b->nbvar();
00217
00218 bool t[dim][dim];
00219
00220 m= eval_poly_matrix( b->template domain<double>() ,J,dim);
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238 ev=det(m,dim);
00239
00240 if ( ev.m*ev.M < 0 )
00241 return false;
00242
00243 for (i=0; i!=dim;++i)
00244 for (j=0; j!=dim;++j){
00245 t[i][j]= miranda_test(SS,i,j);
00246
00247 }
00248
00249 c=0;
00250 for (i=0; i!=dim;++i)
00251 for (j=0; j!=dim;++j)
00252 if (t[i][j]==true)
00253 {
00254 c++; break;
00255 }
00256 if (c<dim) return false;
00257
00258 c=0;
00259 for (i=0; i!=dim;++i)
00260 for (j=0; j!=dim;++j)
00261 if (t[j][i]==true)
00262 {
00263 c++; break;
00264 }
00265 if (c<dim) return false;
00266
00267
00268 return true;
00269 };
00270
00272 TMPL
00273 bool include2(BOX* b,DPOL * J)
00274 {
00275 Interval<double> ev;
00276 Interval<double> * m;
00277
00278
00279 unsigned dim= b->nbvar();
00280
00281
00282 m= eval_poly_matrix( b->template domain<double>() ,J,dim);
00283 ev= det(m,dim);
00284 if ( ev.m*ev.M < 0 )
00285 return false;
00286
00287 int td= topological_degree_2d<POL, double>(b);
00288
00289
00290
00291
00292
00293 return ( (td==-1 || td==1) );
00294 }
00295
00297 template<class POL>
00298 bool include3(BOX* b,DPOL * J, Seq<POL*>& S0)
00299 {
00300 double *fc, *jc, *ijc;
00301 Interval<double> ev, *jb, *m;
00302 Seq<double> c= b->point( b->middle() ) ;
00303 Seq<Interval<double> > bx;
00304 unsigned i,j,n= b->nbvar();
00305 DPOL tmp;
00306
00307 fc= new double[n];
00308 for(i=0; i< n;i++ )
00309 {
00310 let::assign(tmp, *S0[i]) ;
00311 tensor::eval( fc[i] , tmp.rep(), c , n );
00312 }
00313 bx= b->template domain<double>();
00314 jc= eval_poly_matrix( c, J, n);
00315 ijc= new double[n*n];
00316 linear::LUinverse(ijc, jc, n);
00317 jb= eval_poly_matrix( bx, J, n);
00318
00319 m= new Interval<double>[n*n];
00320 for(i=0; i< n;i++ )
00321 for(j=0; j< n;j++ )
00322 m[i*n+j] += ijc[i*n+j]*jb[j*n+i];
00323 free(jb);
00324
00325
00326
00327 for(i=0; i< n;i++ )
00328 {
00329 bx[i]= bx[i]- c[i];
00330 ev= bx[i];
00331 for(j=0; j<n; j++ )
00332 {
00333 ev+= -ijc[i*n+j]*fc[j] - m[i*n+j] *bx[i];
00334 }
00335
00336 if ( ev.m < bx[i].m ||
00337 ev.M > bx[i].M )
00338 {
00339
00340 return false;
00341 }
00342 }
00343
00344 return ( true );
00345 }
00346
00347
00349 template<class POL>
00350 bool exclude1(BOX*b, Seq<POL *>& S0)
00351 {
00352 Interval<double> ev;
00353 Seq<Interval<double> > dom;
00354 dom= b->template domain<double>();
00355 DPOL tmp;
00356
00357 for (unsigned i=0; i!=b->nbpol(); ++i)
00358 {
00359 let::assign(tmp, *S0[i] );
00360 tensor::eval( ev , tmp.rep(),
00361 dom , b->nbvar() );
00362 if ( ev.m*ev.M > 0 )
00363 {
00364
00365
00366
00367
00368
00369 return true;
00370 }
00371
00372 }
00373 return false;
00374 };
00375
00376
00378 TMPL
00379 bool exclude2(BOX* b,DPOL * J)
00380 {
00381 Interval<double> ev;
00382 Interval<double> * m;
00383
00384
00385 unsigned dim= b->nbvar();
00386
00387
00388 m= eval_poly_matrix( b->template domain<double>() ,J,dim);
00389 ev= det(m,dim);
00390 if ( ev.m*ev.M < 0 )
00391 return false;
00392
00393 int td= topological_degree_2d<POL, double>(b);
00394
00395 return ( td==0 );
00396 }
00397
00398 TMPL
00399 inline bool exclude3(BOX* b)
00400 {
00401 Seq<POL> S= b->system();
00402 for (unsigned i=0; i!=b->nbpol(); ++i)
00403 if ( no_variation( S[i]) )
00404 return true;
00405 return false;
00406 }
00407
00408
00409
00410
00411
00412
00414 template<class FT>
00415 inline int sgn(FT a )
00416 {
00417 if (a==FT(0) ) return 0;
00418 return ( (a>0) ? 1 : -1 );
00419 }
00420
00422 template<class POL, class FT>
00423 int inline sgn(POL & f, FT a )
00424 {
00425 if (a==FT(-1) )
00426 return sgn ( f[ f.rep().szs()[0]-1 ]);
00427 else
00428 return sgn ( f(a) );
00429 };
00430
00431
00432 template<class POL, class FT>
00433 int topological_degree_2d(BOX* b)
00434 {
00435
00436 typedef typename POL::scalar_t C;
00437
00438 assert( b->nbvar() ==2 );
00439 Seq<FT> p;
00440 Seq<FT> tmp;
00441 unsigned ev1(0), ev2(0), ev3(0), ev4(0);
00442 FT p1[2],p2[2];
00443 int tdeg;
00444 unsigned i;
00445
00446 POL fl0, fr0,fl1,fr1, s,
00447 gl0, gr0,gl1,gr1 ;
00448
00449 fl0=b->lface(0,1); gl0=b->lface(1,1);
00450 fr0=b->rface(0,1); gr0=b->rface(1,1);
00451 fl1=b->lface(0,0); gl1=b->lface(1,0);
00452 fr1=b->rface(0,0); gr1=b->rface(1,0);
00453
00454 s= fl0 * gl0 ;
00455 p<< FT(0);p << FT(0) ;
00456 tmp= solver<ring<C,MonomialTensor>, CFseparate>::template solve<FT>(s);
00457 for (i=0;i< tmp.size(); ++i )
00458 p<< tmp[i] << FT(0) ;
00459
00460 s= fr1 * gr1;
00461 p<< FT(-1);p << FT(0);
00462 tmp= solver<ring<C,MonomialTensor>, CFseparate>::template solve<FT>(s);
00463 for (i=0;i< tmp.size(); ++i )
00464 p<< FT(-1) << tmp[i];
00465
00466 s= fr0 * gr0 ;
00467 p<< FT(-1);p << FT(-1) ;
00468 tmp= solver<ring<C,MonomialTensor>, CFseparate>::template solve<FT>(s);
00469 for (i=0;i< tmp.size(); ++i )
00470 p<< tmp[i] << FT(-1) ;
00471
00472 s= fl1 * gl1 ;
00473 p<< FT(0);p << FT(-1) ;
00474 tmp= solver<ring<C,MonomialTensor>, CFseparate>::template solve<FT>(s);
00475 for (i=0;i< tmp.size(); ++i )
00476 p<< FT(0) << tmp[i] ;
00477
00478 p<< FT(0);p << FT(0) ;
00479
00480 tdeg= 0;
00481 for (i=0; i<p.size()-3; i+=2)
00482 {
00483 p1[0]= p[i] ;p1[1]=p[i+1];
00484 p2[0]= p[i+2];p2[1]=p[i+3];
00485
00486 if (p1[1]==FT(0) && p2[1]==FT(0)){
00487 ev1= sgn( fl0, p1[0] );
00488 ev3= sgn( fl0, p2[0] );
00489 ev4= sgn( gl0, p1[0] );
00490 ev2= sgn( gl0, p2[0] );
00491 }else
00492 if (p1[0]==FT(-1) && p2[0]==FT(-1)){
00493 ev1= sgn( fr1, p1[1] );
00494 ev3= sgn( fr1, p2[1] );
00495 ev4= sgn( gr1, p1[1] );
00496 ev2= sgn( gr1, p2[1] );
00497 }else
00498 if (p1[1]==FT(-1) && p2[1]==FT(-1)){
00499 ev1= sgn( fr0, p1[0] );
00500 ev3= sgn( fr0, p2[0] );
00501 ev4= sgn( gr0, p1[0] );
00502 ev2= sgn( gr0, p2[0] );
00503 }else
00504 if (p1[0]==FT(0) && p2[0]==FT(0)){
00505 ev1= sgn( fl1, p1[1] );
00506 ev3= sgn( fl1, p2[1] );
00507 ev4= sgn( gl1, p1[1] );
00508 ev2= sgn( gl1, p2[1] );
00509 }
00510 else
00511 std::cout<< "TROUBLE!"<<std::endl;
00512
00513 tdeg+= ( (ev1)*(ev2) - (ev3)*(ev4) );
00514 }
00515 return (tdeg/8);
00516 }
00517
00518
00519 }
00520 }
00521
00522 # undef DPOL
00523