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