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