Skip to content

Commit

Permalink
v.0.9.8 contd
Browse files Browse the repository at this point in the history 
* Polynomial functional composition method
* Find multiple rational roots of Polys if present (checking derivatives as well)
* Fix Poly primitive/content factorisation to use LCM of denominators
* Minor adjustments in limits for different methods of factorisation and testing
* update examples tests
  • Loading branch information
@foo123
foo123 committed Jan 13, 2020
1 parent 0d28129 commit bb3cd98
Show file tree
Hide file tree
Showing 4 changed files with 137 additions and 31 deletions.
82 changes: 63 additions & 19 deletions src/js/Abacus.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1730,12 +1730,12 @@ function is_prime( n )
//n = Arithmetic.abs(/*N(*/n/*)*/);
ndigits = Arithmetic.digits(n).length;
// try to use fastest algorithm based on size of number (number of digits)
if ( ndigits <= 8 )
if ( ndigits <= 6 )
{
// deterministic test
return wheel_trial_div_test(n);
}
else if ( ndigits <= 1000 )
else if ( ndigits <= 20 )
{
// deterministic test
/*
Expand All @@ -1746,11 +1746,11 @@ function is_prime( n )
If n < 3474749660383 is a 2, 3, 5, 7, 11 and 13-SPRP, then n is prime [Jaeschke93].
If n < 341550071728321 is a 2, 3, 5, 7, 11, 13 and 17-SPRP, then n is prime [Jaeschke93].
*/
/*if ( Arithmetic.lt(n, N(1373653)) )
if ( Arithmetic.lt(n, N(1373653)) )
return miller_rabin_test(n, [two, N(3)]);
else if ( Arithmetic.lt(n, N("25326001")) )
return miller_rabin_test(n, [two, N(3), N(5)]);
else*/ if ( Arithmetic.lt(n, N("25000000000")) )
else if ( Arithmetic.lt(n, N("25000000000")) )
return Arithmetic.equ(n, N("3215031751")) ? false : miller_rabin_test(n, [two, N(3), N(5), N(7)]);
else if ( Arithmetic.lt(n, N("2152302898747")) )
return miller_rabin_test(n, [two, N(3), N(5), N(7), N(11)]);
Expand Down Expand Up @@ -2704,14 +2704,13 @@ function polyxgcd( /* args */ )
}
function divisors( n, as_generator )
{
// time+space O(sqrt(n)) to find all distinct divisors of n (including 1 and n itself)
var Arithmetic = Abacus.Arithmetic, O = Arithmetic.O, I = Arithmetic.I,
list = null, D2 = null, D1 = null, L1 = 0, L2 = 0, node, sqrn, i, n_i, next, factors;
//n = Arithmetic.num(n);
n = Arithmetic.abs(n);
if ( true===as_generator )
{
if ( Arithmetic.gte(n, 1000000) )
if ( Arithmetic.gte(n, 100000) )
{
// for very large numbers,
// compute divisors through prime factorisation
Expand All @@ -2727,6 +2726,7 @@ function divisors( n, as_generator )
}
else
{
// time+space O(sqrt(n)) to find all distinct divisors of n (including 1 and n itself)
sqrn = isqrt(n);
i = I; next = null;
// return iterator/generator
Expand Down Expand Up @@ -2766,6 +2766,7 @@ function divisors( n, as_generator )
}
else
{
// time+space O(sqrt(n)) to find all distinct divisors of n (including 1 and n itself)
sqrn = isqrt(n);
for (i=I; Arithmetic.lte(i,sqrn); i=Arithmetic.add(i,I))
{
Expand Down Expand Up @@ -3549,8 +3550,8 @@ function factorial( n, m )
// simple factorial = F(n) = n F(n-1) = n!
if ( 12 >= n ) return 0 > n ? O : NUM(([1,1,2,6,24,120,720,5040,40320,362880,3628800,39916800,479001600 /*MAX: 2147483647*/])[n]);

// for very large factorials, use the prime factorisation of n!
if ( 500 <= n ) return prime_factorial(Nn);
// for large factorials, use the prime factorisation of n!
if ( 600 <= n ) return prime_factorial(Nn);

key = String(n)/*+'!'*/;
if ( null == factorial.mem1[key] )
Expand Down Expand Up @@ -7101,11 +7102,11 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
,primitive: function( and_content ) {
// factorise into content and primitive part
// https://en.wikipedia.org/wiki/Factorization_of_polynomials#Primitive_part%E2%80%93content_factorization
var self = this, Arithmetic = Abacus.Arithmetic, coeff = self.coeff, coeffp, M, content;
var self = this, Arithmetic = Abacus.Arithmetic, coeff = self.coeff, coeffp, LCM, content;
if ( null == self._prim )
{
M = coeff.reduce(function(M, c){return Arithmetic.mul(M, c.den);}, Arithmetic.I);
coeffp = coeff.map(function(c){return c.mul(M).num;});
LCM = lcm(coeff.map(function(c){return c.den;}));
coeffp = coeff.map(function(c){return c.mul(LCM).integer();});
content = gcd(coeffp);
coeffp = coeffp.map(function(c){return Arithmetic.div(c, content);});
// make positive lead
Expand All @@ -7114,7 +7115,7 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
coeffp = coeffp.map(function(c){return Arithmetic.neg(c);});
content = Arithmetic.neg(content);
}
self._prim = [Polynomial(coeffp, self.symbol), Rational(content, M)];
self._prim = [Polynomial(coeffp, self.symbol), Rational(content, LCM)];
}
return true===and_content ? self._prim.slice() : self._prim[0];
}
Expand All @@ -7123,7 +7124,7 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
// https://en.wikipedia.org/wiki/Rational_root_theorem
// https://en.wikipedia.org/wiki/Gauss%27s_lemma_(polynomial)
var self = this, Arithmetic = Abacus.Arithmetic, O = Arithmetic.O,
roots = [], primitive,c, i, d0, dn, iter, comb, root;
roots = [], primitive,c, p, i, d0, dn, iter, comb, root, nroot, found;

if ( 0 === self.deg() ) return roots; // no roots for constant polynomials

Expand All @@ -7137,17 +7138,37 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
if ( i+1<c.length )
{
// try all possible rational divisors of c_0(excluding trivial zero terms) and c_n
d0 = divisors(c[i].num); dn = divisors(c[c.length-1].num);
iter = divisors(c[i].num, true);
d0 = iter.get(); iter.dispose();
iter = divisors(c[c.length-1].num, true);
dn = iter.get(); iter.dispose();

iter = Tensor([d0.length, dn.length]);
while(iter.hasNext())
{
comb = iter.next();
// try positive root
// positive root
root = Rational(d0[comb[0]], dn[comb[1]]);
if ( primitive.valueOf(root).equ(O) ) roots.push(root);
// try negative root
root = root.neg();
if ( primitive.valueOf(root).equ(O) ) roots.push(root);
// negative root
nroot = Rational(Arithmetic.neg(d0[comb[0]]), dn[comb[1]]);
p = primitive; found = true;
while( found && (0<p.deg()) )
{
found = false;
// try positive root
if ( p.valueOf(root).equ(O) )
{
roots.push(root);
found = true;
}
// try negative root
if ( p.valueOf(nroot).equ(O) )
{
roots.push(nroot);
found = true;
}
if ( found ) p = p.d(); // get derivative to check if roots are multiple
}
}
iter.dispose();
}
Expand Down Expand Up @@ -7295,6 +7316,29 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
return pow;
}
}
,compose: function( x ) {
// functionaly compose one polynomial with another. ie result = P(Q(x))
var self = this, Arithmetic = Abacus.Arithmetic, O = Arithmetic.O, px, c, i;
if ( x instanceof Matrix ) return null;
if ( x instanceof Term ) x = Expr(x);
if ( x instanceof Expr ) x = Polynomial.fromExpr(x, self.symbol);
if ( x instanceof Complex ) x = x.real;
if ( Arithmetic.isNumber(x) || (x instanceof Rational) )
{
return Polynomial([self.valueOf(x)], self.symbol);
}
else if ( x instanceof Polynomial )
{
// Composition through variation of Horner's algorithm for fast evaluation
if ( 0 === self.deg() ) return self.clone();
if ( 0 === x.deg() ) return Polynomial([self.valueOf(x.c())], x.symbol);
c = self.coeff; i = c.length-1;
px = Polynomial([c[i]], x.symbol);
while(i--) px = Polynomial([c[i]], x.symbol).add(px.mul(x));
return px;
}
return self;
}
,shift: function( s ) {
// shift <-> equivalent to multiplication/division by a monomial x^s
var self = this, Arithmetic = Abacus.Arithmetic;
Expand Down
2 changes: 1 addition & 1 deletion src/js/Abacus.min.js
View file

Large diffs are not rendered by default.

57 changes: 49 additions & 8 deletions test/polynomials.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,15 +10,27 @@ function check_div( n, d, q, r )
var nn = q.mul(d).add(r);
console.log('('+n.toString()+')/('+d.toString()+')=('+d.toString()+')*('+q.toString()+')+('+r.toString()+')='+nn.toString(), nn.equ(n));
}
function check_xgcd( args, gcd )
function check_xgcd( args, gcd, tex )
{
var out = '', res = Abacus.Polynomial(0);
for(i=0; i<args.length; i++)
if ( tex )
{
out += (out.length ? ' + ' : '') + '('+args[i].toString()+')'+'('+gcd[i+1].toString()+')';
res = res.add(args[i].mul(gcd[i+1]));
for(i=0; i<args.length; i++)
{
out += (out.length ? ' + ' : '') + '('+args[i].toTex()+')'+'('+gcd[i+1].toTex()+')';
res = res.add(args[i].mul(gcd[i+1]));
}
out += ' = '+gcd[0].toTex();
}
else
{
for(i=0; i<args.length; i++)
{
out += (out.length ? ' + ' : '') + '('+args[i].toString()+')'+'('+gcd[i+1].toString()+')';
res = res.add(args[i].mul(gcd[i+1]));
}
out += ' = '+gcd[0].toString();
}
out += ' = '+gcd[0].toString();
console.log(out, res.equ(gcd[0]));
}
var o, d, qr, args;
Expand Down Expand Up @@ -87,6 +99,14 @@ echo('o.shift(1)');
echo(o.shift(1).toString());
echo('o.shift(-1)');
echo(o.shift(-1).toString());
echo('o.compose(Abacus.Polynomial([1]))');
echo(o.compose(Abacus.Polynomial([1])).toString());
echo('o.compose(Abacus.Polynomial([0,1]))');
echo(o.compose(Abacus.Polynomial([0,1])).toString());
echo('o.compose(Abacus.Polynomial([1,1]))');
echo(o.compose(Abacus.Polynomial([1,1])).toString());
echo('Abacus.Polynomial([1,1,1]).compose(o)');
echo(Abacus.Polynomial([1,1,1]).compose(o).toString());
echo('o.pow(0)');
echo(o.pow(0).toString());
echo('o.pow(1)');
Expand Down Expand Up @@ -117,6 +137,15 @@ echo('o.dispose()');
o.dispose();
echo('---');

echo('o=Abacus.Polynomial([6,12])');
o=Abacus.Polynomial([6,12]);
echo('o.toString()');
echo(o.toString());
echo('o.toTex()');
echo(o.toTex());
o.dispose();
echo('---');

echo('o=Abacus.Polynomial([-4,0,-2,1])');
o=Abacus.Polynomial([-4,0,-2,1]);
echo('o.toString()');
Expand Down Expand Up @@ -156,14 +185,20 @@ echo(Abacus.Polynomial([1,1]).roots().map(function(r){return r.toString();}).joi
echo('Abacus.Polynomial([-1,1,0,2]).roots()'); // no rational roots
echo(Abacus.Polynomial([-1,1,0,2]).roots().map(function(r){return r.toString();}).join(', '));

echo('Abacus.Polynomial(6,-7,0,1]).roots()'); // 1,2,-3
echo('Abacus.Polynomial([6,-7,0,1]).roots()'); // 1,2,-3
echo(Abacus.Polynomial([6,-7,0,1]).roots().map(function(r){return r.toString();}).join(', '));

echo('Abacus.Polynomial(6,-7,0,1]).shift(2).roots()'); // 0,0,1,2,-3
echo('Abacus.Polynomial([6,-7,0,1]).shift(2).roots()'); // 0,0,1,2,-3
echo(Abacus.Polynomial([6,-7,0,1]).shift(2).roots().map(function(r){return r.toString();}).join(', '));

echo('Abacus.Polynomial(-2,5,-5,3]).roots()'); // one root
echo('Abacus.Polynomial([-2,5,-5,3]).roots()'); // one root
echo(Abacus.Polynomial([-2,5,-5,3]).roots().map(function(r){return r.toString();}).join(', '));

echo('Abacus.Polynomial([1,1]).pow(2).roots()'); // multiple root
echo(Abacus.Polynomial([1,1]).pow(2).roots().map(function(r){return r.toString();}).join(', '));

echo('Abacus.Polynomial([1,1]).pow(2).mul(Abacus.Polynomial([0,0,1])).roots()'); // multiple roots
echo(Abacus.Polynomial([1,1]).pow(2).mul(Abacus.Polynomial([0,0,1])).roots().map(function(r){return r.toString();}).join(', '));
echo('---');

echo('Polynomial GCD, generalisation of GCD of numbers');
Expand Down Expand Up @@ -195,6 +230,12 @@ echo('---');

echo('Polynomial Extended GCD, generalisation of xGCD of numbers');
echo('---');
echo('Abacus.Math.polyxgcd(Abacus.Polynomial([2,0,1]),Abacus.Polynomial([6,12]))');
args=[Abacus.Polynomial([2,0,1]),Abacus.Polynomial([6,12])];
o=Abacus.Math.polyxgcd(args);
check_xgcd(args, o);
check_xgcd(args, o, true);

echo('Abacus.Math.polyxgcd(Abacus.Polynomial([1,2]),Abacus.Polynomial([1,3,4]))');
args=[Abacus.Polynomial([1,2]),Abacus.Polynomial([1,3,4])];
o=Abacus.Math.polyxgcd(args);
Expand Down
27 changes: 24 additions & 3 deletions test/polynomials.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,14 @@ o.shift(1)
2*x^2+x
o.shift(-1)
2
o.compose(Abacus.Polynomial([1]))
3
o.compose(Abacus.Polynomial([0,1]))
2*x+1
o.compose(Abacus.Polynomial([1,1]))
2*x+3
Abacus.Polynomial([1,1,1]).compose(o)
4*x^2+6*x+3
o.pow(0)
1
o.pow(1)
Expand All @@ -82,6 +90,12 @@ Abacus.Polynomial.fromExpr(o.toExpr())
2*x+1
o.dispose()
---
o=Abacus.Polynomial([6,12])
o.toString()
12*x+6
o.toTex()
12x+6
---
o=Abacus.Polynomial([-4,0,-2,1])
o.toString()
x^3-2*x^2-4
Expand Down Expand Up @@ -109,12 +123,16 @@ Abacus.Polynomial([1,1]).roots()
-1
Abacus.Polynomial([-1,1,0,2]).roots()

Abacus.Polynomial(6,-7,0,1]).roots()
Abacus.Polynomial([6,-7,0,1]).roots()
1, 2, -3
Abacus.Polynomial(6,-7,0,1]).shift(2).roots()
Abacus.Polynomial([6,-7,0,1]).shift(2).roots()
0, 0, 1, 2, -3
Abacus.Polynomial(-2,5,-5,3]).roots()
Abacus.Polynomial([-2,5,-5,3]).roots()
2/3
Abacus.Polynomial([1,1]).pow(2).roots()
-1, -1
Abacus.Polynomial([1,1]).pow(2).mul(Abacus.Polynomial([0,0,1])).roots()
0, 0, -1, -1
---
Polynomial GCD, generalisation of GCD of numbers
---
Expand All @@ -137,6 +155,9 @@ Abacus.Math.polygcd(Abacus.Polynomial([74]),Abacus.Polynomial([32]),Abacus.Polyn
---
Polynomial Extended GCD, generalisation of xGCD of numbers
---
Abacus.Math.polyxgcd(Abacus.Polynomial([2,0,1]),Abacus.Polynomial([6,12]))
(x^2+2)((4/9)) + (12*x+6)(-(1/27)*x+(1/54)) = 1 true
(x^{2}+2)(\frac{4}{9}) + (12x+6)(-\frac{1}{27}x+\frac{1}{54}) = 1 true
Abacus.Math.polyxgcd(Abacus.Polynomial([1,2]),Abacus.Polynomial([1,3,4]))
(2*x+1)(-4*x-1) + (4*x^2+3*x+1)(2) = 1 true
Abacus.Math.polyxgcd(Abacus.Polynomial([1,1,1,1,5]),Abacus.Polynomial([2,1,3]))
Expand Down

0 comments on commit bb3cd98

Please sign in to comment.