My Project
Loading...
Searching...
No Matches
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::*
35vandermonde::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
93poly 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
146number * 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 loeschen
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( ... )
300void 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 )
417bool 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()
437bool rootContainer::solver( const int polishmode )
438{
439 int i;
440
441 // there are maximal tdg roots, so *roots ranges form 0 to tdg-1.
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 )
467bool 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
540theend:
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(...)
550void 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;
613ende:
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
735void 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
754void 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++ ) { // fuer 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++ ) { // fuer 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++ ) { // fuer 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//
972simplex::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
1263void 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
1298void 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
1337void 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
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
mpf_t * _mpfp()
Definition: mpr_complex.h:134
bool isZero() const
Definition: mpr_complex.cc:252
const mpf_t * mpfp() const
Definition: mpr_complex.h:133
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 ad[0..tdg] (generated from the number coeffic...
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
static matrix mu(matrix A, const ring R)
Definition: matpol.cc:2025
#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
char * complexToStr(gmp_complex &c, const unsigned int oprec, const coeffs src)
Definition: mpr_complex.cc:704
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
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
Compatibility 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:75