[35aab3] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
| 4 | |
---|
[e2ad5d] | 5 | /* $Id: mpr_numeric.cc,v 1.8 2007-05-24 17:46:04 Singular Exp $ */ |
---|
[35aab3] | 6 | |
---|
| 7 | /* |
---|
| 8 | * ABSTRACT - multipolynomial resultants - numeric stuff |
---|
| 9 | * ( root finder, vandermonde system solver, simplex ) |
---|
| 10 | */ |
---|
| 11 | |
---|
| 12 | #include "mod2.h" |
---|
| 13 | //#ifdef HAVE_MPR |
---|
| 14 | |
---|
| 15 | //#define mprDEBUG_ALL |
---|
| 16 | |
---|
| 17 | //-> includes |
---|
| 18 | #include "structs.h" |
---|
| 19 | #include "febase.h" |
---|
| 20 | #include "omalloc.h" |
---|
| 21 | #include "numbers.h" |
---|
| 22 | #include "polys.h" |
---|
| 23 | #include "ideals.h" |
---|
| 24 | #include "intvec.h" |
---|
| 25 | #include "matpol.h" |
---|
| 26 | #include "ring.h" |
---|
| 27 | //#include "longrat.h" |
---|
| 28 | |
---|
| 29 | #include <math.h> |
---|
| 30 | #include "mpr_numeric.h" |
---|
| 31 | |
---|
| 32 | extern size_t gmp_output_digits; |
---|
| 33 | //<- |
---|
| 34 | |
---|
| 35 | extern void nPrint(number n); // for debugging output |
---|
| 36 | |
---|
| 37 | //----------------------------------------------------------------------------- |
---|
| 38 | //-------------- vandermonde system solver ------------------------------------ |
---|
| 39 | //----------------------------------------------------------------------------- |
---|
| 40 | |
---|
| 41 | //-> vandermonde::* |
---|
| 42 | vandermonde::vandermonde( const long _cn, const long _n, const long _maxdeg, |
---|
| 43 | number *_p, const bool _homog ) |
---|
| 44 | : n(_n), cn(_cn), maxdeg(_maxdeg), p(_p), homog(_homog) |
---|
| 45 | { |
---|
| 46 | long j; |
---|
| 47 | l= (long)pow((double)maxdeg+1,(int)n); |
---|
| 48 | x= (number *)omAlloc( cn * sizeof(number) ); |
---|
| 49 | for ( j= 0; j < cn; j++ ) x[j]= nInit(1); |
---|
| 50 | init(); |
---|
| 51 | } |
---|
| 52 | |
---|
| 53 | vandermonde::~vandermonde() |
---|
| 54 | { |
---|
| 55 | int j; |
---|
| 56 | for ( j= 0; j < cn; j++ ) nDelete( x+j ); |
---|
| 57 | omFreeSize( (ADDRESS)x, cn * sizeof( number ) ); |
---|
| 58 | } |
---|
| 59 | |
---|
| 60 | void vandermonde::init() |
---|
| 61 | { |
---|
| 62 | int j; |
---|
| 63 | long i,c,sum; |
---|
| 64 | number tmp,tmp1; |
---|
| 65 | |
---|
| 66 | c=0; |
---|
| 67 | sum=0; |
---|
| 68 | |
---|
| 69 | intvec exp( n ); |
---|
| 70 | for ( j= 0; j < n; j++ ) exp[j]=0; |
---|
| 71 | |
---|
| 72 | for ( i= 0; i < l; i++ ) |
---|
| 73 | { |
---|
| 74 | if ( !homog || (sum == maxdeg) ) |
---|
| 75 | { |
---|
| 76 | for ( j= 0; j < n; j++ ) |
---|
| 77 | { |
---|
| 78 | nPower( p[j], exp[j], &tmp ); |
---|
| 79 | tmp1 = nMult( tmp, x[c] ); |
---|
| 80 | x[c]= tmp1; |
---|
| 81 | nDelete( &tmp ); |
---|
| 82 | } |
---|
| 83 | c++; |
---|
| 84 | } |
---|
| 85 | exp[0]++; |
---|
| 86 | sum=0; |
---|
| 87 | for ( j= 0; j < n - 1; j++ ) |
---|
| 88 | { |
---|
| 89 | if ( exp[j] > maxdeg ) |
---|
| 90 | { |
---|
| 91 | exp[j]= 0; |
---|
| 92 | exp[j + 1]++; |
---|
| 93 | } |
---|
| 94 | sum+= exp[j]; |
---|
| 95 | } |
---|
| 96 | sum+= exp[n - 1]; |
---|
| 97 | } |
---|
| 98 | } |
---|
| 99 | |
---|
| 100 | poly vandermonde::numvec2poly( const number * q ) |
---|
| 101 | { |
---|
| 102 | int j; |
---|
| 103 | long i,c,sum; |
---|
| 104 | number tmp; |
---|
| 105 | |
---|
| 106 | poly pnew,pit=NULL; |
---|
| 107 | |
---|
| 108 | c=0; |
---|
| 109 | sum=0; |
---|
| 110 | |
---|
| 111 | int *exp= (int *) omAlloc( (n+1) * sizeof(int) ); |
---|
| 112 | |
---|
| 113 | for ( j= 0; j < n+1; j++ ) exp[j]=0; |
---|
| 114 | |
---|
| 115 | for ( i= 0; i < l; i++ ) |
---|
| 116 | { |
---|
| 117 | if ( (!homog || (sum == maxdeg)) && q[i] && !nIsZero(q[i]) ) |
---|
| 118 | { |
---|
| 119 | pnew= pOne(); |
---|
| 120 | pSetCoeff(pnew,q[i]); |
---|
| 121 | pSetExpV(pnew,exp); |
---|
| 122 | if ( pit ) |
---|
| 123 | { |
---|
| 124 | pNext(pnew)= pit; |
---|
| 125 | pit= pnew; |
---|
| 126 | } |
---|
| 127 | else |
---|
| 128 | { |
---|
| 129 | pit= pnew; |
---|
| 130 | pNext(pnew)= NULL; |
---|
| 131 | } |
---|
| 132 | pSetm(pit); |
---|
| 133 | } |
---|
| 134 | exp[1]++; |
---|
| 135 | sum=0; |
---|
| 136 | for ( j= 1; j < n; j++ ) |
---|
| 137 | { |
---|
| 138 | if ( exp[j] > maxdeg ) |
---|
| 139 | { |
---|
| 140 | exp[j]= 0; |
---|
| 141 | exp[j + 1]++; |
---|
| 142 | } |
---|
| 143 | sum+= exp[j]; |
---|
| 144 | } |
---|
| 145 | sum+= exp[n]; |
---|
| 146 | } |
---|
| 147 | |
---|
| 148 | omFreeSize( (ADDRESS) exp, (n+1) * sizeof(int) ); |
---|
| 149 | |
---|
| 150 | pSortAdd(pit); |
---|
| 151 | return pit; |
---|
| 152 | } |
---|
| 153 | |
---|
| 154 | number * vandermonde::interpolateDense( const number * q ) |
---|
| 155 | { |
---|
| 156 | int i,j,k; |
---|
| 157 | number newnum,tmp1; |
---|
| 158 | number b,t,xx,s; |
---|
| 159 | number *c; |
---|
| 160 | number *w; |
---|
| 161 | |
---|
| 162 | b=t=xx=s=tmp1=NULL; |
---|
| 163 | |
---|
| 164 | w= (number *)omAlloc( cn * sizeof(number) ); |
---|
| 165 | c= (number *)omAlloc( cn * sizeof(number) ); |
---|
| 166 | for ( j= 0; j < cn; j++ ) |
---|
| 167 | { |
---|
| 168 | w[j]= nInit(0); |
---|
| 169 | c[j]= nInit(0); |
---|
| 170 | } |
---|
| 171 | |
---|
| 172 | if ( cn == 1 ) |
---|
| 173 | { |
---|
| 174 | nDelete( &w[0] ); |
---|
| 175 | w[0]= nCopy(q[0]); |
---|
| 176 | } |
---|
| 177 | else |
---|
| 178 | { |
---|
| 179 | nDelete( &c[cn-1] ); |
---|
| 180 | c[cn-1]= nCopy(x[0]); |
---|
| 181 | c[cn-1]= nNeg(c[cn-1]); // c[cn]= -x[1] |
---|
| 182 | |
---|
| 183 | for ( i= 1; i < cn; i++ ) { // i=2; i <= cn |
---|
| 184 | nDelete( &xx ); |
---|
| 185 | xx= nCopy(x[i]); |
---|
| 186 | xx= nNeg(xx); // xx= -x[i] |
---|
| 187 | |
---|
| 188 | for ( j= (cn-i-1); j <= (cn-2); j++) { // j=(cn+1-i); j <= (cn-1) |
---|
| 189 | nDelete( &tmp1 ); |
---|
| 190 | tmp1= nMult( xx, c[j+1] ); // c[j]= c[j] + (xx * c[j+1]) |
---|
| 191 | newnum= nAdd( c[j], tmp1 ); |
---|
| 192 | nDelete( &c[j] ); |
---|
| 193 | c[j]= newnum; |
---|
| 194 | } |
---|
| 195 | |
---|
| 196 | newnum= nAdd( xx, c[cn-1] ); // c[cn-1]= c[cn-1] + xx |
---|
| 197 | nDelete( &c[cn-1] ); |
---|
| 198 | c[cn-1]= newnum; |
---|
| 199 | } |
---|
| 200 | |
---|
| 201 | for ( i= 0; i < cn; i++ ) { // i=1; i <= cn |
---|
| 202 | nDelete( &xx ); |
---|
| 203 | xx= nCopy(x[i]); // xx= x[i] |
---|
| 204 | |
---|
| 205 | nDelete( &t ); |
---|
| 206 | t= nInit( 1 ); // t= b= 1 |
---|
| 207 | nDelete( &b ); |
---|
| 208 | b= nInit( 1 ); |
---|
| 209 | nDelete( &s ); // s= q[cn-1] |
---|
| 210 | s= nCopy( q[cn-1] ); |
---|
| 211 | |
---|
| 212 | for ( k= cn-1; k >= 1; k-- ) { // k=cn; k >= 2 |
---|
| 213 | nDelete( &tmp1 ); |
---|
| 214 | tmp1= nMult( xx, b ); // b= c[k] + (xx * b) |
---|
| 215 | nDelete( &b ); |
---|
| 216 | b= nAdd( c[k], tmp1 ); |
---|
| 217 | |
---|
| 218 | nDelete( &tmp1 ); |
---|
| 219 | tmp1= nMult( q[k-1], b ); // s= s + (q[k-1] * b) |
---|
| 220 | newnum= nAdd( s, tmp1 ); |
---|
| 221 | nDelete( &s ); |
---|
| 222 | s= newnum; |
---|
| 223 | |
---|
| 224 | nDelete( &tmp1 ); |
---|
| 225 | tmp1= nMult( xx, t ); // t= (t * xx) + b |
---|
| 226 | newnum= nAdd( tmp1, b ); |
---|
| 227 | nDelete( &t ); |
---|
| 228 | t= newnum; |
---|
| 229 | } |
---|
| 230 | |
---|
| 231 | if (!nIsZero(t)) |
---|
| 232 | { |
---|
| 233 | nDelete( &w[i] ); // w[i]= s/t |
---|
| 234 | w[i]= nDiv( s, t ); |
---|
| 235 | nNormalize( w[i] ); |
---|
| 236 | } |
---|
| 237 | |
---|
| 238 | mprSTICKYPROT(ST_VANDER_STEP); |
---|
| 239 | } |
---|
| 240 | } |
---|
| 241 | mprSTICKYPROT("\n"); |
---|
| 242 | |
---|
| 243 | // free mem |
---|
| 244 | for ( j= 0; j < cn; j++ ) nDelete( c+j ); |
---|
| 245 | omFreeSize( (ADDRESS)c, cn * sizeof( number ) ); |
---|
| 246 | |
---|
| 247 | nDelete( &tmp1 ); |
---|
| 248 | nDelete( &s ); |
---|
| 249 | nDelete( &t ); |
---|
| 250 | nDelete( &b ); |
---|
| 251 | nDelete( &xx ); |
---|
| 252 | |
---|
| 253 | // makes quotiens smaller |
---|
| 254 | for ( j= 0; j < cn; j++ ) nNormalize( w[j] ); |
---|
| 255 | |
---|
| 256 | return w; |
---|
| 257 | } |
---|
| 258 | //<- |
---|
| 259 | |
---|
| 260 | //----------------------------------------------------------------------------- |
---|
| 261 | //-------------- rootContainer ------------------------------------------------ |
---|
| 262 | //----------------------------------------------------------------------------- |
---|
| 263 | |
---|
| 264 | //-> definitions |
---|
| 265 | #define MR 8 // never change this value |
---|
| 266 | #define MT 5 |
---|
| 267 | #define MMOD (MT*MR) |
---|
| 268 | #define MAXIT (5*MMOD) // max number of iterations in laguer root finder |
---|
| 269 | |
---|
| 270 | |
---|
| 271 | //-> rootContainer::rootContainer() |
---|
| 272 | rootContainer::rootContainer() |
---|
| 273 | { |
---|
| 274 | rt=none; |
---|
| 275 | |
---|
| 276 | coeffs= NULL; |
---|
| 277 | ievpoint= NULL; |
---|
| 278 | theroots= NULL; |
---|
| 279 | |
---|
| 280 | found_roots= false; |
---|
| 281 | } |
---|
| 282 | //<- |
---|
| 283 | |
---|
| 284 | //-> rootContainer::~rootContainer() |
---|
| 285 | rootContainer::~rootContainer() |
---|
| 286 | { |
---|
| 287 | int i; |
---|
| 288 | int n= pVariables; |
---|
| 289 | |
---|
| 290 | // free coeffs, ievpoint |
---|
| 291 | if ( ievpoint != NULL ) |
---|
| 292 | { |
---|
| 293 | for ( i=0; i < anz+2; i++ ) nDelete( ievpoint + i ); |
---|
| 294 | omFreeSize( (ADDRESS)ievpoint, (anz+2) * sizeof( number ) ); |
---|
| 295 | } |
---|
| 296 | |
---|
| 297 | for ( i=0; i <= tdg; i++ ) nDelete( coeffs + i ); |
---|
| 298 | omFreeSize( (ADDRESS)coeffs, (tdg+1) * sizeof( number ) ); |
---|
| 299 | |
---|
| 300 | // theroots löschen |
---|
| 301 | for ( i=0; i < tdg; i++ ) delete theroots[i]; |
---|
| 302 | omFreeSize( (ADDRESS) theroots, (tdg)*sizeof(gmp_complex*) ); |
---|
| 303 | |
---|
[cfbe751] | 304 | //mprPROTnl("~rootContainer()"); |
---|
[35aab3] | 305 | } |
---|
| 306 | //<- |
---|
| 307 | |
---|
| 308 | //-> void rootContainer::fillContainer( ... ) |
---|
| 309 | void rootContainer::fillContainer( number *_coeffs, number *_ievpoint, |
---|
| 310 | const int _var, const int _tdg, |
---|
| 311 | const rootType _rt, const int _anz ) |
---|
| 312 | { |
---|
| 313 | int i; |
---|
| 314 | number nn= nInit(0); |
---|
| 315 | var=_var; |
---|
| 316 | tdg=_tdg; |
---|
| 317 | coeffs=_coeffs; |
---|
| 318 | rt=_rt; |
---|
| 319 | anz=_anz; |
---|
| 320 | |
---|
| 321 | for ( i=0; i <= tdg; i++ ) |
---|
| 322 | { |
---|
| 323 | if ( nEqual(coeffs[i],nn) ) |
---|
| 324 | { |
---|
| 325 | nDelete( &coeffs[i] ); |
---|
| 326 | coeffs[i]=NULL; |
---|
| 327 | } |
---|
| 328 | } |
---|
| 329 | nDelete( &nn ); |
---|
| 330 | |
---|
[e2ad5d] | 331 | if ( rt == cspecialmu && _ievpoint ) // copy ievpoint |
---|
| 332 | { |
---|
[35aab3] | 333 | ievpoint= (number *)omAlloc( (anz+2) * sizeof( number ) ); |
---|
| 334 | for (i=0; i < anz+2; i++) ievpoint[i]= nCopy( _ievpoint[i] ); |
---|
| 335 | } |
---|
| 336 | |
---|
| 337 | theroots= NULL; |
---|
| 338 | found_roots= false; |
---|
| 339 | } |
---|
| 340 | //<- |
---|
| 341 | |
---|
| 342 | //-> poly rootContainer::getPoly() |
---|
| 343 | poly rootContainer::getPoly() |
---|
| 344 | { |
---|
| 345 | int i; |
---|
| 346 | |
---|
| 347 | poly result= NULL; |
---|
| 348 | poly ppos; |
---|
| 349 | |
---|
| 350 | if ( (rt == cspecial) || ( rt == cspecialmu ) ) |
---|
| 351 | { |
---|
| 352 | for ( i= tdg; i >= 0; i-- ) |
---|
| 353 | { |
---|
| 354 | if ( coeffs[i] ) |
---|
| 355 | { |
---|
| 356 | poly p= pOne(); |
---|
| 357 | //pSetExp( p, var+1, i); |
---|
| 358 | pSetExp( p, 1, i); |
---|
| 359 | pSetCoeff( p, nCopy( coeffs[i] ) ); |
---|
| 360 | pSetm( p ); |
---|
| 361 | if (result) |
---|
| 362 | { |
---|
| 363 | ppos->next=p; |
---|
| 364 | ppos=ppos->next; |
---|
| 365 | } |
---|
| 366 | else |
---|
| 367 | { |
---|
| 368 | result=p; |
---|
| 369 | ppos=p; |
---|
| 370 | } |
---|
| 371 | |
---|
| 372 | } |
---|
| 373 | } |
---|
[cfbe751] | 374 | if (result!=NULL) pSetm( result ); |
---|
[35aab3] | 375 | } |
---|
| 376 | |
---|
| 377 | return result; |
---|
| 378 | } |
---|
| 379 | //<- |
---|
| 380 | |
---|
| 381 | //-> const gmp_complex & rootContainer::opterator[] ( const int i ) |
---|
| 382 | // this is now inline! |
---|
| 383 | // gmp_complex & rootContainer::operator[] ( const int i ) |
---|
| 384 | // { |
---|
| 385 | // if ( found_roots && ( i >= 0) && ( i < tdg ) ) |
---|
| 386 | // { |
---|
| 387 | // return *theroots[i]; |
---|
| 388 | // } |
---|
| 389 | // // warning |
---|
| 390 | // Warn("rootContainer::getRoot: Wrong index %d, found_roots %s",i,found_roots?"true":"false"); |
---|
| 391 | // gmp_complex *tmp= new gmp_complex(); |
---|
| 392 | // return *tmp; |
---|
| 393 | // } |
---|
| 394 | //<- |
---|
| 395 | |
---|
| 396 | //-> gmp_complex & rootContainer::evPointCoord( int i ) |
---|
| 397 | gmp_complex & rootContainer::evPointCoord( const int i ) |
---|
| 398 | { |
---|
| 399 | if (! ((i >= 0) && (i < anz+2) ) ) |
---|
| 400 | WarnS("rootContainer::evPointCoord: index out of range"); |
---|
| 401 | if (ievpoint == NULL) |
---|
| 402 | WarnS("rootContainer::evPointCoord: ievpoint == NULL"); |
---|
| 403 | |
---|
[e2ad5d] | 404 | if ( (rt == cspecialmu) && found_roots ) // FIX ME |
---|
| 405 | { |
---|
[35aab3] | 406 | if ( ievpoint[i] != NULL ) |
---|
| 407 | { |
---|
| 408 | gmp_complex *tmp= new gmp_complex(); |
---|
| 409 | *tmp= numberToComplex(ievpoint[i]); |
---|
| 410 | return *tmp; |
---|
| 411 | } |
---|
| 412 | else |
---|
| 413 | { |
---|
| 414 | Warn("rootContainer::evPointCoord: NULL index %d",i); |
---|
| 415 | } |
---|
| 416 | } |
---|
| 417 | |
---|
| 418 | // warning |
---|
| 419 | Warn("rootContainer::evPointCoord: Wrong index %d, found_roots %s",i,found_roots?"true":"false"); |
---|
| 420 | gmp_complex *tmp= new gmp_complex(); |
---|
| 421 | return *tmp; |
---|
| 422 | } |
---|
| 423 | //<- |
---|
| 424 | |
---|
| 425 | //-> bool rootContainer::changeRoots( int from, int to ) |
---|
| 426 | bool rootContainer::swapRoots( const int from, const int to ) |
---|
| 427 | { |
---|
| 428 | if ( found_roots && ( from >= 0) && ( from < tdg ) && ( to >= 0) && ( to < tdg ) ) |
---|
| 429 | { |
---|
| 430 | if ( to != from ) |
---|
| 431 | { |
---|
| 432 | gmp_complex tmp( *theroots[from] ); |
---|
| 433 | *theroots[from]= *theroots[to]; |
---|
| 434 | *theroots[to]= tmp; |
---|
| 435 | } |
---|
| 436 | return true; |
---|
| 437 | } |
---|
| 438 | |
---|
| 439 | // warning |
---|
| 440 | Warn(" rootContainer::changeRoots: Wrong index %d, %d",from,to); |
---|
| 441 | return false; |
---|
| 442 | } |
---|
| 443 | //<- |
---|
| 444 | |
---|
| 445 | //-> void rootContainer::solver() |
---|
| 446 | bool rootContainer::solver( const int polishmode ) |
---|
| 447 | { |
---|
| 448 | int i; |
---|
| 449 | |
---|
| 450 | // there are maximal tdg roots, so *roots ranges form 0 to tdg-1. |
---|
| 451 | theroots= (gmp_complex**)omAlloc( tdg*sizeof(gmp_complex*) ); |
---|
| 452 | for ( i=0; i < tdg; i++ ) theroots[i]= new gmp_complex(); |
---|
| 453 | |
---|
| 454 | // copy the coefficients of type number to type gmp_complex |
---|
| 455 | gmp_complex **ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) ); |
---|
| 456 | for ( i=0; i <= tdg; i++ ) |
---|
| 457 | { |
---|
| 458 | ad[i]= new gmp_complex(); |
---|
| 459 | if ( coeffs[i] ) *ad[i] = numberToComplex( coeffs[i] ); |
---|
| 460 | } |
---|
| 461 | |
---|
| 462 | // now solve |
---|
| 463 | found_roots= laguer_driver( ad, theroots, polishmode != 0 ); |
---|
| 464 | if (!found_roots) |
---|
| 465 | WarnS("rootContainer::solver: No roots found!"); |
---|
| 466 | |
---|
| 467 | // free memory |
---|
| 468 | for ( i=0; i <= tdg; i++ ) delete ad[i]; |
---|
| 469 | omFreeSize( (ADDRESS) ad, (tdg+1)*sizeof(gmp_complex*) ); |
---|
| 470 | |
---|
| 471 | return found_roots; |
---|
| 472 | } |
---|
| 473 | //<- |
---|
| 474 | |
---|
| 475 | //-> gmp_complex* rootContainer::laguer_driver( bool polish ) |
---|
| 476 | bool rootContainer::laguer_driver(gmp_complex ** a, gmp_complex ** roots, bool polish ) |
---|
| 477 | { |
---|
| 478 | int i,j,k,its; |
---|
| 479 | gmp_float zero(0.0); |
---|
| 480 | gmp_complex x(0.0),o(1.0); |
---|
| 481 | bool ret= true, isf=isfloat(a), type=true; |
---|
| 482 | |
---|
| 483 | gmp_complex ** ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) ); |
---|
| 484 | for ( i=0; i <= tdg; i++ ) ad[i]= new gmp_complex( *a[i] ); |
---|
| 485 | |
---|
| 486 | k = 0; |
---|
| 487 | i = tdg; |
---|
| 488 | j = i-1; |
---|
| 489 | while (i>2) |
---|
| 490 | { |
---|
| 491 | // run laguer alg |
---|
| 492 | x = zero; |
---|
| 493 | laguer(ad, i, &x, &its, type); |
---|
| 494 | if ( its > MAXIT ) |
---|
| 495 | { |
---|
| 496 | type = !type; |
---|
| 497 | x = zero; |
---|
| 498 | laguer(ad, i, &x, &its, type); |
---|
| 499 | } |
---|
| 500 | |
---|
| 501 | mprSTICKYPROT(ST_ROOTS_LGSTEP); |
---|
| 502 | if ( its > MAXIT ) |
---|
| 503 | { // error |
---|
[695341] | 504 | WarnS("Laguerre solver: Too many iterations!"); |
---|
[35aab3] | 505 | ret= false; |
---|
| 506 | goto theend; |
---|
| 507 | } |
---|
| 508 | if ( polish ) |
---|
| 509 | { |
---|
| 510 | laguer( a, tdg, &x, &its , type); |
---|
| 511 | if ( its > MAXIT ) |
---|
| 512 | { // error |
---|
[695341] | 513 | WarnS("Laguerre solver: Too many iterations in polish!"); |
---|
[35aab3] | 514 | ret= false; |
---|
| 515 | goto theend; |
---|
| 516 | } |
---|
| 517 | } |
---|
[cfbe751] | 518 | if ((!type)&&(!((x.real()==zero)&&(x.imag()==zero)))) x = o/x; |
---|
[35aab3] | 519 | if (x.imag() == zero) |
---|
| 520 | { |
---|
| 521 | *roots[k] = x; |
---|
| 522 | k++; |
---|
| 523 | divlin(ad,x,i); |
---|
| 524 | i--; |
---|
| 525 | } |
---|
| 526 | else |
---|
| 527 | { |
---|
| 528 | if(isf) |
---|
| 529 | { |
---|
| 530 | *roots[j] = x; |
---|
| 531 | *roots[j-1]= gmp_complex(x.real(),-x.imag()); |
---|
| 532 | j -= 2; |
---|
| 533 | divquad(ad,x,i); |
---|
| 534 | i -= 2; |
---|
| 535 | } |
---|
| 536 | else |
---|
| 537 | { |
---|
| 538 | *roots[j] = x; |
---|
| 539 | j--; |
---|
| 540 | divlin(ad,x,i); |
---|
| 541 | i--; |
---|
| 542 | } |
---|
| 543 | } |
---|
| 544 | type = !type; |
---|
| 545 | } |
---|
| 546 | solvequad(ad,roots,k,j); |
---|
| 547 | sortroots(roots,k,j,isf); |
---|
| 548 | |
---|
| 549 | theend: |
---|
| 550 | mprSTICKYPROT("\n"); |
---|
| 551 | for ( i=0; i <= tdg; i++ ) delete ad[i]; |
---|
| 552 | omFreeSize( (ADDRESS) ad, (tdg+1)*sizeof( gmp_complex* )); |
---|
| 553 | |
---|
| 554 | return ret; |
---|
| 555 | } |
---|
| 556 | //<- |
---|
| 557 | |
---|
| 558 | //-> void rootContainer::laguer(...) |
---|
| 559 | void rootContainer::laguer(gmp_complex ** a, int m, gmp_complex *x, int *its, bool type) |
---|
| 560 | { |
---|
| 561 | int iter,j; |
---|
| 562 | gmp_float zero(0.0),one(1.0),deg(m); |
---|
| 563 | gmp_float abx_g, err_g, fabs; |
---|
| 564 | gmp_complex dx, x1, b, d, f, g, h, sq, gp, gm, g2; |
---|
| 565 | gmp_float frac_g[MR+1] = { 0.0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0 }; |
---|
| 566 | |
---|
| 567 | gmp_float epss(0.1); |
---|
| 568 | mpf_pow_ui(*epss._mpfp(),*epss.mpfp(),getGMPFloatDigits()); |
---|
| 569 | |
---|
| 570 | for ( iter= 1; iter <= MAXIT; iter++ ) |
---|
| 571 | { |
---|
| 572 | mprSTICKYPROT(ST_ROOTS_LG); |
---|
| 573 | *its=iter; |
---|
| 574 | if (type) |
---|
| 575 | computefx(a,*x,m,b,d,f,abx_g,err_g); |
---|
| 576 | else |
---|
| 577 | computegx(a,*x,m,b,d,f,abx_g,err_g); |
---|
| 578 | err_g *= epss; // EPSS; |
---|
| 579 | |
---|
| 580 | fabs = abs(b); |
---|
| 581 | if (fabs <= err_g) |
---|
| 582 | { |
---|
| 583 | if ((fabs==zero) || (abs(d)==zero)) return; |
---|
| 584 | *x -= (b/d); // a last newton-step |
---|
| 585 | goto ende; |
---|
| 586 | } |
---|
| 587 | |
---|
| 588 | g= d / b; |
---|
| 589 | g2 = g * g; |
---|
| 590 | h= g2 - (((f+f) / b )); |
---|
| 591 | sq= sqrt(( ( h * deg ) - g2 ) * (deg - one)); |
---|
| 592 | gp= g + sq; |
---|
| 593 | gm= g - sq; |
---|
| 594 | if (abs(gp)<abs(gm)) |
---|
| 595 | { |
---|
| 596 | dx = deg/gm; |
---|
| 597 | } |
---|
| 598 | else |
---|
| 599 | { |
---|
| 600 | if((gp.real()==zero)&&(gp.imag()==zero)) |
---|
| 601 | { |
---|
| 602 | dx.real(cos((mprfloat)iter)); |
---|
| 603 | dx.imag(sin((mprfloat)iter)); |
---|
| 604 | dx = dx*(one+abx_g); |
---|
| 605 | } |
---|
| 606 | else |
---|
| 607 | { |
---|
| 608 | dx = deg/gp; |
---|
| 609 | } |
---|
| 610 | } |
---|
| 611 | x1= *x - dx; |
---|
| 612 | |
---|
| 613 | if (*x == x1) goto ende; |
---|
| 614 | |
---|
| 615 | j = iter%MMOD; |
---|
| 616 | if (j==0) j=MT; |
---|
| 617 | if ( j % MT ) *x= x1; |
---|
| 618 | else *x -= ( dx * frac_g[ j / MT ] ); |
---|
| 619 | } |
---|
| 620 | |
---|
| 621 | *its= MAXIT+1; |
---|
| 622 | ende: |
---|
| 623 | checkimag(x,epss); |
---|
| 624 | } |
---|
| 625 | |
---|
| 626 | void rootContainer::checkimag(gmp_complex *x, gmp_float &e) |
---|
| 627 | { |
---|
| 628 | if(abs(x->imag())<abs(x->real())*e) |
---|
| 629 | { |
---|
| 630 | x->imag(0.0); |
---|
| 631 | } |
---|
| 632 | } |
---|
| 633 | |
---|
| 634 | bool rootContainer::isfloat(gmp_complex **a) |
---|
| 635 | { |
---|
| 636 | gmp_float z(0.0); |
---|
| 637 | gmp_complex *b; |
---|
| 638 | for (int i=tdg; i >= 0; i-- ) |
---|
| 639 | { |
---|
| 640 | b = &(*a[i]); |
---|
| 641 | if (!(b->imag()==z)) |
---|
| 642 | return false; |
---|
| 643 | } |
---|
| 644 | return true; |
---|
| 645 | } |
---|
| 646 | |
---|
| 647 | void rootContainer::divlin(gmp_complex **a, gmp_complex x, int j) |
---|
| 648 | { |
---|
| 649 | int i; |
---|
| 650 | gmp_float o(1.0); |
---|
| 651 | |
---|
| 652 | if (abs(x)<o) |
---|
| 653 | { |
---|
| 654 | for (i= j-1; i > 0; i-- ) |
---|
| 655 | *a[i] += (*a[i+1]*x); |
---|
| 656 | for (i= 0; i < j; i++ ) |
---|
| 657 | *a[i] = *a[i+1]; |
---|
| 658 | } |
---|
| 659 | else |
---|
| 660 | { |
---|
| 661 | gmp_complex y(o/x); |
---|
| 662 | for (i= 1; i < j; i++) |
---|
| 663 | *a[i] += (*a[i-1]*y); |
---|
| 664 | } |
---|
| 665 | } |
---|
| 666 | |
---|
| 667 | void rootContainer::divquad(gmp_complex **a, gmp_complex x, int j) |
---|
| 668 | { |
---|
| 669 | int i; |
---|
| 670 | gmp_float o(1.0),p(x.real()+x.real()), |
---|
| 671 | q((x.real()*x.real())+(x.imag()*x.imag())); |
---|
| 672 | |
---|
| 673 | if (abs(x)<o) |
---|
| 674 | { |
---|
| 675 | *a[j-1] += (*a[j]*p); |
---|
| 676 | for (i= j-2; i > 1; i-- ) |
---|
| 677 | *a[i] += ((*a[i+1]*p)-(*a[i+2]*q)); |
---|
| 678 | for (i= 0; i < j-1; i++ ) |
---|
| 679 | *a[i] = *a[i+2]; |
---|
| 680 | } |
---|
| 681 | else |
---|
| 682 | { |
---|
| 683 | p = p/q; |
---|
| 684 | q = o/q; |
---|
| 685 | *a[1] += (*a[0]*p); |
---|
| 686 | for (i= 2; i < j-1; i++) |
---|
| 687 | *a[i] += ((*a[i-1]*p)-(*a[i-2]*q)); |
---|
| 688 | } |
---|
| 689 | } |
---|
| 690 | |
---|
| 691 | void rootContainer::solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j) |
---|
| 692 | { |
---|
| 693 | gmp_float zero(0.0); |
---|
| 694 | |
---|
[cfbe751] | 695 | if ((j>k) |
---|
[c30c46] | 696 | &&((!(*a[2]).real().isZero())||(!(*a[2]).imag().isZero()))) |
---|
[35aab3] | 697 | { |
---|
| 698 | gmp_complex sq(zero); |
---|
| 699 | gmp_complex h1(*a[1]/(*a[2] + *a[2])), h2(*a[0] / *a[2]); |
---|
| 700 | gmp_complex disk((h1 * h1) - h2); |
---|
[bcb200] | 701 | if (disk.imag().isZero()) |
---|
[35aab3] | 702 | { |
---|
| 703 | if (disk.real()<zero) |
---|
| 704 | { |
---|
| 705 | sq.real(zero); |
---|
| 706 | sq.imag(sqrt(-disk.real())); |
---|
| 707 | } |
---|
| 708 | else |
---|
| 709 | sq = (gmp_complex)sqrt(disk.real()); |
---|
| 710 | } |
---|
| 711 | else |
---|
| 712 | sq = sqrt(disk); |
---|
| 713 | *r[k+1] = sq - h1; |
---|
| 714 | sq += h1; |
---|
| 715 | *r[k] = (gmp_complex)0.0-sq; |
---|
[bcb200] | 716 | if(sq.imag().isZero()) |
---|
[35aab3] | 717 | { |
---|
| 718 | k = j; |
---|
| 719 | j++; |
---|
| 720 | } |
---|
| 721 | else |
---|
| 722 | { |
---|
| 723 | j = k; |
---|
| 724 | k--; |
---|
| 725 | } |
---|
| 726 | } |
---|
| 727 | else |
---|
| 728 | { |
---|
[bcb200] | 729 | if (((*a[1]).real().isZero()) && ((*a[1]).imag().isZero())) |
---|
[cfbe751] | 730 | { |
---|
| 731 | WerrorS("precision lost, try again with higher precision"); |
---|
| 732 | } |
---|
[35aab3] | 733 | else |
---|
[cfbe751] | 734 | { |
---|
| 735 | *r[k]= (gmp_complex)0.0-(*a[0] / *a[1]); |
---|
[bcb200] | 736 | if(r[k]->imag().isZero()) |
---|
[cfbe751] | 737 | j++; |
---|
| 738 | else |
---|
| 739 | k--; |
---|
| 740 | } |
---|
[35aab3] | 741 | } |
---|
| 742 | } |
---|
| 743 | |
---|
| 744 | void rootContainer::sortroots(gmp_complex **ro, int r, int c, bool isf) |
---|
| 745 | { |
---|
| 746 | int j; |
---|
| 747 | |
---|
| 748 | for (j=0; j<r; j++) // the real roots |
---|
| 749 | sortre(ro, j, r, 1); |
---|
| 750 | if (c>=tdg) return; |
---|
| 751 | if (isf) |
---|
| 752 | { |
---|
| 753 | for (j=c; j+2<tdg; j+=2) // the complex roots for a real poly |
---|
| 754 | sortre(ro, j, tdg-1, 2); |
---|
| 755 | } |
---|
| 756 | else |
---|
| 757 | { |
---|
| 758 | for (j=c; j+1<tdg; j++) // the complex roots for a general poly |
---|
| 759 | sortre(ro, j, tdg-1, 1); |
---|
| 760 | } |
---|
| 761 | } |
---|
| 762 | |
---|
| 763 | void rootContainer::sortre(gmp_complex **r, int l, int u, int inc) |
---|
| 764 | { |
---|
| 765 | int pos,i; |
---|
| 766 | gmp_complex *x,*y; |
---|
| 767 | |
---|
| 768 | pos = l; |
---|
| 769 | x = r[pos]; |
---|
| 770 | for (i=l+inc; i<=u; i+=inc) |
---|
| 771 | { |
---|
| 772 | if (r[i]->real()<x->real()) |
---|
| 773 | { |
---|
| 774 | pos = i; |
---|
| 775 | x = r[pos]; |
---|
| 776 | } |
---|
| 777 | } |
---|
| 778 | if (pos>l) |
---|
| 779 | { |
---|
| 780 | if (inc==1) |
---|
| 781 | { |
---|
| 782 | for (i=pos; i>l; i--) |
---|
| 783 | r[i] = r[i-1]; |
---|
| 784 | r[l] = x; |
---|
| 785 | } |
---|
| 786 | else |
---|
| 787 | { |
---|
| 788 | y = r[pos+1]; |
---|
| 789 | for (i=pos+1; i+1>l; i--) |
---|
| 790 | r[i] = r[i-2]; |
---|
| 791 | if (x->imag()>y->imag()) |
---|
| 792 | { |
---|
| 793 | r[l] = x; |
---|
| 794 | r[l+1] = y; |
---|
| 795 | } |
---|
| 796 | else |
---|
| 797 | { |
---|
| 798 | r[l] = y; |
---|
| 799 | r[l+1] = x; |
---|
| 800 | } |
---|
| 801 | } |
---|
| 802 | } |
---|
| 803 | else if ((inc==2)&&(x->imag()<r[l+1]->imag())) |
---|
| 804 | { |
---|
| 805 | r[l] = r[l+1]; |
---|
| 806 | r[l+1] = x; |
---|
| 807 | } |
---|
| 808 | } |
---|
| 809 | |
---|
| 810 | void rootContainer::computefx(gmp_complex **a, gmp_complex x, int m, |
---|
| 811 | gmp_complex &f0, gmp_complex &f1, gmp_complex &f2, |
---|
| 812 | gmp_float &ex, gmp_float &ef) |
---|
| 813 | { |
---|
| 814 | int k; |
---|
| 815 | |
---|
| 816 | f0= *a[m]; |
---|
| 817 | ef= abs(f0); |
---|
| 818 | f1= gmp_complex(0.0); |
---|
| 819 | f2= f1; |
---|
| 820 | ex= abs(x); |
---|
| 821 | |
---|
| 822 | for ( k= m-1; k >= 0; k-- ) |
---|
| 823 | { |
---|
| 824 | f2 = ( x * f2 ) + f1; |
---|
| 825 | f1 = ( x * f1 ) + f0; |
---|
| 826 | f0 = ( x * f0 ) + *a[k]; |
---|
| 827 | ef = abs( f0 ) + ( ex * ef ); |
---|
| 828 | } |
---|
| 829 | } |
---|
| 830 | |
---|
| 831 | void rootContainer::computegx(gmp_complex **a, gmp_complex x, int m, |
---|
| 832 | gmp_complex &f0, gmp_complex &f1, gmp_complex &f2, |
---|
| 833 | gmp_float &ex, gmp_float &ef) |
---|
| 834 | { |
---|
| 835 | int k; |
---|
| 836 | |
---|
| 837 | f0= *a[0]; |
---|
| 838 | ef= abs(f0); |
---|
| 839 | f1= gmp_complex(0.0); |
---|
| 840 | f2= f1; |
---|
| 841 | ex= abs(x); |
---|
| 842 | |
---|
| 843 | for ( k= 1; k <= m; k++ ) |
---|
| 844 | { |
---|
| 845 | f2 = ( x * f2 ) + f1; |
---|
| 846 | f1 = ( x * f1 ) + f0; |
---|
| 847 | f0 = ( x * f0 ) + *a[k]; |
---|
| 848 | ef = abs( f0 ) + ( ex * ef ); |
---|
| 849 | } |
---|
| 850 | } |
---|
| 851 | |
---|
| 852 | //----------------------------------------------------------------------------- |
---|
| 853 | //-------------- rootArranger ------------------------------------------------- |
---|
| 854 | //----------------------------------------------------------------------------- |
---|
| 855 | |
---|
| 856 | //-> rootArranger::rootArranger(...) |
---|
| 857 | rootArranger::rootArranger( rootContainer ** _roots, |
---|
| 858 | rootContainer ** _mu, |
---|
| 859 | const int _howclean ) |
---|
| 860 | : roots(_roots), mu(_mu), howclean(_howclean) |
---|
| 861 | { |
---|
| 862 | found_roots=false; |
---|
| 863 | } |
---|
| 864 | //<- |
---|
| 865 | |
---|
| 866 | //-> void rootArranger::solve_all() |
---|
| 867 | void rootArranger::solve_all() |
---|
| 868 | { |
---|
| 869 | int i; |
---|
| 870 | found_roots= true; |
---|
| 871 | |
---|
| 872 | // find roots of polys given by coeffs in roots |
---|
| 873 | rc= roots[0]->getAnzElems(); |
---|
| 874 | for ( i= 0; i < rc; i++ ) |
---|
| 875 | if ( !roots[i]->solver( howclean ) ) |
---|
| 876 | { |
---|
| 877 | found_roots= false; |
---|
| 878 | return; |
---|
| 879 | } |
---|
| 880 | // find roots of polys given by coeffs in mu |
---|
| 881 | mc= mu[0]->getAnzElems(); |
---|
| 882 | for ( i= 0; i < mc; i++ ) |
---|
| 883 | if ( ! mu[i]->solver( howclean ) ) |
---|
| 884 | { |
---|
| 885 | found_roots= false; |
---|
| 886 | return; |
---|
| 887 | } |
---|
| 888 | } |
---|
| 889 | //<- |
---|
| 890 | |
---|
| 891 | //-> void rootArranger::arrange() |
---|
| 892 | void rootArranger::arrange() |
---|
| 893 | { |
---|
| 894 | gmp_complex tmp,zwerg; |
---|
| 895 | int anzm= mu[0]->getAnzElems(); |
---|
| 896 | int anzr= roots[0]->getAnzRoots(); |
---|
| 897 | int xkoord, r, rtest, xk, mtest; |
---|
| 898 | bool found; |
---|
| 899 | //gmp_complex mprec(1.0/pow(10,gmp_output_digits-5),1.0/pow(10,gmp_output_digits-5)); |
---|
| 900 | |
---|
| 901 | for ( xkoord= 0; xkoord < anzm; xkoord++ ) { // für x1,x2, x1,x2,x3, x1,x2,...,xn |
---|
[49698f3] | 902 | gmp_float mprec(1.0/pow(10.0,(int)(gmp_output_digits/3))); |
---|
[35aab3] | 903 | for ( r= 0; r < anzr; r++ ) { // für jede Nullstelle |
---|
| 904 | // (x1-koordinate) * evp[1] + (x2-koordinate) * evp[2] + |
---|
| 905 | // ... + (xkoord-koordinate) * evp[xkoord] |
---|
| 906 | tmp= gmp_complex(); |
---|
| 907 | for ( xk =0; xk <= xkoord; xk++ ) |
---|
| 908 | { |
---|
| 909 | tmp -= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1); //xk+1 |
---|
| 910 | } |
---|
| 911 | found= false; |
---|
[49698f3] | 912 | do { // while not found |
---|
| 913 | for ( rtest= r; rtest < anzr; rtest++ ) { // für jede Nullstelle |
---|
| 914 | zwerg = tmp - (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xk+1, xkoord+2 |
---|
| 915 | for ( mtest= 0; mtest < anzr; mtest++ ) |
---|
| 916 | { |
---|
| 917 | // if ( tmp == (*mu[xkoord])[mtest] ) |
---|
| 918 | // { |
---|
| 919 | if ( ((zwerg.real() <= (*mu[xkoord])[mtest].real() + mprec) && |
---|
| 920 | (zwerg.real() >= (*mu[xkoord])[mtest].real() - mprec)) && |
---|
| 921 | ((zwerg.imag() <= (*mu[xkoord])[mtest].imag() + mprec) && |
---|
| 922 | (zwerg.imag() >= (*mu[xkoord])[mtest].imag() - mprec)) ) |
---|
| 923 | { |
---|
| 924 | roots[xk]->swapRoots( r, rtest ); |
---|
| 925 | found= true; |
---|
| 926 | break; |
---|
| 927 | } |
---|
| 928 | } |
---|
| 929 | } // rtest |
---|
| 930 | if (!found) |
---|
[35aab3] | 931 | { |
---|
[49698f3] | 932 | WarnS("rootArranger::arrange: precision lost"); |
---|
| 933 | mprec*=10; |
---|
| 934 | } |
---|
| 935 | } while(!found); |
---|
| 936 | #if 0 |
---|
[35aab3] | 937 | if ( !found ) |
---|
| 938 | { |
---|
| 939 | Warn("rootArranger::arrange: No match? coord %d, root %d.",xkoord,r); |
---|
[49698f3] | 940 | //#ifdef mprDEBUG_PROT |
---|
[35aab3] | 941 | WarnS("One of these ..."); |
---|
| 942 | for ( rtest= r; rtest < anzr; rtest++ ) |
---|
| 943 | { |
---|
| 944 | tmp= gmp_complex(); |
---|
| 945 | for ( xk =0; xk <= xkoord; xk++ ) |
---|
| 946 | { |
---|
| 947 | tmp-= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1); |
---|
| 948 | } |
---|
| 949 | tmp-= (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xkoord+2 |
---|
| 950 | Warn(" %s",complexToStr(tmp,gmp_output_digits+1),rtest); |
---|
| 951 | } |
---|
| 952 | WarnS(" ... must match to one of these:"); |
---|
| 953 | for ( mtest= 0; mtest < anzr; mtest++ ) |
---|
| 954 | { |
---|
| 955 | Warn(" %s",complexToStr((*mu[xkoord])[mtest],gmp_output_digits+1)); |
---|
| 956 | } |
---|
[49698f3] | 957 | //#endif |
---|
[35aab3] | 958 | } |
---|
[49698f3] | 959 | #endif |
---|
[35aab3] | 960 | } // r |
---|
| 961 | } // xkoord |
---|
| 962 | } |
---|
| 963 | //<- |
---|
| 964 | |
---|
| 965 | //----------------------------------------------------------------------------- |
---|
| 966 | //-------------- simplex ----- ------------------------------------------------ |
---|
| 967 | //----------------------------------------------------------------------------- |
---|
| 968 | |
---|
| 969 | // #ifdef mprDEBUG_PROT |
---|
| 970 | // #define error(a) a |
---|
| 971 | // #else |
---|
| 972 | // #define error(a) |
---|
| 973 | // #endif |
---|
| 974 | |
---|
| 975 | #define error(a) a |
---|
| 976 | |
---|
| 977 | #define MAXPOINTS 1000 |
---|
| 978 | |
---|
| 979 | //-> simplex::* |
---|
| 980 | // |
---|
| 981 | simplex::simplex( int rows, int cols ) |
---|
| 982 | : LiPM_cols(cols), LiPM_rows(rows) |
---|
| 983 | { |
---|
| 984 | int i; |
---|
| 985 | |
---|
| 986 | LiPM_rows=LiPM_rows+3; |
---|
| 987 | LiPM_cols=LiPM_cols+2; |
---|
| 988 | |
---|
| 989 | LiPM = (mprfloat **)omAlloc( LiPM_rows * sizeof(mprfloat *) ); // LP matrix |
---|
| 990 | for( i= 0; i < LiPM_rows; i++ ) |
---|
| 991 | { |
---|
| 992 | // Mem must be allocated aligned, also for type double! |
---|
| 993 | LiPM[i] = (mprfloat *)omAlloc0Aligned( LiPM_cols * sizeof(mprfloat) ); |
---|
| 994 | } |
---|
| 995 | |
---|
| 996 | iposv = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) ); |
---|
| 997 | izrov = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) ); |
---|
| 998 | |
---|
| 999 | m=n=m1=m2=m3=icase=0; |
---|
| 1000 | |
---|
| 1001 | #ifdef mprDEBUG_ALL |
---|
| 1002 | Print("LiPM size: %d, %d\n",LiPM_rows,LiPM_cols); |
---|
| 1003 | #endif |
---|
| 1004 | } |
---|
| 1005 | |
---|
| 1006 | simplex::~simplex() |
---|
| 1007 | { |
---|
| 1008 | // clean up |
---|
| 1009 | int i; |
---|
| 1010 | for( i= 0; i < LiPM_rows; i++ ) |
---|
| 1011 | { |
---|
| 1012 | omFreeSize( (ADDRESS) LiPM[i], LiPM_cols * sizeof(mprfloat) ); |
---|
| 1013 | } |
---|
| 1014 | omFreeSize( (ADDRESS) LiPM, LiPM_rows * sizeof(mprfloat *) ); |
---|
| 1015 | |
---|
| 1016 | omFreeSize( (ADDRESS) iposv, 2*LiPM_rows*sizeof(int) ); |
---|
| 1017 | omFreeSize( (ADDRESS) izrov, 2*LiPM_rows*sizeof(int) ); |
---|
| 1018 | } |
---|
| 1019 | |
---|
| 1020 | BOOLEAN simplex::mapFromMatrix( matrix m ) |
---|
| 1021 | { |
---|
| 1022 | int i,j; |
---|
| 1023 | // if ( MATROWS( m ) > LiPM_rows || MATCOLS( m ) > LiPM_cols ) { |
---|
| 1024 | // WarnS(""); |
---|
| 1025 | // return FALSE; |
---|
| 1026 | // } |
---|
| 1027 | |
---|
| 1028 | number coef; |
---|
| 1029 | for ( i= 1; i <= MATROWS( m ); i++ ) |
---|
| 1030 | { |
---|
| 1031 | for ( j= 1; j <= MATCOLS( m ); j++ ) |
---|
| 1032 | { |
---|
| 1033 | if ( MATELEM(m,i,j) != NULL ) |
---|
| 1034 | { |
---|
| 1035 | coef= pGetCoeff( MATELEM(m,i,j) ); |
---|
| 1036 | if ( coef != NULL && !nIsZero(coef) ) |
---|
| 1037 | LiPM[i][j]= (double)(*(gmp_float*)coef); |
---|
| 1038 | //#ifdef mpr_DEBUG_PROT |
---|
| 1039 | //Print("%f ",LiPM[i][j]); |
---|
| 1040 | //#endif |
---|
| 1041 | } |
---|
| 1042 | } |
---|
| 1043 | // PrintLn(); |
---|
| 1044 | } |
---|
| 1045 | |
---|
| 1046 | return TRUE; |
---|
| 1047 | } |
---|
| 1048 | |
---|
| 1049 | matrix simplex::mapToMatrix( matrix m ) |
---|
| 1050 | { |
---|
| 1051 | int i,j; |
---|
| 1052 | // if ( MATROWS( m ) < LiPM_rows-3 || MATCOLS( m ) < LiPM_cols-2 ) { |
---|
| 1053 | // WarnS(""); |
---|
| 1054 | // return NULL; |
---|
| 1055 | // } |
---|
| 1056 | |
---|
| 1057 | //Print(" %d x %d\n",MATROWS( m ),MATCOLS( m )); |
---|
| 1058 | |
---|
| 1059 | number coef; |
---|
| 1060 | gmp_float * bla; |
---|
| 1061 | for ( i= 1; i <= MATROWS( m ); i++ ) |
---|
| 1062 | { |
---|
| 1063 | for ( j= 1; j <= MATCOLS( m ); j++ ) |
---|
| 1064 | { |
---|
| 1065 | pDelete( &(MATELEM(m,i,j)) ); |
---|
| 1066 | MATELEM(m,i,j)= NULL; |
---|
| 1067 | //Print(" %3.0f ",LiPM[i][j]); |
---|
| 1068 | if ( LiPM[i][j] != 0.0 ) |
---|
| 1069 | { |
---|
| 1070 | bla= new gmp_float(LiPM[i][j]); |
---|
| 1071 | coef= (number)bla; |
---|
| 1072 | MATELEM(m,i,j)= pOne(); |
---|
| 1073 | pSetCoeff( MATELEM(m,i,j), coef ); |
---|
| 1074 | } |
---|
| 1075 | } |
---|
| 1076 | //PrintLn(); |
---|
| 1077 | } |
---|
| 1078 | |
---|
| 1079 | return m; |
---|
| 1080 | } |
---|
| 1081 | |
---|
| 1082 | intvec * simplex::posvToIV() |
---|
| 1083 | { |
---|
| 1084 | int i; |
---|
| 1085 | intvec * iv = new intvec( m ); |
---|
| 1086 | for ( i= 1; i <= m; i++ ) |
---|
| 1087 | { |
---|
| 1088 | IMATELEM(*iv,i,1)= iposv[i]; |
---|
| 1089 | } |
---|
| 1090 | return iv; |
---|
| 1091 | } |
---|
| 1092 | |
---|
| 1093 | intvec * simplex::zrovToIV() |
---|
| 1094 | { |
---|
| 1095 | int i; |
---|
| 1096 | intvec * iv = new intvec( n ); |
---|
| 1097 | for ( i= 1; i <= n; i++ ) |
---|
| 1098 | { |
---|
| 1099 | IMATELEM(*iv,i,1)= izrov[i]; |
---|
| 1100 | } |
---|
| 1101 | return iv; |
---|
| 1102 | } |
---|
| 1103 | |
---|
| 1104 | void simplex::compute() |
---|
| 1105 | { |
---|
| 1106 | int i,ip,ir,is,k,kh,kp,m12,nl1,nl2; |
---|
| 1107 | int *l1,*l2,*l3; |
---|
| 1108 | mprfloat q1, bmax; |
---|
| 1109 | |
---|
| 1110 | if ( m != (m1+m2+m3) ) |
---|
| 1111 | { |
---|
| 1112 | // error: bad input |
---|
| 1113 | error(WarnS("simplex::compute: Bad input constraint counts!");) |
---|
| 1114 | icase=-2; |
---|
| 1115 | return; |
---|
| 1116 | } |
---|
| 1117 | |
---|
| 1118 | l1= (int *) omAlloc0( (n+1) * sizeof(int) ); |
---|
| 1119 | l2= (int *) omAlloc0( (m+1) * sizeof(int) ); |
---|
| 1120 | l3= (int *) omAlloc0( (m+1) * sizeof(int) ); |
---|
| 1121 | |
---|
| 1122 | nl1= n; |
---|
| 1123 | for ( k=1; k<=n; k++ ) l1[k]=izrov[k]=k; |
---|
| 1124 | nl2=m; |
---|
| 1125 | for ( i=1; i<=m; i++ ) |
---|
| 1126 | { |
---|
| 1127 | if ( LiPM[i+1][1] < 0.0 ) |
---|
| 1128 | { |
---|
| 1129 | // error: bad input |
---|
| 1130 | error(WarnS("simplex::compute: Bad input tableau!");) |
---|
| 1131 | error(Warn("simplex::compute: in input Matrix row %d, column 1, value %f",i+1,LiPM[i+1][1]);) |
---|
| 1132 | icase=-2; |
---|
| 1133 | // free mem l1,l2,l3; |
---|
| 1134 | omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) ); |
---|
| 1135 | omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) ); |
---|
| 1136 | omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) ); |
---|
| 1137 | return; |
---|
| 1138 | } |
---|
| 1139 | l2[i]= i; |
---|
| 1140 | iposv[i]= n+i; |
---|
| 1141 | } |
---|
| 1142 | for ( i=1; i<=m2; i++) l3[i]= 1; |
---|
| 1143 | ir= 0; |
---|
| 1144 | if (m2+m3) |
---|
| 1145 | { |
---|
| 1146 | ir=1; |
---|
| 1147 | for ( k=1; k <= (n+1); k++ ) |
---|
| 1148 | { |
---|
| 1149 | q1=0.0; |
---|
| 1150 | for ( i=m1+1; i <= m; i++ ) q1+= LiPM[i+1][k]; |
---|
| 1151 | LiPM[m+2][k]= -q1; |
---|
| 1152 | } |
---|
| 1153 | |
---|
| 1154 | do |
---|
| 1155 | { |
---|
| 1156 | simp1(LiPM,m+1,l1,nl1,0,&kp,&bmax); |
---|
| 1157 | if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] < -SIMPLEX_EPS ) |
---|
| 1158 | { |
---|
| 1159 | icase= -1; // no solution found |
---|
| 1160 | // free mem l1,l2,l3; |
---|
| 1161 | omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) ); |
---|
| 1162 | omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) ); |
---|
| 1163 | omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) ); |
---|
| 1164 | return; |
---|
| 1165 | } |
---|
| 1166 | else if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] <= SIMPLEX_EPS ) |
---|
| 1167 | { |
---|
| 1168 | m12= m1+m2+1; |
---|
| 1169 | if ( m12 <= m ) |
---|
| 1170 | { |
---|
| 1171 | for ( ip= m12; ip <= m; ip++ ) |
---|
| 1172 | { |
---|
| 1173 | if ( iposv[ip] == (ip+n) ) |
---|
| 1174 | { |
---|
| 1175 | simp1(LiPM,ip,l1,nl1,1,&kp,&bmax); |
---|
| 1176 | if ( fabs(bmax) >= SIMPLEX_EPS) |
---|
| 1177 | goto one; |
---|
| 1178 | } |
---|
| 1179 | } |
---|
| 1180 | } |
---|
| 1181 | ir= 0; |
---|
| 1182 | --m12; |
---|
| 1183 | if ( m1+1 <= m12 ) |
---|
| 1184 | for ( i=m1+1; i <= m12; i++ ) |
---|
| 1185 | if ( l3[i-m1] == 1 ) |
---|
| 1186 | for ( k=1; k <= n+1; k++ ) |
---|
| 1187 | LiPM[i+1][k] = -(LiPM[i+1][k]); |
---|
| 1188 | break; |
---|
| 1189 | } |
---|
| 1190 | //#if DEBUG |
---|
| 1191 | //print_bmat( a, m+2, n+3); |
---|
| 1192 | //#endif |
---|
| 1193 | simp2(LiPM,n,l2,nl2,&ip,kp,&q1); |
---|
| 1194 | if ( ip == 0 ) |
---|
| 1195 | { |
---|
| 1196 | icase = -1; // no solution found |
---|
| 1197 | // free mem l1,l2,l3; |
---|
| 1198 | omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) ); |
---|
| 1199 | omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) ); |
---|
| 1200 | omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) ); |
---|
| 1201 | return; |
---|
| 1202 | } |
---|
| 1203 | one: simp3(LiPM,m+1,n,ip,kp); |
---|
| 1204 | // #if DEBUG |
---|
| 1205 | // print_bmat(a,m+2,n+3); |
---|
| 1206 | // #endif |
---|
| 1207 | if ( iposv[ip] >= (n+m1+m2+1)) |
---|
| 1208 | { |
---|
| 1209 | for ( k= 1; k <= nl1; k++ ) |
---|
| 1210 | if ( l1[k] == kp ) break; |
---|
| 1211 | --nl1; |
---|
| 1212 | for ( is=k; is <= nl1; is++ ) l1[is]= l1[is+1]; |
---|
| 1213 | ++(LiPM[m+2][kp+1]); |
---|
| 1214 | for ( i= 1; i <= m+2; i++ ) LiPM[i][kp+1] = -(LiPM[i][kp+1]); |
---|
| 1215 | } |
---|
| 1216 | else |
---|
| 1217 | { |
---|
| 1218 | if ( iposv[ip] >= (n+m1+1) ) |
---|
| 1219 | { |
---|
| 1220 | kh= iposv[ip]-m1-n; |
---|
| 1221 | if ( l3[kh] ) |
---|
| 1222 | { |
---|
| 1223 | l3[kh]= 0; |
---|
| 1224 | ++(LiPM[m+2][kp+1]); |
---|
| 1225 | for ( i=1; i<= m+2; i++ ) |
---|
| 1226 | LiPM[i][kp+1] = -(LiPM[i][kp+1]); |
---|
| 1227 | } |
---|
| 1228 | } |
---|
| 1229 | } |
---|
| 1230 | is= izrov[kp]; |
---|
| 1231 | izrov[kp]= iposv[ip]; |
---|
| 1232 | iposv[ip]= is; |
---|
| 1233 | } while (ir); |
---|
| 1234 | } |
---|
| 1235 | /* end of phase 1, have feasible sol, now optimize it */ |
---|
| 1236 | loop |
---|
| 1237 | { |
---|
| 1238 | // #if DEBUG |
---|
| 1239 | // print_bmat( a, m+1, n+5); |
---|
| 1240 | // #endif |
---|
| 1241 | simp1(LiPM,0,l1,nl1,0,&kp,&bmax); |
---|
| 1242 | if (bmax <= /*SIMPLEX_EPS*/0.0) |
---|
| 1243 | { |
---|
| 1244 | icase=0; // finite solution found |
---|
| 1245 | // free mem l1,l2,l3 |
---|
| 1246 | omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) ); |
---|
| 1247 | omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) ); |
---|
| 1248 | omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) ); |
---|
| 1249 | return; |
---|
| 1250 | } |
---|
| 1251 | simp2(LiPM,n,l2,nl2,&ip,kp,&q1); |
---|
| 1252 | if (ip == 0) |
---|
| 1253 | { |
---|
| 1254 | //printf("Unbounded:"); |
---|
| 1255 | // #if DEBUG |
---|
| 1256 | // print_bmat( a, m+1, n+1); |
---|
| 1257 | // #endif |
---|
| 1258 | icase=1; /* unbounded */ |
---|
| 1259 | // free mem |
---|
| 1260 | omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) ); |
---|
| 1261 | omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) ); |
---|
| 1262 | omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) ); |
---|
| 1263 | return; |
---|
| 1264 | } |
---|
| 1265 | simp3(LiPM,m,n,ip,kp); |
---|
| 1266 | is= izrov[kp]; |
---|
| 1267 | izrov[kp]= iposv[ip]; |
---|
| 1268 | iposv[ip]= is; |
---|
| 1269 | }/*for ;;*/ |
---|
| 1270 | } |
---|
| 1271 | |
---|
| 1272 | void simplex::simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax ) |
---|
| 1273 | { |
---|
| 1274 | int k; |
---|
| 1275 | mprfloat test; |
---|
| 1276 | |
---|
| 1277 | if( nll <= 0) |
---|
| 1278 | { /* init'tion: fixed */ |
---|
| 1279 | *bmax = 0.0; |
---|
| 1280 | return; |
---|
| 1281 | } |
---|
| 1282 | *kp=ll[1]; |
---|
| 1283 | *bmax=a[mm+1][*kp+1]; |
---|
| 1284 | for (k=2;k<=nll;k++) |
---|
| 1285 | { |
---|
| 1286 | if (iabf == 0) |
---|
| 1287 | { |
---|
| 1288 | test=a[mm+1][ll[k]+1]-(*bmax); |
---|
| 1289 | if (test > 0.0) |
---|
| 1290 | { |
---|
| 1291 | *bmax=a[mm+1][ll[k]+1]; |
---|
| 1292 | *kp=ll[k]; |
---|
| 1293 | } |
---|
| 1294 | } |
---|
| 1295 | else |
---|
| 1296 | { /* abs values: have fixed it */ |
---|
| 1297 | test=fabs(a[mm+1][ll[k]+1])-fabs(*bmax); |
---|
| 1298 | if (test > 0.0) |
---|
| 1299 | { |
---|
| 1300 | *bmax=a[mm+1][ll[k]+1]; |
---|
| 1301 | *kp=ll[k]; |
---|
| 1302 | } |
---|
| 1303 | } |
---|
| 1304 | } |
---|
| 1305 | } |
---|
| 1306 | |
---|
| 1307 | void simplex::simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 ) |
---|
| 1308 | { |
---|
| 1309 | int k,ii,i; |
---|
| 1310 | mprfloat qp,q0,q; |
---|
| 1311 | |
---|
| 1312 | *ip= 0; |
---|
| 1313 | for ( i=1; i <= nl2; i++ ) |
---|
| 1314 | { |
---|
| 1315 | if ( a[l2[i]+1][kp+1] < -SIMPLEX_EPS ) |
---|
| 1316 | { |
---|
| 1317 | *q1= -a[l2[i]+1][1] / a[l2[i]+1][kp+1]; |
---|
| 1318 | *ip= l2[i]; |
---|
| 1319 | for ( i= i+1; i <= nl2; i++ ) |
---|
| 1320 | { |
---|
| 1321 | ii= l2[i]; |
---|
| 1322 | if (a[ii+1][kp+1] < -SIMPLEX_EPS) |
---|
| 1323 | { |
---|
| 1324 | q= -a[ii+1][1] / a[ii+1][kp+1]; |
---|
| 1325 | if (q - *q1 < -SIMPLEX_EPS) |
---|
| 1326 | { |
---|
| 1327 | *ip=ii; |
---|
| 1328 | *q1=q; |
---|
| 1329 | } |
---|
| 1330 | else if (q - *q1 < SIMPLEX_EPS) |
---|
| 1331 | { |
---|
| 1332 | for ( k=1; k<= n; k++ ) |
---|
| 1333 | { |
---|
| 1334 | qp= -a[*ip+1][k+1]/a[*ip+1][kp+1]; |
---|
| 1335 | q0= -a[ii+1][k+1]/a[ii+1][kp+1]; |
---|
| 1336 | if ( q0 != qp ) break; |
---|
| 1337 | } |
---|
| 1338 | if ( q0 < qp ) *ip= ii; |
---|
| 1339 | } |
---|
| 1340 | } |
---|
| 1341 | } |
---|
| 1342 | } |
---|
| 1343 | } |
---|
| 1344 | } |
---|
| 1345 | |
---|
| 1346 | void simplex::simp3( mprfloat **a, int i1, int k1, int ip, int kp ) |
---|
| 1347 | { |
---|
| 1348 | int kk,ii; |
---|
| 1349 | mprfloat piv; |
---|
| 1350 | |
---|
| 1351 | piv= 1.0 / a[ip+1][kp+1]; |
---|
| 1352 | for ( ii=1; ii <= i1+1; ii++ ) |
---|
| 1353 | { |
---|
| 1354 | if ( ii -1 != ip ) |
---|
| 1355 | { |
---|
| 1356 | a[ii][kp+1] *= piv; |
---|
| 1357 | for ( kk=1; kk <= k1+1; kk++ ) |
---|
| 1358 | if ( kk-1 != kp ) |
---|
| 1359 | a[ii][kk] -= a[ip+1][kk] * a[ii][kp+1]; |
---|
| 1360 | } |
---|
| 1361 | } |
---|
| 1362 | for ( kk=1; kk<= k1+1; kk++ ) |
---|
| 1363 | if ( kk-1 != kp ) a[ip+1][kk] *= -piv; |
---|
| 1364 | a[ip+1][kp+1]= piv; |
---|
| 1365 | } |
---|
| 1366 | //<- |
---|
| 1367 | |
---|
| 1368 | //----------------------------------------------------------------------------- |
---|
| 1369 | |
---|
| 1370 | //#endif // HAVE_MPR |
---|
| 1371 | |
---|
| 1372 | // local Variables: *** |
---|
| 1373 | // folded-file: t *** |
---|
| 1374 | // compile-command-1: "make installg" *** |
---|
| 1375 | // compile-command-2: "make install" *** |
---|
| 1376 | // End: *** |
---|