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