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