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++ )
158 for (
j= 0;
j <
cn;
j++ )
175 for (
i= 1;
i <
cn;
i++ ) {
180 for (
j= (
cn-
i-1);
j <= (
cn-2);
j++) {
193 for (
i= 0;
i <
cn;
i++ ) {
204 for (
k=
cn-1;
k >= 1;
k-- ) {
260#define MAXIT (5*MMOD)
287 for (
i=0;
i <=
tdg;
i++ )
312 for (
i=0;
i <=
tdg;
i++ )
343 for (
i=
tdg;
i >= 0;
i-- )
390 if (! ((
i >= 0) && (
i <
anz+2) ) )
391 WarnS(
"rootContainer::evPointCoord: index out of range");
393 WarnS(
"rootContainer::evPointCoord: ievpoint == NULL");
405 Warn(
"rootContainer::evPointCoord: NULL index %d",
i);
410 Warn(
"rootContainer::evPointCoord: Wrong index %d, found_roots %s",
i,
found_roots?
"true":
"false");
419 if (
found_roots && ( from >= 0) && ( from <
tdg ) && ( to >= 0) && ( to <
tdg ) )
431 Warn(
" rootContainer::changeRoots: Wrong index %d, %d",from,to);
447 for (
i=0;
i <=
tdg;
i++ )
456 WarnS(
"rootContainer::solver: No roots found!");
495 WarnS(
"Laguerre solver: Too many iterations!");
504 WarnS(
"Laguerre solver: Too many iterations in polish!");
509 if ((!type)&&(!((
x.real()==zero)&&(
x.imag()==zero))))
x = o/
x;
510 if (
x.imag() == zero)
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 };
574 if ((
fabs==zero) || (
abs(d)==zero))
return;
582 sq=
sqrt(( (
h * deg ) -
g2 ) * (deg - one));
591 if((
gp.real()==zero)&&(
gp.imag()==zero))
629 for (
int i=
tdg;
i >= 0;
i-- )
645 for (
i=
j-1;
i > 0;
i-- )
646 *a[
i] += (*a[
i+1]*
x);
647 for (
i= 0;
i <
j;
i++ )
653 for (
i= 1;
i <
j;
i++)
654 *a[
i] += (*a[
i-1]*
y);
662 q((
x.real()*
x.real())+(
x.imag()*
x.imag()));
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++ )
677 for (
i= 2;
i <
j-1;
i++)
678 *a[
i] += ((*a[
i-1]*
p)-(*a[
i-2]*q));
687 &&((!(*a[2]).real().
isZero())||(!(*a[2]).imag().
isZero())))
692 if (
disk.imag().isZero())
694 if (
disk.real()<zero)
707 if(
sq.imag().isZero())
720 if (((*a[1]).real().
isZero()) && ((*a[1]).imag().isZero()))
722 WerrorS(
"precision lost, try again with higher precision");
761 for (
i=
l+inc;
i<=u;
i+=inc)
763 if (r[
i]->real()<
x->real())
773 for (
i=pos;
i>
l;
i--)
780 for (
i=pos+1;
i+1>
l;
i--)
782 if (
x->imag()>
y->imag())
794 else if ((inc==2)&&(
x->imag()<r[
l+1]->
imag()))
813 for (
k=
m-1;
k >= 0;
k-- )
815 f2 = (
x * f2 ) + f1;
816 f1 = (
x * f1 ) +
f0;
834 for (
k= 1;
k <=
m;
k++ )
836 f2 = (
x * f2 ) + f1;
837 f1 = (
x * f1 ) +
f0;
865 for (
i= 0;
i <
rc;
i++ )
873 for (
i= 0;
i <
mc;
i++ )
894 for ( r= 0; r <
anzr; r++ ) {
923 WarnS(
"rootArranger::arrange: precision lost");
930 Warn(
"rootArranger::arrange: No match? coord %d, root %d.",
xkoord,r);
932 WarnS(
"One of these ...");
943 WarnS(
" ... must match to one of these:");
968#define MAXPOINTS 1000
973 : LiPM_cols(cols), LiPM_rows(rows)
1077 for (
i= 1;
i <=
m;
i++ )
1088 for (
i= 1;
i <=
n;
i++ )
1104 error(
WarnS(
"simplex::compute: Bad input constraint counts!");)
1116 for (
i=1;
i<=
m;
i++ )
1118 if (
LiPM[
i+1][1] < 0.0 )
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]);)
1138 for (
k=1;
k <= (
n+1);
k++ )
1176 if (
l3[
i-
m1] == 1 )
1177 for (
k=1;
k <=
n+1;
k++ )
1200 for (
k= 1;
k <=
nl1;
k++ )
1201 if (
l1[
k] ==
kp )
break;
1203 for ( is=
k; is <=
nl1; is++ )
l1[is]=
l1[is+1];
1216 for (
i=1;
i<=
m+2;
i++ )
1304 for (
i=1;
i <=
nl2;
i++ )
1315 q= -a[
ii+1][1] / a[
ii+1][
kp+1];
1323 for (
k=1;
k<=
nn;
k++ )
1327 if (
q0 !=
qp )
break;
1342 piv= 1.0 / a[
ip+1][
kp+1];
1354 if (
kk-1 !=
kp ) a[
ip+1][
kk] *= -piv;
Rational pow(const Rational &a, int e)
Rational abs(const Rational &a)
gmp_complex numbers based on
rootArranger(rootContainer **_roots, rootContainer **_mu, const int _howclean=PM_CORRUPT)
complex root finder for univariate polynomials based on laguers algorithm
void sortre(gmp_complex **r, int l, int u, int inc)
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...
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)
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
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 divlin(gmp_complex **a, gmp_complex x, int j)
void sortroots(gmp_complex **roots, int r, int c, bool isf)
void divquad(gmp_complex **a, gmp_complex x, int j)
bool swapRoots(const int from, const int to)
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)
void solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j)
gmp_complex & evPointCoord(const int i)
bool isfloat(gmp_complex **a)
void checkimag(gmp_complex *x, gmp_float &e)
bool solver(const int polishmode=PM_NONE)
BOOLEAN mapFromMatrix(matrix m)
simplex(int rows, int cols)
#rows should be >= m+2, #cols >= n+1
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)
poly numvec2poly(const number *q)
vandermonde(const long _cn, const long _n, const long _maxdeg, number *_p, const bool _homog=true)
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.
const CanonicalForm int s
const CanonicalForm int const CFList const Variable & y
bool isZero(const CFArray &A)
checks if entries of A are zero
void WerrorS(const char *s)
#define IMATELEM(M, I, J)
static matrix mu(matrix A, const ring R)
#define MATELEM(mat, i, j)
1-based access to matrix
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
EXTERN_VAR size_t gmp_output_digits
gmp_float sin(const gmp_float &a)
char * complexToStr(gmp_complex &c, const unsigned int oprec, const coeffs src)
gmp_float sqrt(const gmp_float &a)
gmp_float exp(const gmp_float &a)
gmp_float cos(const gmp_float &a)
gmp_complex numberToComplex(number num, const coeffs r)
#define mprSTICKYPROT(msg)
The main handler for Singular numbers which are suitable for Singular polynomials.
#define nPower(a, b, res)
#define omFreeSize(addr, size)
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Compatibility layer for legacy polynomial operations (over currRing)
#define pSortAdd(p)
sorts p, p may have equal monomials
#define pSetCoeff(p, n)
deletes old coeff before setting the new one