v.0.9.10

* faster sparse polynomial division and multiplication using Maxheap * Abacus.Heap implementation * System of diophantine equations check additional rows for solution if over-determined system * update examples, tests
foo123 · Jan 28, 2020 · 7ef6b91 · 7ef6b91
1 parent ecc08fe
commit 7ef6b91
Show file tree

Hide file tree

Showing 44 changed files with 818 additions and 727 deletions.
diff --git a/README.md b/README.md
@@ -139,7 +139,7 @@ A variety of combinatorial algorithms, number theory algorithms &amp; statistics
 * [Numerical algorithms for the computation of the Smith normal form of integral matrices, C. Koukouvinos, M. Mitrouli, J. Seberry](https://ro.uow.edu.au/cgi/viewcontent.cgi?referer=https://www.google.com/&httpsredir=1&article=2173&context=infopapers)
 * [Fraction-free matrix factors: new forms for LU and QR factors, Wenqin ZHOU, David J. JEFFREY](http://ftp.cecm.sfu.ca/personal/pborwein/MITACS/papers/FFMatFacs08.pdf)
 * [Algorithms and Data Structures for Sparse Polynomial Arithmetic, M. Asadi, A. Brandt, R. H. C. Moir, M. M. Maza](https://www.researchgate.net/publication/333182217_Algorithms_and_Data_Structures_for_Sparse_Polynomial_Arithmetic)
-* [High Performance Sparse Multivariate Polynomials: Fundamental Data Structures and Algorithms, Alex Brandt (MSc thesis)](https://pdfs.semanticscholar.org/0f93/5035aefc4afce591cd52c507cd7fa35eb061.pdf?_ga=2.149338422.1292196222.1579944113-1940599073.1579944113)
+* [High Performance Sparse Multivariate Polynomials: Fundamental Data Structures and Algorithms, Alex Brandt (MSc thesis)](https://www.semanticscholar.org/paper/High-Performance-Sparse-Multivariate-Polynomials%3A-Brandt/016a97690ecaed04d7a60c1dbf27eb5a96de2dc1)
 * [Algorithms for Normal Forms for Matrices of Polynomials and Ore Polynomials, Howard Cheng (PhD thesis)](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.9.4150&rep=rep1&type=pdf)
 
 ### Example API
@@ -605,7 +605,7 @@ o.dispose()
 * support efficient primality tests and prime sieves **[DONE PARTIALY]**
 * support efficient integer factorization algorithms **[DONE PARTIALY]**
 * support solutions of (systems of) **linear diophantine and linear congruence equations** (with one or many variables) **[DONE]**
-* use sparse representation for polynomials (univariate and multivariate) instead of the, in general, inefficient dense representation **[DONE PARTIALY]**
+* use sparse representation for polynomials (univariate and multivariate) instead of the, in general, inefficient dense representation (and optimise associated arithmetic operations) **[DONE]**
 * support (univariate) polynomial (partial) factorisation, (rational) root finding **[DONE PARTIALY]**
 * support multivariate polynomial, multivariate division/reduction, groebner basis computations (TODO)
 * implement `LLL` algorithm (TODO)

diff --git a/src/js/Abacus.js b/src/js/Abacus.js
@@ -3325,6 +3325,12 @@ function solvediophs( a, b, with_param, with_free_vars )
         }));
     });
     free_vars.symbol = symbol;
+
+    // if over-determined system (m > k)
+    // check if additional rows are satisfied by solution as well
+    for(i=k; i<m; i++)
+        if ( !Expr(solutions.map(function(xj){return xj.mul(a.val[i][j]);})).equ(b[i]) )
+            return null; // no solution
 
     /*
     // solve by successive substitution
@@ -4640,6 +4646,7 @@ function addition_sparse( a, b, do_subtraction, gt )
     // merge terms by efficient merging and produce already sorted order c
     // eg http://www.cecm.sfu.ca/~mmonagan/teaching/TopicsinCA11/johnson.pdf
     // and https://www.researchgate.net/publication/333182217_Algorithms_and_Data_Structures_for_Sparse_Polynomial_Arithmetic
+    // and https://www.semanticscholar.org/paper/High-Performance-Sparse-Multivariate-Polynomials%3A-Brandt/016a97690ecaed04d7a60c1dbf27eb5a96de2dc1
     gt = gt || function(a, b){return a>b;};
     var i = 0, j = 0, k = 0, n1 = a.length, n2 = b.length, c = new Array(n1+n2), res;
     while( i<n1 && j<n2 )
@@ -4681,7 +4688,8 @@ function multiplication_sparse( a, b )
     // merge terms by efficient merging and produce already sorted order c
     // eg http://www.cecm.sfu.ca/~mmonagan/teaching/TopicsinCA11/johnson.pdf
     // and https://www.researchgate.net/publication/333182217_Algorithms_and_Data_Structures_for_Sparse_Polynomial_Arithmetic
-    var k, l, s, t, n1, n2, c, f, max, res, heap;
+    // and https://www.semanticscholar.org/paper/High-Performance-Sparse-Multivariate-Polynomials%3A-Brandt/016a97690ecaed04d7a60c1dbf27eb5a96de2dc1
+    var k, t, n1, n2, c, f, max, heap;
     if ( a.length > b.length ){ t=a; a=b; b=t;} // swap to achieve better performance
     n1 = a.length; n2 = b.length; c = new Array(n1*n2);
     if ( 0<n1 && 0<n2 )
@@ -4691,17 +4699,22 @@ function multiplication_sparse( a, b )
         heap = Heap(array(n1, function(i){
             return [a[i].e+b[0].e, a[i].c.mul(b[0].c), i];
         }), "max", function(a, b){
-            return a[0]<b[0] ? -1 : (a[0]>b[0] ? 1 : 0);
+            return a[0]-b[0];
         });
         f = array(n1, 0);
         while( max=heap.peek() )
         {
-            if ( c[k].e!==max[0] && !c[k].c.equ(0) ) c[++k] = Coeff(0, max[0]);
-            heap.pop();
+            if ( c[k].e!==max[0] )
+            {
+                if ( !c[k].c.equ(0) ) c[++k] = Coeff(0, -1);
+                c[k].e = max[0];
+            }
             c[k].c = c[k].c.add(max[1]);
             f[max[2]]++;
-            if ( f[max[2]] < n2 ) heap.push([a[max[2]].e+b[f[max[2]]].e, a[max[2]].c.mul(b[f[max[2]]].c), max[2]]);
+            if ( f[max[2]] < n2 ) heap.replace([a[max[2]].e+b[f[max[2]]].e, a[max[2]].c.mul(b[f[max[2]]].c), max[2]]);
+            else heap.pop();
         }
+        heap.dispose();
         if ( c.length > k+1 ) c.length = k+1; // truncate if needed
     }
     return c;
@@ -8495,7 +8508,7 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
     ,div: function( x, q_and_r ) {
         var self = this, Arithmetic = Abacus.Arithmetic,
             O = Arithmetic.O, I = Arithmetic.I,
-            q, r, d, diff, diff0;
+            q, r, d, diff, diff0, k, res, Q, a, na, b, nb, heap;
 
         if ( x instanceof Complex ) x = x.real;
         else if ( Arithmetic.isNumber(x) ) x = Rational(x);
@@ -8518,9 +8531,10 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
                 }), self.symbol);
                 return true===q_and_r ? [q, Polynomial([], self.symbol)] : q;
             }
+
             // polynomial long division
             // TODO: make it faster
-            r = Polynomial(self);
+            /*r = Polynomial(self);
             diff = r.deg()-x.deg();
             if ( 0 <= diff )
             {
@@ -8539,7 +8553,56 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
             {
                 q = [];
             }
+            q = Polynomial(q, self.symbol);*/
+
+            // sparse polynomial long division
+            // https://www.semanticscholar.org/paper/High-Performance-Sparse-Multivariate-Polynomials%3A-Brandt/016a97690ecaed04d7a60c1dbf27eb5a96de2dc1
+            a = self.coeff; na = a.length; b = x.coeff; nb = b.length;
+            heap = Heap([], "max", function(a,b){return a.coeff.e-b.coeff.e;});
+            q = []; r = []; k = 0;
+            while( (d=heap.peek()) || k<na )
+            {
+                if ( (null == d) || (k<na && d.coeff.e<a[k].e) )
+                {
+                    res = a[k].clone();
+                    k++;
+                }
+                else if ( k<na && d.coeff.e===a[k].e )
+                {
+                    res = Coeff(a[k].c.sub(d.coeff.c), d.coeff.e);
+                    if ( nb>d.n )
+                        heap.replace({coeff:Coeff(d.Q.c.mul(b[d.n].c), d.Q.e+b[d.n].e), n:d.n+1, Q:d.Q});
+                    else
+                        heap.pop();
+                    k++;
+
+                    if ( res.c.equ(O) ) continue; // zero coefficient, skip
+                }
+                else
+                {
+                    res = Coeff(d.coeff.c.neg(), d.coeff.e);
+                    if ( nb>d.n )
+                        heap.replace({coeff:Coeff(d.Q.c.mul(b[d.n].c), d.Q.e+b[d.n].e), n:d.n+1, Q:d.Q});
+                    else
+                        heap.pop();
+                }
+
+                if ( b[0].e<=res.e ) // if b[0] divides res
+                {
+                    Q = Coeff(res.c.div(b[0].c), res.e-b[0].e);
+                    q = addition_sparse(q, [Q]);
+                    if ( nb>1 )
+                        heap.push({coeff:Coeff(Q.c.mul(b[1].c), Q.e+b[1].e), n:2, Q:Q});
+                }
+                else
+                {
+                    r = addition_sparse(r, [res]);
+                }
+            }
+            heap.dispose();
             q = Polynomial(q, self.symbol);
+            r = Polynomial(r, self.symbol);
+
             // return both quotient and remainder if requested
             return true===q_and_r ? [q, r] : q;
         }

diff --git a/src/js/Abacus.min.js b/src/js/Abacus.min.js
diff --git a/test/combination_subsets.txt b/test/combination_subsets.txt
@@ -1,4 +1,4 @@
-Abacus.Subsets via Abacus.Combinations (VERSION = 0.9.9)
+Abacus.Subsets via Abacus.Combinations (VERSION = 0.9.10)
 ---
 o = Abacus.CombinatorialIterator([Abacus.Combination(5,0),Abacus.Combination(5,1),Abacus.Combination(5,2),Abacus.Combination(5,3),Abacus.Combination(5,4),Abacus.Combination(5,5)])
 o.total()
@@ -243,38 +243,38 @@ o.order("colex,reversed")
 [ 0 ]
 []
 o.order("random")
-[ 1 ]
-[ 0, 2 ]
-[ 2, 1 ]
-[ 0 ]
+[ 0, 1, 2, 3 ]
+[ 4 ]
+[ 2, 3, 4 ]
+[ 1, 3, 4 ]
+[ 0, 1, 4 ]
+[ 3 ]
 [ 0, 1, 2 ]
 [ 1 ]
 [ 3, 1 ]
-[ 1, 4 ]
-[ 1, 3 ]
-[ 0, 1, 2, 4 ]
-[ 1, 0 ]
-[ 2 ]
-[ 2 ]
+[ 3, 1 ]
+[ 1, 2, 3, 4 ]
 [ 0, 1, 4 ]
-[ 0, 1, 3, 4 ]
+[ 1, 2 ]
+[ 4 ]
+[ 0, 3 ]
+[ 1, 2 ]
 [ 1, 2, 4 ]
-[ 0, 1, 2, 3 ]
 [ 0, 2, 3, 4 ]
-[ 1 ]
-[ 1, 2 ]
-[ 2, 0 ]
-[ 1, 3, 4 ]
-[ 0, 3, 4 ]
+[ 3 ]
+[ 0, 1, 2, 3, 4 ]
 [ 0, 2, 3, 4 ]
 [ 1, 3 ]
+[ 2, 3 ]
+[ 1, 2 ]
 [ 0, 1, 3, 4 ]
-[ 0, 1 ]
-[ 2, 3, 4 ]
-[ 0, 1, 2, 4 ]
-[ 4, 3 ]
-[ 1 ]
-[ 0, 1, 2, 4 ]
+[ 1, 3, 4 ]
+[ 2 ]
+[ 1, 2, 3, 4 ]
+[ 0, 1, 2, 3, 4 ]
+[ 0, 2, 3, 4 ]
+[ 0, 2, 3, 4 ]
+[ 0, 2, 3, 4 ]
 o.random()
 [ 3 ]
 o.dispose()
diff --git a/test/combinations.txt b/test/combinations.txt
@@ -1,4 +1,4 @@
-Abacus.Combinations (VERSION = 0.9.9)
+Abacus.Combinations (VERSION = 0.9.10)
 ---
 o = Abacus.Combination(6,3)
 o.total()
@@ -159,28 +159,28 @@ o.order("colex,reversed")
 [ 0, 1, 3 ]
 [ 0, 1, 2 ]
 o.order("random")
+[ 0, 2, 5 ]
 [ 1, 4, 5 ]
-[ 0, 2, 4 ]
 [ 2, 3, 4 ]
-[ 1, 2, 5 ]
-[ 0, 2, 3 ]
-[ 0, 2, 5 ]
-[ 0, 3, 4 ]
-[ 0, 3, 5 ]
-[ 2, 3, 5 ]
-[ 0, 1, 4 ]
-[ 1, 2, 3 ]
-[ 0, 1, 5 ]
 [ 0, 4, 5 ]
+[ 0, 2, 4 ]
+[ 1, 3, 5 ]
 [ 2, 4, 5 ]
+[ 2, 3, 5 ]
 [ 0, 1, 3 ]
-[ 1, 3, 5 ]
 [ 1, 3, 4 ]
+[ 0, 3, 5 ]
 [ 1, 2, 4 ]
+[ 0, 3, 4 ]
+[ 0, 1, 5 ]
+[ 1, 2, 5 ]
 [ 0, 1, 2 ]
 [ 3, 4, 5 ]
+[ 0, 2, 3 ]
+[ 0, 1, 4 ]
+[ 1, 2, 3 ]
 o.random()
-[ 2, 3, 4 ]
+[ 0, 1, 3 ]
 o.order("colex").range(-5, -1)
 [ 2, 3, 5 ]
 [ 0, 4, 5 ]

diff --git a/test/combinations_repeats.txt b/test/combinations_repeats.txt
@@ -1,4 +1,4 @@
-Abacus.CombinationRepeats (VERSION = 0.9.9)
+Abacus.CombinationRepeats (VERSION = 0.9.10)
 ---
 o = Abacus.Combination(6,3,{type:"repeated"})
 o.total()
@@ -411,64 +411,64 @@ o.order("colex,reversed")
 [ 0, 0, 1 ]
 [ 0, 0, 0 ]
 o.order("random")
-[ 2, 3, 5 ]
-[ 2, 3, 3 ]
-[ 3, 3, 4 ]
-[ 4, 5, 5 ]
-[ 2, 2, 3 ]
-[ 2, 3, 4 ]
+[ 0, 0, 3 ]
 [ 1, 1, 5 ]
-[ 0, 3, 3 ]
-[ 1, 1, 1 ]
-[ 0, 0, 4 ]
-[ 2, 4, 5 ]
-[ 1, 3, 4 ]
-[ 1, 2, 2 ]
 [ 2, 2, 5 ]
-[ 0, 3, 4 ]
+[ 2, 3, 3 ]
+[ 0, 1, 1 ]
+[ 0, 2, 2 ]
+[ 0, 1, 5 ]
+[ 4, 4, 4 ]
+[ 0, 1, 4 ]
+[ 0, 1, 2 ]
+[ 3, 5, 5 ]
+[ 4, 4, 5 ]
+[ 0, 5, 5 ]
+[ 2, 3, 5 ]
+[ 0, 3, 3 ]
 [ 0, 2, 3 ]
-[ 0, 0, 5 ]
-[ 1, 2, 4 ]
+[ 3, 3, 5 ]
 [ 2, 4, 4 ]
+[ 3, 4, 5 ]
+[ 0, 0, 2 ]
 [ 2, 2, 2 ]
-[ 3, 3, 5 ]
-[ 0, 2, 2 ]
+[ 0, 2, 4 ]
 [ 1, 2, 3 ]
-[ 0, 0, 3 ]
-[ 0, 4, 5 ]
+[ 1, 5, 5 ]
+[ 0, 0, 0 ]
 [ 1, 4, 5 ]
-[ 1, 1, 4 ]
+[ 1, 2, 4 ]
+[ 2, 2, 3 ]
 [ 0, 1, 3 ]
-[ 0, 1, 1 ]
-[ 1, 1, 3 ]
+[ 3, 4, 4 ]
+[ 1, 1, 4 ]
+[ 0, 0, 4 ]
+[ 0, 2, 5 ]
 [ 2, 2, 4 ]
+[ 1, 2, 2 ]
+[ 5, 5, 5 ]
+[ 2, 3, 4 ]
+[ 0, 0, 5 ]
+[ 3, 3, 4 ]
+[ 0, 3, 4 ]
 [ 0, 0, 1 ]
+[ 0, 3, 5 ]
+[ 1, 4, 4 ]
+[ 1, 3, 3 ]
+[ 4, 5, 5 ]
+[ 0, 4, 4 ]
+[ 1, 1, 3 ]
+[ 0, 4, 5 ]
+[ 1, 1, 1 ]
+[ 3, 3, 3 ]
 [ 1, 2, 5 ]
+[ 1, 3, 4 ]
 [ 2, 5, 5 ]
-[ 3, 4, 5 ]
-[ 3, 4, 4 ]
-[ 3, 3, 3 ]
-[ 1, 3, 3 ]
-[ 0, 1, 5 ]
-[ 1, 5, 5 ]
+[ 2, 4, 5 ]
 [ 1, 1, 2 ]
-[ 3, 5, 5 ]
-[ 0, 5, 5 ]
-[ 4, 4, 4 ]
-[ 1, 4, 4 ]
-[ 4, 4, 5 ]
-[ 0, 2, 5 ]
-[ 0, 0, 2 ]
-[ 0, 1, 4 ]
 [ 1, 3, 5 ]
-[ 0, 4, 4 ]
-[ 5, 5, 5 ]
-[ 0, 1, 2 ]
-[ 0, 0, 0 ]
-[ 0, 3, 5 ]
-[ 0, 2, 4 ]
 o.random()
-[ 4, 4, 4 ]
+[ 0, 3, 3 ]
 o.order("colex").range(-5, -1)
 [ 1, 5, 5 ]
 [ 2, 5, 5 ]

diff --git a/test/complex.txt b/test/complex.txt
@@ -1,4 +1,4 @@
-Abacus.Complex (VERSION = 0.9.9)
+Abacus.Complex (VERSION = 0.9.10)
 ---
 o=Abacus.Complex()
 o.toString()