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