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