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