49 number *_p,
const bool _homog )
50 : n(_n), cn(_cn), maxdeg(_maxdeg),
p(_p), homog(_homog)
78 for (
i= 0;
i <
l;
i++ )
82 for (
j= 0;
j <
n;
j++ )
93 for (
j= 0;
j <
n - 1;
j++ )
118 for (
j= 0;
j <
n+1;
j++ )
exp[
j]=0;
120 for (
i= 0;
i <
l;
i++ )
141 for (
j= 1;
j <
n;
j++ )
170 c= (number *)
omAlloc(
cn *
sizeof(number) );
171 for (
j= 0;
j <
cn;
j++ )
188 for (
i= 1;
i <
cn;
i++ ) {
193 for (
j= (
cn-
i-1);
j <= (
cn-2);
j++) {
201 newnum=
nAdd( xx, c[
cn-1] );
206 for (
i= 0;
i <
cn;
i++ ) {
217 for (
k=
cn-1;
k >= 1;
k-- ) {
270 #define MR 8 // never change this value 273 #define MAXIT (5*MMOD) // max number of iterations in laguer root finder 313 const int _var,
const int _tdg,
314 const rootType _rt,
const int _anz )
324 for (
i=0;
i <=
tdg;
i++ )
355 for (
i=
tdg;
i >= 0;
i-- )
402 if (! ((
i >= 0) && (
i <
anz+2) ) )
403 WarnS(
"rootContainer::evPointCoord: index out of range");
405 WarnS(
"rootContainer::evPointCoord: ievpoint == NULL");
417 Warn(
"rootContainer::evPointCoord: NULL index %d",
i);
422 Warn(
"rootContainer::evPointCoord: Wrong index %d, found_roots %s",
i,
found_roots?
"true":
"false");
431 if (
found_roots && ( from >= 0) && ( from <
tdg ) && ( to >= 0) && ( to <
tdg ) )
443 Warn(
" rootContainer::changeRoots: Wrong index %d, %d",from,to);
459 for (
i=0;
i <=
tdg;
i++ )
468 WarnS(
"rootContainer::solver: No roots found!");
471 for (
i=0;
i <=
tdg;
i++ )
delete ad[
i];
484 bool ret=
true, isf=
isfloat(
a), type=
true;
507 WarnS(
"Laguerre solver: Too many iterations!");
516 WarnS(
"Laguerre solver: Too many iterations in polish!");
521 if ((!type)&&(!((
x.real()==zero)&&(
x.imag()==zero))))
x = o/
x;
522 if (
x.imag() == zero)
554 for (
i=0;
i <=
tdg;
i++ )
delete ad[
i];
567 gmp_complex dx, x1,
b, d,
f,
g,
h, sq,
gp, gm, g2;
568 gmp_float frac_g[
MR+1] = { 0.0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0 };
586 if ((fabs==zero) || (
abs(d)==zero))
return;
593 h= g2 - (((
f+
f) /
b ));
594 sq=
sqrt(( (
h * deg ) - g2 ) * (deg - one));
603 if((
gp.real()==zero)&&(
gp.imag()==zero))
616 if (*
x == x1)
goto ende;
620 if (
j %
MT ) *
x= x1;
621 else *
x -= ( dx * frac_g[
j /
MT ] );
641 for (
int i=
tdg;
i >= 0;
i-- )
657 for (
i=
j-1;
i > 0;
i-- )
658 *
a[
i] += (*
a[
i+1]*
x);
659 for (
i= 0;
i <
j;
i++ )
665 for (
i= 1;
i <
j;
i++)
666 *
a[
i] += (*
a[
i-1]*
y);
674 q((
x.real()*
x.real())+(
x.imag()*
x.imag()));
678 *
a[
j-1] += (*
a[
j]*
p);
679 for (
i=
j-2;
i > 1;
i-- )
680 *
a[
i] += ((*
a[
i+1]*
p)-(*
a[
i+2]*q));
681 for (
i= 0;
i <
j-1;
i++ )
689 for (
i= 2;
i <
j-1;
i++)
690 *
a[
i] += ((*
a[
i-1]*
p)-(*
a[
i-2]*q));
706 if (disk.
real()<zero)
732 if (((*
a[1]).real().
isZero()) && ((*
a[1]).imag().isZero()))
734 WerrorS(
"precision lost, try again with higher precision");
739 if(
r[
k]->imag().isZero())
773 for (
i=
l+inc;
i<=u;
i+=inc)
775 if (
r[
i]->real()<
x->real())
785 for (
i=pos;
i>
l;
i--)
792 for (
i=pos+1;
i+1>
l;
i--)
794 if (
x->imag()>
y->imag())
806 else if ((inc==2)&&(
x->imag()<
r[
l+1]->imag()))
825 for (
k=
m-1;
k >= 0;
k-- )
827 f2 = (
x * f2 ) + f1;
828 f1 = (
x * f1 ) + f0;
829 f0 = (
x * f0 ) + *
a[
k];
830 ef =
abs( f0 ) + ( ex * ef );
846 for (
k= 1;
k <=
m;
k++ )
848 f2 = (
x * f2 ) + f1;
849 f1 = (
x * f1 ) + f0;
850 f0 = (
x * f0 ) + *
a[
k];
851 ef =
abs( f0 ) + ( ex * ef );
862 const int _howclean )
863 : roots(_roots),
mu(_mu), howclean(_howclean)
877 for (
i= 0;
i <
rc;
i++ )
885 for (
i= 0;
i <
mc;
i++ )
900 int xkoord,
r, rtest, xk, mtest;
904 for ( xkoord= 0; xkoord < anzm; xkoord++ ) {
906 for (
r= 0;
r < anzr;
r++ ) {
910 for ( xk =0; xk <= xkoord; xk++ )
912 tmp -= (*
roots[xk])[
r] *
mu[xkoord]->evPointCoord(xk+1);
916 for ( rtest=
r; rtest < anzr; rtest++ ) {
917 zwerg = tmp - (*
roots[xk])[rtest] *
mu[xkoord]->evPointCoord(xk+1);
918 for ( mtest= 0; mtest < anzr; mtest++ )
922 if ( ((zwerg.
real() <= (*
mu[xkoord])[mtest].real() + mprec) &&
923 (zwerg.
real() >= (*
mu[xkoord])[mtest].real() - mprec)) &&
924 ((zwerg.
imag() <= (*
mu[xkoord])[mtest].imag() + mprec) &&
925 (zwerg.
imag() >= (*
mu[xkoord])[mtest].imag() - mprec)) )
935 WarnS(
"rootArranger::arrange: precision lost");
942 Warn(
"rootArranger::arrange: No match? coord %d, root %d.",xkoord,
r);
944 WarnS(
"One of these ...");
945 for ( rtest=
r; rtest < anzr; rtest++ )
948 for ( xk =0; xk <= xkoord; xk++ )
950 tmp-= (*
roots[xk])[
r] *
mu[xkoord]->evPointCoord(xk+1);
952 tmp-= (*
roots[xk])[rtest] *
mu[xkoord]->evPointCoord(xk+1);
955 WarnS(
" ... must match to one of these:");
956 for ( mtest= 0; mtest < anzr; mtest++ )
980 #define MAXPOINTS 1000 985 : LiPM_cols(cols), LiPM_rows(rows)
1089 for (
i= 1;
i <=
m;
i++ )
1100 for (
i= 1;
i <=
n;
i++ )
1109 int i,ip,ir,is,
k,kh,kp,m12,nl1,nl2;
1116 error(
WarnS(
"simplex::compute: Bad input constraint counts!");)
1121 l1= (
int *)
omAlloc0( (
n+1) *
sizeof(int) );
1122 l2= (
int *)
omAlloc0( (
m+1) *
sizeof(int) );
1123 l3= (
int *)
omAlloc0( (
m+1) *
sizeof(int) );
1128 for (
i=1;
i<=
m;
i++ )
1130 if (
LiPM[
i+1][1] < 0.0 )
1133 error(
WarnS(
"simplex::compute: Bad input tableau!");)
1134 error(
Warn(
"simplex::compute: in input Matrix row %d, column 1, value %f",
i+1,
LiPM[
i+1][1]);)
1145 for (
i=1;
i<=
m2;
i++) l3[
i]= 1;
1150 for (
k=1;
k <= (
n+1);
k++ )
1174 for ( ip= m12; ip <=
m; ip++ )
1176 if (
iposv[ip] == (ip+
n) )
1187 for (
i=
m1+1;
i <= m12;
i++ )
1188 if ( l3[
i-
m1] == 1 )
1189 for (
k=1;
k <=
n+1;
k++ )
1212 for (
k= 1;
k <= nl1;
k++ )
1213 if ( l1[
k] == kp )
break;
1215 for ( is=
k; is <= nl1; is++ ) l1[is]= l1[is+1];
1216 ++(
LiPM[
m+2][kp+1]);
1227 ++(
LiPM[
m+2][kp+1]);
1228 for (
i=1;
i<=
m+2;
i++ )
1286 *bmax=
a[mm+1][*kp+1];
1287 for (
k=2;
k<=nll;
k++)
1291 test=
a[mm+1][ll[
k]+1]-(*bmax);
1294 *bmax=
a[mm+1][ll[
k]+1];
1300 test=fabs(
a[mm+1][ll[
k]+1])-fabs(*bmax);
1303 *bmax=
a[mm+1][ll[
k]+1];
1316 for (
i=1;
i <= nl2;
i++ )
1320 *q1= -
a[l2[
i]+1][1] /
a[l2[
i]+1][kp+1];
1322 for (
i=
i+1;
i <= nl2;
i++ )
1327 q= -
a[ii+1][1] /
a[ii+1][kp+1];
1335 for (
k=1;
k<= nn;
k++ )
1337 qp= -
a[*ip+1][
k+1]/
a[*ip+1][kp+1];
1338 q0= -
a[ii+1][
k+1]/
a[ii+1][kp+1];
1339 if ( q0 != qp )
break;
1341 if ( q0 < qp ) *ip= ii;
1354 piv= 1.0 /
a[ip+1][kp+1];
1355 for ( ii=1; ii <= i1+1; ii++ )
1360 for ( kk=1; kk <= k1+1; kk++ )
1362 a[ii][kk] -=
a[ip+1][kk] *
a[ii][kp+1];
1365 for ( kk=1; kk<= k1+1; kk++ )
1366 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)
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)
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
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
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)
vandermonde(const long _cn, const long _n, const long _maxdeg, number *_p, const bool _homog=true)