36 number *_p,
const bool _homog )
37 : n(_n), cn(_cn), maxdeg(_maxdeg),
p(_p), homog(_homog)
42 for ( j= 0; j <
cn; j++ )
x[j]=
nInit(1);
63 for ( j= 0; j <
n; j++ ) exp[j]=0;
65 for ( i= 0; i <
l; i++ )
69 for ( j= 0; j <
n; j++ )
80 for ( j= 0; j < n - 1; j++ )
105 for ( j= 0; j <
n+1; j++ ) exp[j]=0;
107 for ( i= 0; i <
l; i++ )
128 for ( j= 1; j <
n; j++ )
140 omFreeSize( (
void *) exp, (n+1) *
sizeof(
int) );
156 w= (number *)
omAlloc(
cn *
sizeof(number) );
157 c= (number *)
omAlloc(
cn *
sizeof(number) );
158 for ( j= 0; j <
cn; j++ )
175 for ( i= 1; i <
cn; i++ ) {
180 for ( j= (cn-i-1); j <= (cn-2); j++) {
182 tmp1=
nMult( xx, c[j+1] );
183 newnum=
nAdd( c[j], tmp1 );
188 newnum=
nAdd( xx, c[cn-1] );
193 for ( i= 0; i <
cn; i++ ) {
204 for ( k= cn-1; k >= 1; k-- ) {
206 tmp1=
nMult( xx, b );
208 b=
nAdd( c[k], tmp1 );
211 tmp1=
nMult( q[k-1], b );
212 newnum=
nAdd( s, tmp1 );
217 tmp1=
nMult( xx, t );
218 newnum=
nAdd( tmp1, b );
236 for ( j= 0; j <
cn; j++ )
nDelete( c+j );
237 omFreeSize( (
void *)c, cn *
sizeof( number ) );
257 #define MR 8 // never change this value 260 #define MAXIT (5*MMOD) // max number of iterations in laguer root finder 281 if ( ievpoint !=
NULL )
283 for ( i=0; i < anz+2; i++ )
nDelete( ievpoint + i );
284 omFreeSize( (
void *)ievpoint, (anz+2) *
sizeof( number ) );
291 for ( i=0; i < tdg; i++ )
delete theroots[i];
300 const int _var,
const int _tdg,
301 const rootType _rt,
const int _anz )
311 for ( i=0; i <= tdg; i++ )
321 if ( rt == cspecialmu && _ievpoint )
323 ievpoint= (number *)
omAlloc( (anz+2) *
sizeof( number ) );
324 for (i=0; i < anz+2; i++) ievpoint[i]=
nCopy( _ievpoint[i] );
340 if ( (rt == cspecial) || ( rt == cspecialmu ) )
342 for ( i= tdg; i >= 0; i-- )
389 if (! ((i >= 0) && (i < anz+2) ) )
390 WarnS(
"rootContainer::evPointCoord: index out of range");
391 if (ievpoint ==
NULL)
392 WarnS(
"rootContainer::evPointCoord: ievpoint == NULL");
394 if ( (rt == cspecialmu) && found_roots )
396 if ( ievpoint[i] !=
NULL )
404 Warn(
"rootContainer::evPointCoord: NULL index %d",i);
409 Warn(
"rootContainer::evPointCoord: Wrong index %d, found_roots %s",i,found_roots?
"true":
"false");
418 if ( found_roots && ( from >= 0) && ( from < tdg ) && ( to >= 0) && ( to < tdg ) )
423 *theroots[from]= *theroots[to];
430 Warn(
" rootContainer::changeRoots: Wrong index %d, %d",from,to);
442 for ( i=0; i < tdg; i++ ) theroots[i]=
new gmp_complex();
446 for ( i=0; i <= tdg; i++ )
453 found_roots= laguer_driver( ad, theroots, polishmode != 0 );
455 WarnS(
"rootContainer::solver: No roots found!");
458 for ( i=0; i <= tdg; i++ )
delete ad[i];
471 bool ret=
true, isf=isfloat(a), type=
true;
474 for ( i=0; i <= tdg; i++ ) ad[i]=
new gmp_complex( *a[i] );
483 laguer(ad, i, &
x, &its, type);
488 laguer(ad, i, &
x, &its, type);
494 WarnS(
"Laguerre solver: Too many iterations!");
500 laguer( a, tdg, &
x, &its , type);
503 WarnS(
"Laguerre solver: Too many iterations in polish!");
508 if ((!type)&&(!((
x.real()==zero)&&(
x.imag()==zero))))
x = o/
x;
509 if (
x.imag() == zero)
536 solvequad(ad,roots,k,j);
537 sortroots(roots,k,j,isf);
541 for ( i=0; i <= tdg; i++ )
delete ad[i];
554 gmp_complex dx, x1,
b, d,
f,
g,
h, sq,
gp, gm, g2;
555 gmp_float frac_g[
MR+1] = { 0.0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0 };
560 for ( iter= 1; iter <=
MAXIT; iter++ )
565 computefx(a,*x,m,b,d,f,abx_g,err_g);
567 computegx(a,*x,m,b,d,f,abx_g,err_g);
573 if ((fabs==zero) || (
abs(d)==zero))
return;
580 h= g2 - (((f+
f) / b ));
581 sq=
sqrt(( ( h * deg ) - g2 ) * (deg - one));
590 if((gp.
real()==zero)&&(gp.
imag()==zero))
603 if (*x == x1)
goto ende;
607 if ( j %
MT ) *x= x1;
608 else *x -= ( dx * frac_g[ j /
MT ] );
628 for (
int i=tdg;
i >= 0;
i-- )
644 for (i= j-1; i > 0; i-- )
645 *a[i] += (*a[i+1]*x);
646 for (i= 0; i <
j; i++ )
652 for (i= 1; i <
j; i++)
653 *a[i] += (*a[i-1]*y);
665 *a[j-1] += (*a[
j]*
p);
666 for (i= j-2; i > 1; i-- )
667 *a[i] += ((*a[i+1]*p)-(*a[i+2]*q));
668 for (i= 0; i < j-1; i++ )
676 for (i= 2; i < j-1; i++)
677 *a[i] += ((*a[i-1]*p)-(*a[i-2]*q));
686 &&((!(*a[2]).real().
isZero())||(!(*a[2]).imag().
isZero())))
689 gmp_complex h1(*a[1]/(*a[2] + *a[2])), h2(*a[0] / *a[2]);
693 if (disk.
real()<zero)
719 if (((*a[1]).real().
isZero()) && ((*a[1]).imag().isZero()))
721 WerrorS(
"precision lost, try again with higher precision");
743 for (j=c; j+2<tdg; j+=2)
744 sortre(ro, j, tdg-1, 2);
748 for (j=c; j+1<tdg; j++)
749 sortre(ro, j, tdg-1, 1);
760 for (i=l+inc; i<=u; i+=inc)
762 if (r[i]->real()<x->
real())
772 for (i=pos; i>
l; i--)
779 for (i=pos+1; i+1>
l; i--)
793 else if ((inc==2)&&(x->
imag()<r[l+1]->
imag()))
812 for ( k= m-1; k >= 0; k-- )
814 f2 = ( x * f2 ) + f1;
815 f1 = ( x * f1 ) + f0;
816 f0 = ( x * f0 ) + *a[k];
817 ef =
abs( f0 ) + ( ex * ef );
833 for ( k= 1; k <=
m; k++ )
835 f2 = ( x * f2 ) + f1;
836 f1 = ( x * f1 ) + f0;
837 f0 = ( x * f0 ) + *a[k];
838 ef =
abs( f0 ) + ( ex * ef );
849 const int _howclean )
850 : roots(_roots),
mu(_mu), howclean(_howclean)
864 for ( i= 0; i <
rc; i++ )
872 for ( i= 0; i <
mc; i++ )
887 int xkoord, r, rtest, xk, mtest;
891 for ( xkoord= 0; xkoord < anzm; xkoord++ ) {
893 for ( r= 0; r < anzr; r++ ) {
897 for ( xk =0; xk <= xkoord; xk++ )
899 tmp -= (*
roots[xk])[r] *
mu[xkoord]->evPointCoord(xk+1);
903 for ( rtest= r; rtest < anzr; rtest++ ) {
904 zwerg = tmp - (*
roots[xk])[rtest] *
mu[xkoord]->evPointCoord(xk+1);
905 for ( mtest= 0; mtest < anzr; mtest++ )
909 if ( ((zwerg.
real() <= (*
mu[xkoord])[mtest].real() + mprec) &&
910 (zwerg.
real() >= (*
mu[xkoord])[mtest].real() - mprec)) &&
911 ((zwerg.
imag() <= (*
mu[xkoord])[mtest].imag() + mprec) &&
912 (zwerg.
imag() >= (*
mu[xkoord])[mtest].imag() - mprec)) )
922 WarnS(
"rootArranger::arrange: precision lost");
929 Warn(
"rootArranger::arrange: No match? coord %d, root %d.",xkoord,r);
931 WarnS(
"One of these ...");
932 for ( rtest= r; rtest < anzr; rtest++ )
935 for ( xk =0; xk <= xkoord; xk++ )
937 tmp-= (*
roots[xk])[r] *
mu[xkoord]->evPointCoord(xk+1);
939 tmp-= (*
roots[xk])[rtest] *
mu[xkoord]->evPointCoord(xk+1);
942 WarnS(
" ... must match to one of these:");
943 for ( mtest= 0; mtest < anzr; mtest++ )
967 #define MAXPOINTS 1000 972 : LiPM_cols(cols), LiPM_rows(rows)
1019 for ( i= 1; i <=
MATROWS( mm ); i++ )
1021 for ( j= 1; j <=
MATCOLS( mm ); j++ )
1051 for ( i= 1; i <=
MATROWS( mm ); i++ )
1053 for ( j= 1; j <=
MATCOLS( mm ); j++ )
1058 if (
LiPM[i][j] != 0.0 )
1076 for ( i= 1; i <=
m; i++ )
1087 for ( i= 1; i <=
n; i++ )
1096 int i,ip,ir,is,
k,kh,kp,m12,nl1,nl2;
1103 error(
WarnS(
"simplex::compute: Bad input constraint counts!");)
1108 l1= (
int *)
omAlloc0( (
n+1) *
sizeof(int) );
1109 l2= (
int *)
omAlloc0( (
m+1) *
sizeof(int) );
1110 l3= (
int *)
omAlloc0( (
m+1) *
sizeof(int) );
1113 for ( k=1; k<=
n; k++ ) l1[k]=
izrov[k]=k;
1115 for ( i=1; i<=
m; i++ )
1117 if (
LiPM[i+1][1] < 0.0 )
1120 error(
WarnS(
"simplex::compute: Bad input tableau!");)
1121 error(
Warn(
"simplex::compute: in input Matrix row %d, column 1, value %f",i+1,
LiPM[i+1][1]);)
1132 for ( i=1; i<=
m2; i++) l3[i]= 1;
1137 for ( k=1; k <= (
n+1); k++ )
1140 for ( i=
m1+1; i <=
m; i++ ) q1+=
LiPM[i+1][k];
1161 for ( ip= m12; ip <=
m; ip++ )
1163 if (
iposv[ip] == (ip+
n) )
1174 for ( i=
m1+1; i <= m12; i++ )
1175 if ( l3[i-
m1] == 1 )
1176 for ( k=1; k <=
n+1; k++ )
1199 for ( k= 1; k <= nl1; k++ )
1200 if ( l1[k] == kp )
break;
1202 for ( is=k; is <= nl1; is++ ) l1[is]= l1[is+1];
1203 ++(
LiPM[
m+2][kp+1]);
1204 for ( i= 1; i <=
m+2; i++ )
LiPM[i][kp+1] = -(
LiPM[i][kp+1]);
1214 ++(
LiPM[
m+2][kp+1]);
1215 for ( i=1; i<=
m+2; i++ )
1273 *bmax=a[mm+1][*kp+1];
1274 for (k=2;k<=nll;k++)
1278 test=a[mm+1][ll[
k]+1]-(*bmax);
1281 *bmax=a[mm+1][ll[
k]+1];
1287 test=fabs(a[mm+1][ll[k]+1])-fabs(*bmax);
1290 *bmax=a[mm+1][ll[
k]+1];
1303 for ( i=1; i <= nl2; i++ )
1307 *q1= -a[l2[
i]+1][1] / a[l2[
i]+1][kp+1];
1309 for ( i= i+1; i <= nl2; i++ )
1314 q= -a[ii+1][1] / a[ii+1][kp+1];
1322 for ( k=1; k<= nn; k++ )
1324 qp= -a[*ip+1][k+1]/a[*ip+1][kp+1];
1325 q0= -a[ii+1][k+1]/a[ii+1][kp+1];
1326 if ( q0 != qp )
break;
1328 if ( q0 < qp ) *ip= ii;
1341 piv= 1.0 / a[ip+1][kp+1];
1342 for ( ii=1; ii <= i1+1; ii++ )
1347 for ( kk=1; kk <= k1+1; kk++ )
1349 a[ii][kk] -= a[ip+1][kk] * a[ii][kp+1];
1352 for ( kk=1; kk<= k1+1; kk++ )
1353 if ( kk-1 != kp ) a[ip+1][kk] *= -piv;
#define mprSTICKYPROT(msg)
complex root finder for univariate polynomials based on laguers algorithm
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)
const CanonicalForm int s
const CanonicalForm int const CFList const Variable & y
number * interpolateDense(const number *q)
Solves the Vandermode linear system {i=1}^{n} x_i^k-1 w_i = q_k, k=1,..,n.
matrix mapToMatrix(matrix m)
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)
bool swapRoots(const int from, const int to)
void sortre(gmp_complex **r, int l, int u, int inc)
void mu(int **points, int sizePoints)
gmp_float exp(const gmp_float &a)
Compatiblity layer for legacy polynomial operations (over currRing)
#define nPower(a, b, res)
#define omFreeSize(addr, size)
void checkimag(gmp_complex *x, gmp_float &e)
bool isfloat(gmp_complex **a)
bool solver(const int polishmode=PM_NONE)
void WerrorS(const char *s)
gmp_complex numbers based on
simplex(int rows, int cols)
#rows should be >= m+2, #cols >= n+1
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
poly numvec2poly(const number *q)
gmp_complex & evPointCoord(const int i)
Rational abs(const Rational &a)
BOOLEAN mapFromMatrix(matrix m)
gmp_complex numberToComplex(number num, const coeffs r)
The main handler for Singular numbers which are suitable for Singular polynomials.
void divquad(gmp_complex **a, gmp_complex x, int j)
void sortroots(gmp_complex **roots, int r, int c, bool isf)
gmp_float sqrt(const gmp_float &a)
const mpf_t * mpfp() const
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 ...
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
void simp2(mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1)
gmp_float cos(const gmp_float &a)
#define pSortAdd(p)
sorts p, p may have equal monomials
rootArranger(rootContainer **_roots, rootContainer **_mu, const int _howclean=PM_CORRUPT)
void solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j)
bool isZero(const CFArray &A)
checks if entries of A are zero
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
char * complexToStr(gmp_complex &c, const unsigned int oprec, const coeffs src)
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...
void simp1(mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax)
Rational pow(const Rational &a, int e)
#define IMATELEM(M, I, J)
gmp_float sin(const gmp_float &a)
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
void simp3(mprfloat **a, int i1, int k1, int ip, int kp)
void divlin(gmp_complex **a, gmp_complex x, int j)
#define MATELEM(mat, i, j)
1-based access to matrix
vandermonde(const long _cn, const long _n, const long _maxdeg, number *_p, const bool _homog=true)