Skip to content

Commit bb3cd98

Browse files
foo123foo123
committed
v.0.9.8 contd
* 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
1 parent 0d28129 commit bb3cd98

File tree

4 files changed

+137
-31
lines changed

4 files changed

+137
-31
lines changed

‎src/js/Abacus.js

+63-19
Original file line numberDiff line numberDiff line change
@@ -1730,12 +1730,12 @@ function is_prime( n )
17301730
//n = Arithmetic.abs(/*N(*/n/*)*/);
17311731
ndigits = Arithmetic.digits(n).length;
17321732
// try to use fastest algorithm based on size of number (number of digits)
1733-
if ( ndigits <= 8 )
1733+
if ( ndigits <= 6 )
17341734
{
17351735
// deterministic test
17361736
return wheel_trial_div_test(n);
17371737
}
1738-
else if ( ndigits <= 1000 )
1738+
else if ( ndigits <= 20 )
17391739
{
17401740
// deterministic test
17411741
/*
@@ -1746,11 +1746,11 @@ function is_prime( n )
17461746
If n < 3474749660383 is a 2, 3, 5, 7, 11 and 13-SPRP, then n is prime [Jaeschke93].
17471747
If n < 341550071728321 is a 2, 3, 5, 7, 11, 13 and 17-SPRP, then n is prime [Jaeschke93].
17481748
*/
1749-
/*if ( Arithmetic.lt(n, N(1373653)) )
1749+
if ( Arithmetic.lt(n, N(1373653)) )
17501750
return miller_rabin_test(n, [two, N(3)]);
17511751
else if ( Arithmetic.lt(n, N("25326001")) )
17521752
return miller_rabin_test(n, [two, N(3), N(5)]);
1753-
else*/ if ( Arithmetic.lt(n, N("25000000000")) )
1753+
else if ( Arithmetic.lt(n, N("25000000000")) )
17541754
return Arithmetic.equ(n, N("3215031751")) ? false : miller_rabin_test(n, [two, N(3), N(5), N(7)]);
17551755
else if ( Arithmetic.lt(n, N("2152302898747")) )
17561756
return miller_rabin_test(n, [two, N(3), N(5), N(7), N(11)]);
@@ -2704,14 +2704,13 @@ function polyxgcd( /* args */ )
27042704
}
27052705
function divisors( n, as_generator )
27062706
{
2707-
// time+space O(sqrt(n)) to find all distinct divisors of n (including 1 and n itself)
27082707
var Arithmetic = Abacus.Arithmetic, O = Arithmetic.O, I = Arithmetic.I,
27092708
list = null, D2 = null, D1 = null, L1 = 0, L2 = 0, node, sqrn, i, n_i, next, factors;
27102709
//n = Arithmetic.num(n);
27112710
n = Arithmetic.abs(n);
27122711
if ( true===as_generator )
27132712
{
2714-
if ( Arithmetic.gte(n, 1000000) )
2713+
if ( Arithmetic.gte(n, 100000) )
27152714
{
27162715
// for very large numbers,
27172716
// compute divisors through prime factorisation
@@ -2727,6 +2726,7 @@ function divisors( n, as_generator )
27272726
}
27282727
else
27292728
{
2729+
// time+space O(sqrt(n)) to find all distinct divisors of n (including 1 and n itself)
27302730
sqrn = isqrt(n);
27312731
i = I; next = null;
27322732
// return iterator/generator
@@ -2766,6 +2766,7 @@ function divisors( n, as_generator )
27662766
}
27672767
else
27682768
{
2769+
// time+space O(sqrt(n)) to find all distinct divisors of n (including 1 and n itself)
27692770
sqrn = isqrt(n);
27702771
for (i=I; Arithmetic.lte(i,sqrn); i=Arithmetic.add(i,I))
27712772
{
@@ -3549,8 +3550,8 @@ function factorial( n, m )
35493550
// simple factorial = F(n) = n F(n-1) = n!
35503551
if ( 12 >= n ) return 0 > n ? O : NUM(([1,1,2,6,24,120,720,5040,40320,362880,3628800,39916800,479001600 /*MAX: 2147483647*/])[n]);
35513552

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

35553556
key = String(n)/*+'!'*/;
35563557
if ( null == factorial.mem1[key] )
@@ -7101,11 +7102,11 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
71017102
,primitive: function( and_content ) {
71027103
// factorise into content and primitive part
71037104
// https://en.wikipedia.org/wiki/Factorization_of_polynomials#Primitive_part%E2%80%93content_factorization
7104-
var self = this, Arithmetic = Abacus.Arithmetic, coeff = self.coeff, coeffp, M, content;
7105+
var self = this, Arithmetic = Abacus.Arithmetic, coeff = self.coeff, coeffp, LCM, content;
71057106
if ( null == self._prim )
71067107
{
7107-
M = coeff.reduce(function(M, c){return Arithmetic.mul(M, c.den);}, Arithmetic.I);
7108-
coeffp = coeff.map(function(c){return c.mul(M).num;});
7108+
LCM = lcm(coeff.map(function(c){return c.den;}));
7109+
coeffp = coeff.map(function(c){return c.mul(LCM).integer();});
71097110
content = gcd(coeffp);
71107111
coeffp = coeffp.map(function(c){return Arithmetic.div(c, content);});
71117112
// make positive lead
@@ -7114,7 +7115,7 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
71147115
coeffp = coeffp.map(function(c){return Arithmetic.neg(c);});
71157116
content = Arithmetic.neg(content);
71167117
}
7117-
self._prim = [Polynomial(coeffp, self.symbol), Rational(content, M)];
7118+
self._prim = [Polynomial(coeffp, self.symbol), Rational(content, LCM)];
71187119
}
71197120
return true===and_content ? self._prim.slice() : self._prim[0];
71207121
}
@@ -7123,7 +7124,7 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
71237124
// https://en.wikipedia.org/wiki/Rational_root_theorem
71247125
// https://en.wikipedia.org/wiki/Gauss%27s_lemma_(polynomial)
71257126
var self = this, Arithmetic = Abacus.Arithmetic, O = Arithmetic.O,
7126-
roots = [], primitive,c, i, d0, dn, iter, comb, root;
7127+
roots = [], primitive,c, p, i, d0, dn, iter, comb, root, nroot, found;
71277128

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

@@ -7137,17 +7138,37 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
71377138
if ( i+1<c.length )
71387139
{
71397140
// try all possible rational divisors of c_0(excluding trivial zero terms) and c_n
7140-
d0 = divisors(c[i].num); dn = divisors(c[c.length-1].num);
7141+
iter = divisors(c[i].num, true);
7142+
d0 = iter.get(); iter.dispose();
7143+
iter = divisors(c[c.length-1].num, true);
7144+
dn = iter.get(); iter.dispose();
7145+
71417146
iter = Tensor([d0.length, dn.length]);
71427147
while(iter.hasNext())
71437148
{
71447149
comb = iter.next();
7145-
// try positive root
7150+
// positive root
71467151
root = Rational(d0[comb[0]], dn[comb[1]]);
7147-
if ( primitive.valueOf(root).equ(O) ) roots.push(root);
7148-
// try negative root
7149-
root = root.neg();
7150-
if ( primitive.valueOf(root).equ(O) ) roots.push(root);
7152+
// negative root
7153+
nroot = Rational(Arithmetic.neg(d0[comb[0]]), dn[comb[1]]);
7154+
p = primitive; found = true;
7155+
while( found && (0<p.deg()) )
7156+
{
7157+
found = false;
7158+
// try positive root
7159+
if ( p.valueOf(root).equ(O) )
7160+
{
7161+
roots.push(root);
7162+
found = true;
7163+
}
7164+
// try negative root
7165+
if ( p.valueOf(nroot).equ(O) )
7166+
{
7167+
roots.push(nroot);
7168+
found = true;
7169+
}
7170+
if ( found ) p = p.d(); // get derivative to check if roots are multiple
7171+
}
71517172
}
71527173
iter.dispose();
71537174
}
@@ -7295,6 +7316,29 @@ Polynomial = Abacus.Polynomial = Class(INumber, {
72957316
return pow;
72967317
}
72977318
}
7319+
,compose: function( x ) {
7320+
// functionaly compose one polynomial with another. ie result = P(Q(x))
7321+
var self = this, Arithmetic = Abacus.Arithmetic, O = Arithmetic.O, px, c, i;
7322+
if ( x instanceof Matrix ) return null;
7323+
if ( x instanceof Term ) x = Expr(x);
7324+
if ( x instanceof Expr ) x = Polynomial.fromExpr(x, self.symbol);
7325+
if ( x instanceof Complex ) x = x.real;
7326+
if ( Arithmetic.isNumber(x) || (x instanceof Rational) )
7327+
{
7328+
return Polynomial([self.valueOf(x)], self.symbol);
7329+
}
7330+
else if ( x instanceof Polynomial )
7331+
{
7332+
// Composition through variation of Horner's algorithm for fast evaluation
7333+
if ( 0 === self.deg() ) return self.clone();
7334+
if ( 0 === x.deg() ) return Polynomial([self.valueOf(x.c())], x.symbol);
7335+
c = self.coeff; i = c.length-1;
7336+
px = Polynomial([c[i]], x.symbol);
7337+
while(i--) px = Polynomial([c[i]], x.symbol).add(px.mul(x));
7338+
return px;
7339+
}
7340+
return self;
7341+
}
72987342
,shift: function( s ) {
72997343
// shift <-> equivalent to multiplication/division by a monomial x^s
73007344
var self = this, Arithmetic = Abacus.Arithmetic;

‎src/js/Abacus.min.js

+1-1
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

‎test/polynomials.js

+49-8
Original file line numberDiff line numberDiff line change
@@ -10,15 +10,27 @@ function check_div( n, d, q, r )
1010
var nn = q.mul(d).add(r);
1111
console.log('('+n.toString()+')/('+d.toString()+')=('+d.toString()+')*('+q.toString()+')+('+r.toString()+')='+nn.toString(), nn.equ(n));
1212
}
13-
function check_xgcd( args, gcd )
13+
function check_xgcd( args, gcd, tex )
1414
{
1515
var out = '', res = Abacus.Polynomial(0);
16-
for(i=0; i<args.length; i++)
16+
if ( tex )
1717
{
18-
out += (out.length ? ' + ' : '') + '('+args[i].toString()+')'+'('+gcd[i+1].toString()+')';
19-
res = res.add(args[i].mul(gcd[i+1]));
18+
for(i=0; i<args.length; i++)
19+
{
20+
out += (out.length ? ' + ' : '') + '('+args[i].toTex()+')'+'('+gcd[i+1].toTex()+')';
21+
res = res.add(args[i].mul(gcd[i+1]));
22+
}
23+
out += ' = '+gcd[0].toTex();
24+
}
25+
else
26+
{
27+
for(i=0; i<args.length; i++)
28+
{
29+
out += (out.length ? ' + ' : '') + '('+args[i].toString()+')'+'('+gcd[i+1].toString()+')';
30+
res = res.add(args[i].mul(gcd[i+1]));
31+
}
32+
out += ' = '+gcd[0].toString();
2033
}
21-
out += ' = '+gcd[0].toString();
2234
console.log(out, res.equ(gcd[0]));
2335
}
2436
var o, d, qr, args;
@@ -87,6 +99,14 @@ echo('o.shift(1)');
8799
echo(o.shift(1).toString());
88100
echo('o.shift(-1)');
89101
echo(o.shift(-1).toString());
102+
echo('o.compose(Abacus.Polynomial([1]))');
103+
echo(o.compose(Abacus.Polynomial([1])).toString());
104+
echo('o.compose(Abacus.Polynomial([0,1]))');
105+
echo(o.compose(Abacus.Polynomial([0,1])).toString());
106+
echo('o.compose(Abacus.Polynomial([1,1]))');
107+
echo(o.compose(Abacus.Polynomial([1,1])).toString());
108+
echo('Abacus.Polynomial([1,1,1]).compose(o)');
109+
echo(Abacus.Polynomial([1,1,1]).compose(o).toString());
90110
echo('o.pow(0)');
91111
echo(o.pow(0).toString());
92112
echo('o.pow(1)');
@@ -117,6 +137,15 @@ echo('o.dispose()');
117137
o.dispose();
118138
echo('---');
119139

140+
echo('o=Abacus.Polynomial([6,12])');
141+
o=Abacus.Polynomial([6,12]);
142+
echo('o.toString()');
143+
echo(o.toString());
144+
echo('o.toTex()');
145+
echo(o.toTex());
146+
o.dispose();
147+
echo('---');
148+
120149
echo('o=Abacus.Polynomial([-4,0,-2,1])');
121150
o=Abacus.Polynomial([-4,0,-2,1]);
122151
echo('o.toString()');
@@ -156,14 +185,20 @@ echo(Abacus.Polynomial([1,1]).roots().map(function(r){return r.toString();}).joi
156185
echo('Abacus.Polynomial([-1,1,0,2]).roots()'); // no rational roots
157186
echo(Abacus.Polynomial([-1,1,0,2]).roots().map(function(r){return r.toString();}).join(', '));
158187

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

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

165-
echo('Abacus.Polynomial(-2,5,-5,3]).roots()'); // one root
194+
echo('Abacus.Polynomial([-2,5,-5,3]).roots()'); // one root
166195
echo(Abacus.Polynomial([-2,5,-5,3]).roots().map(function(r){return r.toString();}).join(', '));
196+
197+
echo('Abacus.Polynomial([1,1]).pow(2).roots()'); // multiple root
198+
echo(Abacus.Polynomial([1,1]).pow(2).roots().map(function(r){return r.toString();}).join(', '));
199+
200+
echo('Abacus.Polynomial([1,1]).pow(2).mul(Abacus.Polynomial([0,0,1])).roots()'); // multiple roots
201+
echo(Abacus.Polynomial([1,1]).pow(2).mul(Abacus.Polynomial([0,0,1])).roots().map(function(r){return r.toString();}).join(', '));
167202
echo('---');
168203

169204
echo('Polynomial GCD, generalisation of GCD of numbers');
@@ -195,6 +230,12 @@ echo('---');
195230

196231
echo('Polynomial Extended GCD, generalisation of xGCD of numbers');
197232
echo('---');
233+
echo('Abacus.Math.polyxgcd(Abacus.Polynomial([2,0,1]),Abacus.Polynomial([6,12]))');
234+
args=[Abacus.Polynomial([2,0,1]),Abacus.Polynomial([6,12])];
235+
o=Abacus.Math.polyxgcd(args);
236+
check_xgcd(args, o);
237+
check_xgcd(args, o, true);
238+
198239
echo('Abacus.Math.polyxgcd(Abacus.Polynomial([1,2]),Abacus.Polynomial([1,3,4]))');
199240
args=[Abacus.Polynomial([1,2]),Abacus.Polynomial([1,3,4])];
200241
o=Abacus.Math.polyxgcd(args);

‎test/polynomials.txt

+24-3
Original file line numberDiff line numberDiff line change
@@ -60,6 +60,14 @@ o.shift(1)
6060
2*x^2+x
6161
o.shift(-1)
6262
2
63+
o.compose(Abacus.Polynomial([1]))
64+
3
65+
o.compose(Abacus.Polynomial([0,1]))
66+
2*x+1
67+
o.compose(Abacus.Polynomial([1,1]))
68+
2*x+3
69+
Abacus.Polynomial([1,1,1]).compose(o)
70+
4*x^2+6*x+3
6371
o.pow(0)
6472
1
6573
o.pow(1)
@@ -82,6 +90,12 @@ Abacus.Polynomial.fromExpr(o.toExpr())
8290
2*x+1
8391
o.dispose()
8492
---
93+
o=Abacus.Polynomial([6,12])
94+
o.toString()
95+
12*x+6
96+
o.toTex()
97+
12x+6
98+
---
8599
o=Abacus.Polynomial([-4,0,-2,1])
86100
o.toString()
87101
x^3-2*x^2-4
@@ -109,12 +123,16 @@ Abacus.Polynomial([1,1]).roots()
109123
-1
110124
Abacus.Polynomial([-1,1,0,2]).roots()
111125

112-
Abacus.Polynomial(6,-7,0,1]).roots()
126+
Abacus.Polynomial([6,-7,0,1]).roots()
113127
1, 2, -3
114-
Abacus.Polynomial(6,-7,0,1]).shift(2).roots()
128+
Abacus.Polynomial([6,-7,0,1]).shift(2).roots()
115129
0, 0, 1, 2, -3
116-
Abacus.Polynomial(-2,5,-5,3]).roots()
130+
Abacus.Polynomial([-2,5,-5,3]).roots()
117131
2/3
132+
Abacus.Polynomial([1,1]).pow(2).roots()
133+
-1, -1
134+
Abacus.Polynomial([1,1]).pow(2).mul(Abacus.Polynomial([0,0,1])).roots()
135+
0, 0, -1, -1
118136
---
119137
Polynomial GCD, generalisation of GCD of numbers
120138
---
@@ -137,6 +155,9 @@ Abacus.Math.polygcd(Abacus.Polynomial([74]),Abacus.Polynomial([32]),Abacus.Polyn
137155
---
138156
Polynomial Extended GCD, generalisation of xGCD of numbers
139157
---
158+
Abacus.Math.polyxgcd(Abacus.Polynomial([2,0,1]),Abacus.Polynomial([6,12]))
159+
(x^2+2)((4/9)) + (12*x+6)(-(1/27)*x+(1/54)) = 1 true
160+
(x^{2}+2)(\frac{4}{9}) + (12x+6)(-\frac{1}{27}x+\frac{1}{54}) = 1 true
140161
Abacus.Math.polyxgcd(Abacus.Polynomial([1,2]),Abacus.Polynomial([1,3,4]))
141162
(2*x+1)(-4*x-1) + (4*x^2+3*x+1)(2) = 1 true
142163
Abacus.Math.polyxgcd(Abacus.Polynomial([1,1,1,1,5]),Abacus.Polynomial([2,1,3]))

0 commit comments

Comments
 (0)