46 #define SM_MIN_LENGTH_BUCKET MIN_LENGTH_BUCKET - 5 48 #define SM_MIN_LENGTH_BUCKET MAX_INT_VAL 51 typedef struct smprec sm_prec;
92 p = currRing->p_Procs->pp_Mult_Coeff_mm_DivSelectMult(p, m, a, b,
154 void smSparseHomog();
172 int smCheckNormalize();
182 void smNewBareiss(
int,
int);
183 void smToIntvec(
intvec *);
202 al = di*
sizeof(long);
204 bl = ra*
sizeof(long);
206 for (i=di-1;i>=0;i--)
214 for (j=
rVar(currRing);j>0;j--)
258 for (i=pos;i<d;i++) c[i] = c[i+1];
267 int *ord=(
int*)
omAlloc0(3*
sizeof(
int));
268 int *block0=(
int*)
omAlloc(3*
sizeof(
int));
269 int *block1=(
int*)
omAlloc(3*
sizeof(
int));
278 tmpR->bitmask = 2*
bound;
279 tmpR->wvhdl = (
int **)
omAlloc0((3) *
sizeof(
int*));
285 if (origR->qideal!=
NULL)
290 Print(
"[%ld:%d]", (
long) tmpR->bitmask, tmpR->ExpL_Size);
361 if (I->ncols != I->rank)
363 Werror(
"det of %ld x %d module (matrix)",I->rank,I->ncols);
464 wrw = (
float *)
omAlloc(
sizeof(
float)*
i);
466 wcl = (
float *)
omAlloc(
sizeof(
float)*
i);
488 if (m_act ==
NULL)
return;
525 for (i=v->
rows()-1; i>=0; i--)
543 if (
act != 0) res = m_act[1]->m;
567 if (
act != 0) res = m_act[1]->m;
594 if (
act != 0) res = m_act[1]->m;
608 if ((x > 0) && (x <
nrows))
674 for (i=tored;
i; i--) wrw[i] = 0.0;
717 if ((wr<0.25) || (wc<0.25))
728 wp = w*(wpoints-wcl[
i]-wr);
747 m_act[
act] = m_act[cpiv];
758 float wc, wp,
w, hp = piv->f;
763 for (i=tored;
i; i--) wrw[i] = 0.0;
777 if (f) w /= m_res[
f]->f;
796 float wopt = 1.0e30, hp = piv->f;
799 int i, copt, ropt,
f,
e = crd;
801 this->smNewWeights();
814 if (f) w /= m_res[
f]->f;
818 if ((wr<0.25) || (wc<0.25))
829 wp = w*(wpoints-wcl[
i]-wr);
848 m_act[
act] = m_act[cpiv];
895 else if (a->pos > b->pos)
929 m_act[r->pos] = dumm->n;
937 poly hp = this->smMultPoly(piv);
968 x =
SM_MULT(b->m, hr, m_res[ir]->m,_R);
970 if(ir)
SM_DIV(x, m_res[ir]->
m,_R);
982 else if (a->pos > b->pos)
985 x =
SM_MULT(b->m, hr, m_res[ir]->m,_R);
987 if(ir)
SM_DIV(x, m_res[ir]->
m,_R);
1000 x =
SM_MULT(ha, m_res[ir]->
m, m_res[ia]->
m,_R);
1003 if (ia)
SM_DIV(ha, m_res[ia]->
m,_R);
1006 x =
SM_MULT(ha, gp, m_res[ia]->
m,_R);
1008 y =
SM_MULT(b->m, hr, m_res[ia]->m,_R);
1014 x =
SM_MULT(ha, m_res[crd]->
m, m_res[ia]->
m,_R);
1022 y =
SM_MULT(y, m_res[ir]->
m, m_res[ip]->
m,_R);
1023 if (ip)
SM_DIV(y, m_res[ip]->
m,_R);
1026 x =
SM_MULT(ha, y, m_res[ia]->
m,_R);
1029 y =
SM_MULT(b->m, hr, m_res[ia]->m,_R);
1033 x =
SM_MULT(hr, m_res[ia]->
m, m_res[ir]->
m,_R);
1034 if (ir)
SM_DIV(x, m_res[ir]->
m,_R);
1035 y =
SM_MULT(b->m, x, m_res[ia]->m,_R);
1037 x =
SM_MULT(ha, gp, m_res[ia]->
m,_R);
1043 if (ia)
SM_DIV(ha, m_res[ia]->
m,_R);
1063 m_act[r->pos] = dumm->n;
1065 }
while (r !=
NULL);
1094 }
while (a->pos < rpiv);
1101 for (i=1; i<
act; i++)
1110 if ((a ==
NULL) || (a->pos > rpiv))
1115 a->m =
p_Neg(a->m,_R);
1122 else if (a->pos == rpiv)
1125 a->m =
p_Neg(a->m,_R);
1147 h->n = m_row[h->pos];
1199 if (i >
act)
return;
1200 if (m_act[i] ==
NULL)
break;
1207 if (m_act[j] !=
NULL)
1209 m_act[
i] = m_act[
j];
1229 if (i >
act)
return;
1230 if (m_act[i]->
pos > tored)
1232 m_res[inred] = m_act[
i];
1242 if (m_act[j]->
pos > tored)
1244 m_res[inred] = m_act[
j];
1249 m_act[
i] = m_act[
j];
1272 perm[crd+
i] = a->pos;
1274 }
while ((a !=
NULL) && (a->pos <= tored));
1281 if (
perm[crd+k] >= a->pos)
1283 if (
perm[crd+k] > a->pos)
1285 for (l=i;l>=
k;l--)
perm[crd+l+1] =
perm[crd+l];
1286 perm[crd+
k] = a->pos;
1290 if ((a ==
NULL) || (a->pos > tored))
break;
1293 if ((k > i) && (a->pos <= tored))
1298 perm[crd+
i] = a->pos;
1300 }
while ((a !=
NULL) && (a->pos <= tored));
1310 while ((a !=
NULL) && (a->pos <= tored))
1312 if (
perm[crd+k] == a->pos)
1324 if (m_row[j] !=
NULL)
1350 m_res[crd] = m_act[
act];
1353 for (i=1;i<=tored;i++)
1355 if(m_row[i] !=
NULL)
1380 for (i=tored+1;i<=
nrows;i++)
1382 if(m_row[i] !=
NULL)
1405 while (inred <
ncols)
1409 m_res[crd] = m_res[inred];
1428 ha =
SM_MULT(a->m, m_res[e]->m, m_res[f]->m,_R);
1430 if (f)
SM_DIV(ha, m_res[f]->
m,_R);
1456 ha =
SM_MULT(a->m, m_res[e]->m, m_res[f]->m, _R);
1458 if (f)
SM_DIV(ha, m_res[f]->
m, _R);
1463 }
while (a !=
NULL);
1482 }
while (a !=
NULL);
1503 }
while (a !=
NULL);
1518 h =
SM_MULT(h, m_res[crd]->
m, m_res[f]->
m, _R);
1519 if (f)
SM_DIV(h, m_res[f]->
m, _R);
1543 }
while (a !=
NULL);
1629 for (i=
rVar(R);
i; i--)
1637 }
while (a !=
NULL);
1643 for (i=
rVar(R);
i; i--)
1656 for (i=
rVar(R);
i; i--)
1661 }
while (t !=
NULL);
1676 poly smMultDiv(poly a, poly b, const poly c)
1682 if ((c == NULL) || pLmIsConstantComp(c))
1684 return pp_Mult_qq(a, b);
1687 pqLength(a, b, la, lb, SM_MIN_LENGTH_BUCKET);
1689 // we iter over b, so, make sure b is the shortest
1690 // such that we minimize our iterations
1705 pSetCoeff0(e,pGetCoeff(b));
1706 if (smIsNegQuot(e, b, c))
1708 lead = pLmDivisibleByNoComp(e, a);
1709 r = smSelectCopy_ExpMultDiv(a, e, b, c);
1714 r = pp_Mult__mm(a, e);
1720 smFindRef(&pa, &res, r);
1730 res = p_Add_q(res, r);
1739 if (!TEST_OPT_NOT_BUCKETS && lb >= SM_MIN_LENGTH_BUCKET)
1741 // use buckets if minimum length is smaller than threshold
1744 // find the last monomial before pa
1748 pNext(append) = res;
1753 while (pNext(append) != pa)
1755 assume(pNext(append) != NULL);
1759 kBucket_pt bucket = kBucketCreate(currRing);
1760 kBucketInit(bucket, pNext(append), 0);
1763 pSetCoeff0(e,pGetCoeff(b));
1764 if (smIsNegQuot(e, b, c))
1767 r = pp_Mult_Coeff_mm_DivSelect_MultDiv(a, lr, e, b, c);
1768 if (pLmDivisibleByNoComp(e, a))
1769 append = kBucket_ExtractLarger_Add_q(bucket, append, r, &lr);
1771 kBucket_Add_q(bucket, r, &lr);
1775 r = pp_Mult__mm(a, e);
1776 append = kBucket_ExtractLarger_Add_q(bucket, append, r, &la);
1779 } while (b != NULL);
1780 pNext(append) = kBucketClear(bucket);
1781 kBucketDestroy(&bucket);
1789 pSetCoeff0(e,pGetCoeff(b));
1790 if (smIsNegQuot(e, b, c))
1792 r = smSelectCopy_ExpMultDiv(a, e, b, c);
1793 if (pLmDivisibleByNoComp(e, a))
1794 smCombineChain(&pa, r);
1800 r = pp_Mult__mm(a, e);
1801 smCombineChain(&pa, r);
1804 } while (b != NULL);
1885 }
while (b !=
NULL);
1974 for (i=
rVar(R); i>0; i--)
2015 }
while (a !=
NULL);
2144 memcpy(r, a,
sizeof(
smprec));
2216 q =
pNext(pp) = a->m;
2239 for(i=
rVar(R); i>0; i--)
2241 if (
p_GetExp(p,i,R) != 0)
return res+1.0;
2256 return res+(float)i;
2291 if (!sw)
return res;
2312 typedef struct smnrec sm_nrec;
2353 void smZeroToredElim();
2360 void smTriangular();
2362 ideal smRes2Ideal();
2380 WerrorS(
"symbol in equation");
2396 WerrorS(
"singular problem for linsolv");
2415 tored =
nrows = smat->rank;
2418 m_row = (smnumber *)
omAlloc0(
sizeof(smnumber)*
i);
2419 wrw = (
int *)
omAlloc(
sizeof(
int)*
i);
2421 wcl = (
int *)
omAlloc(
sizeof(
int)*
i);
2422 m_act = (smnumber *)
omAlloc(
sizeof(smnumber)*
i);
2423 m_res = (smnumber *)
omAlloc0(
sizeof(smnumber)*
i);
2457 this->smZeroToredElim();
2458 if (sing != 0)
return;
2461 this->smRealPivot();
2468 this->smZeroToredElim();
2469 if (sing != 0)
return;
2489 smnumber
s, d,
r = m_row[
nrows];
2492 sol = (number *)
omAlloc0(
sizeof(number)*(crd+1));
2504 sol[
i] =
n_Div(x, m_res[i]->
m,_R->cf);
2518 z =
n_Mult(sol[j], s->m,_R->cf);
2522 x =
n_Sub(y, z,_R->cf);
2535 y =
n_Add(x, sol[i],_R->cf);
2551 sol[
i] =
n_Div(x, d->m,_R->cf);
2591 while ((a!=
NULL) && (a->pos<=tored))
2623 m_act[
act] = m_act[copt];
2635 smnumber c = m_act[
act];
2651 w =
n_Mult(r->m, p,_R->cf);
2661 res->m =
n_Mult(b->m, w,_R->cf);
2663 }
while (b !=
NULL);
2666 if (a->pos < b->pos)
2671 else if (a->pos > b->pos)
2674 res->m =
n_Mult(b->m, w,_R->cf);
2679 hb =
n_Mult(b->m, w,_R->cf);
2680 ha =
n_Add(a->m, hb,_R->cf);
2701 m_act[r->pos] = dumm->n;
2703 }
while (r !=
NULL);
2732 }
while (a->pos < rpiv);
2739 for (i=1; i<
act; i++)
2748 if ((a ==
NULL) || (a->pos > rpiv))
2760 else if (a->pos == rpiv)
2778 smnumber c = m_act[
act];
2785 h->n = m_row[h->pos];
2798 smnumber
r = m_row[rpiv];
2838 if ((a==
NULL) || (a->pos>tored))
2883 smnumber
a = *
r,
b = a->n;
2893 memcpy(r, a,
sizeof(
smnrec));
2947 if ((i == 0) || (i != I->rank-1))
2949 WerrorS(
"wrong dimensions for linsolv");
2954 if(I->m[i-1] ==
NULL)
2956 WerrorS(
"singular input for linsolv");
static FORCE_INLINE number n_Sub(number a, number b, const coeffs r)
return the difference of 'a' and 'b', i.e., a-b
static poly sm_Smpoly2Poly(smpoly, const ring)
static FORCE_INLINE BOOLEAN n_Greater(number a, number b, const coeffs r)
ordered fields: TRUE iff 'a' is larger than 'b'; in Z/pZ: TRUE iff la > lb, where la and lb are the l...
sparse_number_mat(ideal, const ring)
const CanonicalForm int s
ring sm_RingChange(const ring origR, long bound)
const CanonicalForm int const CFList const Variable & y
void sm_SpecialPolyDiv(poly a, poly b, const ring R)
static void sm_ExpMultDiv(poly, const poly, const poly, const ring)
static BOOLEAN p_LmIsConstantComp(const poly p, const ring r)
static poly normalize(poly next_p, ideal add_generators, syStrategy syzstr, int *g_l, int *p_l, int crit_comp)
void kBucketInit(kBucket_pt bucket, poly lm, int length)
static CanonicalForm bound(const CFMatrix &M)
void sm_PolyDiv(poly a, poly b, const ring R)
static poly pp_Mult_Coeff_mm_DivSelect_MultDiv(poly p, int &lp, poly m, poly a, poly b, const ring currRing)
static BOOLEAN p_LmDivisibleByNoComp(poly a, poly b, const ring r)
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
poly sm_CallDet(ideal I, const ring R)
poly prMoveR(poly &p, ring src_r, ring dest_r)
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
#define omFreeSize(addr, size)
static short rVar(const ring r)
#define rVar(r) (r->N)
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
static BOOLEAN smSmaller(poly, poly)
static poly pp_Mult_mm(poly p, poly m, const ring r)
BOOLEAN id_IsConstant(ideal id, const ring r)
test if the ideal has only constant polynomials NOTE: zero ideal/module is also constant ...
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
void WerrorS(const char *s)
static number sm_Cleardenom(ideal, const ring)
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
static void sm_NumberDelete(smnumber *, const ring R)
long sm_ExpBound(ideal m, int di, int ra, int t, const ring currRing)
static number p_SetCoeff(poly p, number n, ring r)
static void p_LmFree(poly p, ring)
static int pLength(poly a)
poly kBucketExtractLm(kBucket_pt bucket)
static long p_SubExp(poly p, int v, long ee, ring r)
#define TEST_OPT_NOT_BUCKETS
static FORCE_INLINE number n_Mult(number a, number b, const coeffs r)
return the product of 'a' and 'b', i.e., a*b
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
static poly sm_Smnumber2Poly(number, const ring)
void smNewBareiss(int, int)
void kBucketDestroy(kBucket_pt *bucket_pt)
long id_RankFreeModule(ideal s, ring lmRing, ring tailRing)
return the maximal component number found in any polynomial in s
static smnumber smNumberCopy(smnumber)
BOOLEAN rComplete(ring r, int force)
this needs to be called whenever a new ring is created: new fields in ring are created (like VarOffse...
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent : the integer VarOffset encodes:
static BOOLEAN p_IsConstant(const poly p, const ring r)
static FORCE_INLINE number n_Add(number a, number b, const coeffs r)
return the sum of 'a' and 'b', i.e., a+b
static void sm_FindRef(poly *, poly *, poly, const ring)
static smpoly sm_Poly2Smpoly(poly, const ring)
ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
static BOOLEAN smCheckSolv(ideal)
static poly pp_Mult_qq(poly p, poly q, const ring r)
BOOLEAN sm_CheckDet(ideal I, int d, BOOLEAN sw, const ring r)
static BOOLEAN sm_HaveDenom(poly, const ring)
ideal idrMoveR(ideal &id, ring src_r, ring dest_r)
static int p_LmCmp(poly p, poly q, const ring r)
static FORCE_INLINE number n_Invers(number a, const coeffs r)
return the multiplicative inverse of 'a'; raise an error if 'a' is not invertible ...
#define SM_MIN_LENGTH_BUCKET
static FORCE_INLINE number n_InpNeg(number n, const coeffs r)
in-place negation of n MUST BE USED: n = n_InpNeg(n) (no copy is returned)
static int si_max(const int a, const int b)
void PrintS(const char *s)
static poly p_Mult_nn(poly p, number n, const ring r)
static BOOLEAN rField_is_Q(const ring r)
static void p_ExpVectorAdd(poly p1, poly p2, const ring r)
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff 'n' represents the zero element.
static void sm_PolyDivN(poly, const number, const ring)
void p_Normalize(poly p, const ring r)
static void p_Delete(poly *p, const ring r)
#define omGetSpecBin(size)
ideal idInit(int idsize, int rank)
initialise an ideal / module
const Variable & v
< [in] a sqrfree bivariate poly
static void p_ExpVectorDiff(poly pr, poly p1, poly p2, const ring r)
static unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
set a single variable exponent : VarOffset encodes the position in p->exp
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
void sm_KillModifiedRing(ring r)
ideal sm_CallSolv(ideal I, const ring R)
static void sm_CombineChain(poly *, poly, const ring)
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
static smnumber sm_Poly2Smnumber(poly, const ring)
static smpoly smElemCopy(smpoly)
sparse_mat(ideal, const ring)
static BOOLEAN p_IsConstantPoly(const poly p, const ring r)
static FORCE_INLINE number n_GetDenom(number &n, const coeffs r)
return the denominator of n (if elements of r are by nature not fractional, result is 1) ...
static FORCE_INLINE BOOLEAN n_Equal(number a, number b, const coeffs r)
TRUE iff 'a' and 'b' represent the same number; they may have different representations.
void rKillModifiedRing(ring r)
ideal idrCopyR(ideal id, ring src_r, ring dest_r)
static void p_Setm(poly p, const ring r)
poly sm_MultDiv(poly a, poly b, const poly c, const ring R)
static void smMinSelect(long *, int, int)
void sm_CallBareiss(ideal I, int x, int y, ideal &M, intvec **iv, const ring R)
static void sm_ExactPolyDiv(poly, poly, const ring)
static poly p_Neg(poly p, const ring r)
static void p_LmDelete(poly p, const ring r)
static float sm_PolyWeight(smpoly, const ring)
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
static FORCE_INLINE BOOLEAN n_GreaterZero(number n, const coeffs r)
ordered fields: TRUE iff 'n' is positive; in Z/pZ: TRUE iff 0 < m <= roundedBelow(p/2), where m is the long representing n in C: TRUE iff (Im(n) != 0 and Im(n) >= 0) or (Im(n) == 0 and Re(n) >= 0) in K(a)/<p(a)>: TRUE iff (n != 0 and (LC(n) > 0 or deg(n) > 0)) in K(t_1, ..., t_n): TRUE iff (LC(numerator(n) is a constant and > 0) or (LC(numerator(n) is not a constant) in Z/2^kZ: TRUE iff 0 < n <= 2^(k-1) in Z/mZ: TRUE iff the internal mpz is greater than zero in Z: TRUE iff n > 0
static poly p_Add_q(poly p, poly q, const ring r)
static poly sm_SelectCopy_ExpMultDiv(poly p, poly m, poly a, poly b, const ring currRing)
#define omFreeBin(addr, bin)
static BOOLEAN rOrd_is_Comp_dp(const ring r)
kBucket_pt kBucketCreate(const ring bucket_ring)
Creation/Destruction of buckets.
static poly p_New(const ring, omBin bin)
static poly p_Init(const ring r, omBin bin)
poly p_Cleardenom(poly p, const ring r)
static void sm_ElemDelete(smpoly *, const ring)
ideal idrCopyR_NoSort(ideal id, ring src_r, ring dest_r)
static FORCE_INLINE int n_Size(number n, const coeffs r)
return a non-negative measure for the complexity of n; return 0 only when n represents zero; (used fo...
void smToIntvec(intvec *)
void Werror(const char *fmt,...)
static poly pp_Mult_Coeff_mm_DivSelect(poly p, const poly m, const ring r)
static long p_MaxComp(poly p, ring lmRing, ring tailRing)
void kBucket_Add_q(kBucket_pt bucket, poly q, int *l)
Add to Bucket a poly ,i.e. Bpoly == q+Bpoly.
static BOOLEAN sm_IsNegQuot(poly, const poly, const poly, const ring)