Skip to content

Commit 47502eb

Browse files
author
jkuipers
committed
removed hashing of terms in polynomial division, since this trick
doesn't work there properly.
1 parent f75f8f6 commit 47502eb

2 files changed

Lines changed: 49 additions & 129 deletions

File tree

‎sources/poly.cc‎

Lines changed: 48 additions & 128 deletions
Original file line numberDiff line numberDiff line change
@@ -1487,11 +1487,15 @@ void poly::divmod_univar (const poly &a, const poly &b, poly &q, poly &r, int va
14871487
* A heap element h is formatted as follows:
14881488
* - h[0] = index in a
14891489
* - h[1] = index in b
1490-
* - h[2] = hash code (-1 if no hash is used)
1490+
* - h[2] = -1 (no hash is used)
14911491
* - h[3] = length of coefficient with sign
14921492
* - h[4...4+AN.poly_num_vars-1] = powers
14931493
* - h[4+AN.poly_num_vars...4+h[3]-1] = coefficient
14941494
*
1495+
* Note: the hashing trick as in multiplication cannot be used
1496+
* easily, since there is no tight upperbound on the exponents in
1497+
* the answer.
1498+
*
14951499
* For details, see M. Monagan, "Polynomial Division using Dynamic
14961500
* Array, Heaps, and Packed Exponent Vectors"
14971501
*/
@@ -1533,42 +1537,7 @@ void poly::divmod_heap (const poly &a, const poly &b, poly &q, poly &r, bool onl
15331537
// allocate heap
15341538
int nb=b.number_of_terms();
15351539

1536-
// determine maximum power in variables
1537-
WORD *maxpower = AT.WorkPointer;
1538-
AT.WorkPointer += AN.poly_num_vars;
1539-
WORD *maxpowera = AT.WorkPointer;
1540-
AT.WorkPointer += AN.poly_num_vars;
1541-
WORD *maxpowerb = AT.WorkPointer;
1542-
AT.WorkPointer += AN.poly_num_vars;
1543-
1544-
for (int i=0; i<AN.poly_num_vars; i++)
1545-
maxpowera[i] = maxpowerb[i] = 0;
1546-
1547-
for (int ai=1; ai<a[0]; ai+=a[ai])
1548-
for (int j=0; j<AN.poly_num_vars; j++)
1549-
maxpowera[j] = MaX(maxpowera[j], a[ai+1+j]);
1550-
1551-
for (int bi=1; bi<b[0]; bi+=b[bi])
1552-
for (int j=0; j<AN.poly_num_vars; j++)
1553-
maxpowerb[j] = MaX(maxpowerb[j], b[bi+1+j]);
1554-
1555-
for (int i=0; i<AN.poly_num_vars; i++)
1556-
maxpower[i] = MaX(maxpowera[i],maxpowerb[i]);
1557-
1558-
// if PROD(max.power) small, allocate hash table
1559-
bool use_hash = true;
1560-
int nhash = 1;
1561-
1562-
for (int i=0; i<AN.poly_num_vars; i++) {
1563-
if (nhash > POLY_MAX_HASH_SIZE / (maxpower[i]+1)) {
1564-
nhash = 1;
1565-
use_hash = false;
1566-
break;
1567-
}
1568-
nhash *= maxpower[i]+1;
1569-
}
1570-
1571-
WantAddPointers(nb+nhash);
1540+
WantAddPointers(nb);
15721541
WORD **heap = AT.pWorkSpace + AT.pWorkPointer;
15731542

15741543
for (int i=0; i<nb; i++)
@@ -1579,10 +1548,6 @@ void poly::divmod_heap (const poly &a, const poly &b, poly &q, poly &r, bool onl
15791548
heap[0][2] = -1;
15801549
WCOPY(&heap[0][3], &a[1], a[1]);
15811550
heap[0][3] = a[a[1]];
1582-
1583-
WORD **hash = AT.pWorkSpace + AT.pWorkPointer + nb;
1584-
for (int i=0; i<nhash; i++)
1585-
hash[i] = NULL;
15861551

15871552
int qi=1, ri=1;
15881553

@@ -1613,8 +1578,6 @@ void poly::divmod_heap (const poly &a, const poly &b, poly &q, poly &r, bool onl
16131578
pop_heap(BHEAD heap, nheap--);
16141579
p = heap[nheap];
16151580

1616-
if (p[2]!=-1) hash[p[2]] = NULL;
1617-
16181581
if (t[0] == -1) {
16191582
WCOPY(t, p, (5+ABS(p[3])+AN.poly_num_vars));
16201583
}
@@ -1650,97 +1613,55 @@ void poly::divmod_heap (const poly &a, const poly &b, poly &q, poly &r, bool onl
16501613
insert.pop_back();
16511614
}
16521615

1653-
// add elements to the heap
1654-
while (true) {
1655-
// prepare the element
1656-
if (p[1]==0) {
1657-
p[0] += a[p[0]];
1658-
if (p[0]==a[0]) break;
1659-
WCOPY(&p[3], &a[p[0]], a[p[0]]);
1660-
p[3] = p[2+p[3]];
1661-
}
1662-
else {
1663-
if (!this_insert)
1664-
p[1] += q[p[1]];
1665-
this_insert = false;
1666-
1667-
if (p[1]==qi) { s++; break; }
1668-
1669-
for (int i=0; i<AN.poly_num_vars; i++)
1670-
p[4+i] = b[p[0]+1+i] + q[p[1]+1+i];
1671-
1672-
// if both polynomials are modulo p^1, use integer calculus
1673-
if (both_mod_small) {
1674-
p[4+AN.poly_num_vars] = ((LONG)b[p[0]+1+AN.poly_num_vars]*b[p[0]+b[p[0]]-1]*
1675-
q[p[1]+1+AN.poly_num_vars]*q[p[1]+q[p[1]]-1]) % q.modp;
1676-
if (p[4+AN.poly_num_vars]==0)
1677-
p[3]=0;
1678-
else {
1679-
if (p[4+AN.poly_num_vars] > +q.modp/2) p[4+AN.poly_num_vars] -= q.modp;
1680-
if (p[4+AN.poly_num_vars] < -q.modp/2) p[4+AN.poly_num_vars] += q.modp;
1681-
p[3] = SGN(p[4+AN.poly_num_vars]);
1682-
p[4+AN.poly_num_vars] = ABS(p[4+AN.poly_num_vars]);
1683-
}
1684-
}
1616+
// prepare the element to add to the heap
1617+
if (p[1]==0) {
1618+
p[0] += a[p[0]];
1619+
if (p[0]==a[0]) break;
1620+
WCOPY(&p[3], &a[p[0]], a[p[0]]);
1621+
p[3] = p[2+p[3]];
1622+
}
1623+
else {
1624+
if (!this_insert)
1625+
p[1] += q[p[1]];
1626+
this_insert = false;
1627+
1628+
if (p[1]==qi) { s++; break; }
1629+
1630+
for (int i=0; i<AN.poly_num_vars; i++)
1631+
p[4+i] = b[p[0]+1+i] + q[p[1]+1+i];
1632+
1633+
// if both polynomials are modulo p^1, use integer calculus
1634+
if (both_mod_small) {
1635+
p[4+AN.poly_num_vars] = ((LONG)b[p[0]+1+AN.poly_num_vars]*b[p[0]+b[p[0]]-1]*
1636+
q[p[1]+1+AN.poly_num_vars]*q[p[1]+q[p[1]]-1]) % q.modp;
1637+
if (p[4+AN.poly_num_vars]==0)
1638+
p[3]=0;
16851639
else {
1686-
// otherwise, use form long calculus
1687-
MulLong((UWORD *)&b[p[0]+1+AN.poly_num_vars], b[p[0]+b[p[0]]-1],
1688-
(UWORD *)&q[p[1]+1+AN.poly_num_vars], q[p[1]+q[p[1]]-1],
1689-
(UWORD *)&p[4+AN.poly_num_vars], &p[3]);
1690-
if (q.modp!=0) TakeNormalModulus((UWORD *)&p[4+AN.poly_num_vars], &p[3],
1691-
modq, nmodq, NOUNPACK);
1640+
if (p[4+AN.poly_num_vars] > +q.modp/2) p[4+AN.poly_num_vars] -= q.modp;
1641+
if (p[4+AN.poly_num_vars] < -q.modp/2) p[4+AN.poly_num_vars] += q.modp;
1642+
p[3] = SGN(p[4+AN.poly_num_vars]);
1643+
p[4+AN.poly_num_vars] = ABS(p[4+AN.poly_num_vars]);
16921644
}
1693-
1694-
p[3] *= -1;
1695-
}
1696-
1697-
// with hashing, calculate hash value
1698-
if (use_hash) {
1699-
p[2] = 0;
1700-
for (int i=0; i<AN.poly_num_vars; i++)
1701-
p[2] = (maxpower[i]+1)*p[2] + p[4+i];
17021645
}
17031646
else {
1704-
p[2] = -1;
1647+
// otherwise, use form long calculus
1648+
MulLong((UWORD *)&b[p[0]+1+AN.poly_num_vars], b[p[0]+b[p[0]]-1],
1649+
(UWORD *)&q[p[1]+1+AN.poly_num_vars], q[p[1]+q[p[1]]-1],
1650+
(UWORD *)&p[4+AN.poly_num_vars], &p[3]);
1651+
if (q.modp!=0) TakeNormalModulus((UWORD *)&p[4+AN.poly_num_vars], &p[3],
1652+
modq, nmodq, NOUNPACK);
17051653
}
1706-
1707-
// add it to a heap element if possible, otherwise push it
17081654

1709-
if (!use_hash || (use_hash && hash[p[2]] == NULL)) {
1710-
if (use_hash) hash[p[2]] = p;
1711-
swap (heap[nheap],p);
1712-
push_heap(BHEAD heap, ++nheap);
1713-
break;
1714-
}
1715-
else {
1716-
WORD *h = hash[p[2]];
1717-
// if both polynomials are modulo p^1, use integer calculus
1718-
if (both_mod_small) {
1719-
h[4+AN.poly_num_vars] = ((LONG)p[4+AN.poly_num_vars]*p[3] + h[4+AN.poly_num_vars]*h[3]) % q.modp;
1720-
if (h[4+AN.poly_num_vars]==0)
1721-
h[3]=0;
1722-
else {
1723-
if (h[4+AN.poly_num_vars] > +q.modp/2) h[4+AN.poly_num_vars] -= q.modp;
1724-
if (h[4+AN.poly_num_vars] < -q.modp/2) h[4+AN.poly_num_vars] += q.modp;
1725-
h[3] = SGN(h[4+AN.poly_num_vars]);
1726-
h[4+AN.poly_num_vars] = ABS(h[4+AN.poly_num_vars]);
1727-
}
1728-
}
1729-
else {
1730-
// otherwise, use form long calculus
1731-
AddLong ((UWORD *)&p[4+AN.poly_num_vars], p[3],
1732-
(UWORD *)&h[4+AN.poly_num_vars], h[3],
1733-
(UWORD *)&h[4+AN.poly_num_vars], &h[3]);
1734-
if (q.modp!=0) TakeNormalModulus((UWORD *)&h[4+AN.poly_num_vars], &h[3],
1735-
modq, nmodq, NOUNPACK);
1736-
}
1737-
1738-
if (h[1]<p[1]) {
1739-
swap(h[0],p[0]);
1740-
swap(h[1],p[1]);
1741-
}
1742-
}
1655+
p[3] *= -1;
17431656
}
1657+
1658+
// no hashing
1659+
p[2] = -1;
1660+
1661+
// add it to a heap element if possible, otherwise push it
1662+
1663+
swap (heap[nheap],p);
1664+
push_heap(BHEAD heap, ++nheap);
17441665
}
17451666
while (t[0]==-1 || (nheap>0 && monomial_compare(BHEAD heap[0]+3, t+3)==0));
17461667

@@ -1826,7 +1747,6 @@ void poly::divmod_heap (const poly &a, const poly &b, poly &q, poly &r, bool onl
18261747
NumberFree(t,"poly::div_heap");
18271748

18281749
if (q.modp!=0) NumberFree(ltbinv,"poly::div_heap");
1829-
AT.WorkPointer -= AN.poly_num_vars;
18301750
}
18311751

18321752
/*

‎sources/poly.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -48,7 +48,7 @@ extern "C" {
4848

4949
// maximum size of the hash table used for multiplication and division
5050

51-
const int POLY_MAX_HASH_SIZE = MiN(1<<20, 1<<(BITSINWORD-2));
51+
const int POLY_MAX_HASH_SIZE = MiN(1<<20, MAXPOSITIVE);
5252

5353
class poly {
5454

0 commit comments

Comments
 (0)