Skip to content

Commit b21831d

Browse files
pboyerpboyer
committed
Add cubic deflation. Don't use Math.pow in Polynomial.evaluate. 16x faster.
1 parent c95076e commit b21831d

File tree

6 files changed

+90
-22
lines changed

6 files changed

+90
-22
lines changed

‎README.md

+8-1
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,13 @@ const epsilon = 1e-6;
4141
const roots = polynomialRoots(f, startSearchInterval, endSearchInterval, epsilon);
4242
```
4343

44-
## Paper
44+
## Resources
45+
46+
This is a great video intro: https://youtu.be/CoGo_3C7xR0?t=7092
4547

4648
Cem Yuksel. 2022. High-Performance Polynomial Root Finding for Graphics. Proc. ACM Comput. Graph. Interact. Tech. 5, 3, Article 7 (July 2022), 15 pages.
49+
50+
## Changelog
51+
52+
- 1.0.11 - Add cubic deflation. Don't use Math.pow in Polynomial.evaluate. 16x faster!
53+
- 1.0.10 - First public release

‎dist/nomial.js

+25-6
Original file line numberDiff line numberDiff line change
@@ -16,9 +16,11 @@ class Polynomial {
1616
this.degree = coefficients.length - 1;
1717
}
1818
evaluate(x) {
19+
let xx = 1;
1920
let result = 0;
2021
for (let i = 0, l = this.coefficients.length; i < l; i++) {
21-
result += this.coefficients[i] * Math.pow(x, i);
22+
result += this.coefficients[i] * xx;
23+
xx *= x;
2224
}
2325
return result;
2426
}
@@ -48,16 +50,33 @@ function polynomialRoots(f, startInterval = -1000, endInterval = 1000, epsilon =
4850
const rootsOfDerivative = polynomialRoots(derivative, startInterval, endInterval, epsilon);
4951
rootsOfDerivative.push(endInterval);
5052
let a = startInterval;
53+
let fa = f.evaluate(startInterval);
5154
const roots = [];
5255
for (let i = 0; i < rootsOfDerivative.length; i++) {
53-
let b = rootsOfDerivative[i];
54-
if (Math.sign(f.evaluate(a)) !== Math.sign(f.evaluate(b))) {
55-
roots.push(findRoot(f, derivative, a, b, epsilon));
56+
const b = rootsOfDerivative[i];
57+
const fb = f.evaluate(b);
58+
if (Math.sign(fa) !== Math.sign(fb)) {
59+
const r = findRoot(f, derivative, a, b, fa, epsilon);
60+
if (f.degree === 3) {
61+
const deflated = deflate(f, r);
62+
return [r, ...findQuadraticRoots(deflated, startInterval, endInterval)];
63+
}
64+
roots.push(r);
5665
}
5766
a = b;
67+
fa = fb;
5868
}
5969
return roots;
6070
}
71+
function deflate(f, xr) {
72+
const c = f.coefficients[1];
73+
const b = f.coefficients[2];
74+
const a = f.coefficients[3];
75+
const ap = a;
76+
const bp = b + ap * xr;
77+
const cp = c + bp * xr;
78+
return new Polynomial([cp, bp, ap]);
79+
}
6180
function findQuadraticRoots(f, startInterval, endInterval) {
6281
const c = f.coefficients[0];
6382
const b = f.coefficients[1];
@@ -85,7 +104,7 @@ function multSign(v, sign) {
85104
return v * (sign < 0 ? -1 : 1);
86105
}
87106
// Following http://www.cemyuksel.com/research/polynomials/polynomial_roots_hpg2022_supplemental.pdf
88-
function findRoot(f, deriv, x1, x2, epsilon) {
107+
function findRoot(f, deriv, x1, x2, fx1, epsilon) {
89108
let xr = (x1 + x2) / 2;
90109
if (Math.abs(x2 - x1) <= 2 * epsilon) {
91110
return xr;
@@ -104,7 +123,7 @@ function findRoot(f, deriv, x1, x2, epsilon) {
104123
xr = (x1 + x2) / 2;
105124
}
106125
}
107-
const y1 = f.evaluate(x1);
126+
const y1 = fx1;
108127
let yr = f.evaluate(xr);
109128
while (true) {
110129
if (Math.sign(yr) === Math.sign(y1)) {

‎dist/nomial.min.js

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

‎examples/index.html

+2-2
Original file line numberDiff line numberDiff line change
@@ -377,12 +377,12 @@
377377
for (const pt of pts) {
378378
ctx.fillStyle = "white";
379379
ctx.beginPath();
380-
ctx.arc(pt.x, pt.y, 3 * ptradius, 0, 2 * Math.PI);
380+
ctx.arc(pt.x, pt.y, 2.5 * ptradius, 0, 2 * Math.PI);
381381
ctx.fill();
382382

383383
ctx.fillStyle = "black";
384384
ctx.beginPath();
385-
ctx.arc(pt.x, pt.y, 2 * ptradius, 0, 2 * Math.PI);
385+
ctx.arc(pt.x, pt.y, 1.5 * ptradius, 0, 2 * Math.PI);
386386
ctx.fill();
387387
}
388388
}

‎examples/nomial.js

+25-6
Original file line numberDiff line numberDiff line change
@@ -16,9 +16,11 @@ class Polynomial {
1616
this.degree = coefficients.length - 1;
1717
}
1818
evaluate(x) {
19+
let xx = 1;
1920
let result = 0;
2021
for (let i = 0, l = this.coefficients.length; i < l; i++) {
21-
result += this.coefficients[i] * Math.pow(x, i);
22+
result += this.coefficients[i] * xx;
23+
xx *= x;
2224
}
2325
return result;
2426
}
@@ -48,16 +50,33 @@ function polynomialRoots(f, startInterval = -1000, endInterval = 1000, epsilon =
4850
const rootsOfDerivative = polynomialRoots(derivative, startInterval, endInterval, epsilon);
4951
rootsOfDerivative.push(endInterval);
5052
let a = startInterval;
53+
let fa = f.evaluate(startInterval);
5154
const roots = [];
5255
for (let i = 0; i < rootsOfDerivative.length; i++) {
53-
let b = rootsOfDerivative[i];
54-
if (Math.sign(f.evaluate(a)) !== Math.sign(f.evaluate(b))) {
55-
roots.push(findRoot(f, derivative, a, b, epsilon));
56+
const b = rootsOfDerivative[i];
57+
const fb = f.evaluate(b);
58+
if (Math.sign(fa) !== Math.sign(fb)) {
59+
const r = findRoot(f, derivative, a, b, fa, epsilon);
60+
if (f.degree === 3) {
61+
const deflated = deflate(f, r);
62+
return [r, ...findQuadraticRoots(deflated, startInterval, endInterval)];
63+
}
64+
roots.push(r);
5665
}
5766
a = b;
67+
fa = fb;
5868
}
5969
return roots;
6070
}
71+
function deflate(f, xr) {
72+
const c = f.coefficients[1];
73+
const b = f.coefficients[2];
74+
const a = f.coefficients[3];
75+
const ap = a;
76+
const bp = b + ap * xr;
77+
const cp = c + bp * xr;
78+
return new Polynomial([cp, bp, ap]);
79+
}
6180
function findQuadraticRoots(f, startInterval, endInterval) {
6281
const c = f.coefficients[0];
6382
const b = f.coefficients[1];
@@ -85,7 +104,7 @@ function multSign(v, sign) {
85104
return v * (sign < 0 ? -1 : 1);
86105
}
87106
// Following http://www.cemyuksel.com/research/polynomials/polynomial_roots_hpg2022_supplemental.pdf
88-
function findRoot(f, deriv, x1, x2, epsilon) {
107+
function findRoot(f, deriv, x1, x2, fx1, epsilon) {
89108
let xr = (x1 + x2) / 2;
90109
if (Math.abs(x2 - x1) <= 2 * epsilon) {
91110
return xr;
@@ -104,7 +123,7 @@ function findRoot(f, deriv, x1, x2, epsilon) {
104123
xr = (x1 + x2) / 2;
105124
}
106125
}
107-
const y1 = f.evaluate(x1);
126+
const y1 = fx1;
108127
let yr = f.evaluate(xr);
109128
while (true) {
110129
if (Math.sign(yr) === Math.sign(y1)) {

‎src/index.ts

+29-6
Original file line numberDiff line numberDiff line change
@@ -12,9 +12,11 @@ export class Polynomial {
1212
}
1313

1414
evaluate(x: number): number {
15+
let xx = 1;
1516
let result = 0;
1617
for (let i = 0, l = this.coefficients.length; i < l; i++) {
17-
result += this.coefficients[i] * Math.pow(x, i);
18+
result += this.coefficients[i] * xx;
19+
xx *= x;
1820
}
1921
return result;
2022
}
@@ -51,20 +53,41 @@ export function polynomialRoots(f: Polynomial, startInterval: number = -1000, en
5153
rootsOfDerivative.push(endInterval);
5254

5355
let a = startInterval;
56+
let fa = f.evaluate(startInterval);
5457

5558
const roots: number[] = [];
5659

5760
for (let i = 0; i < rootsOfDerivative.length; i++) {
58-
let b = rootsOfDerivative[i];
59-
if (Math.sign(f.evaluate(a)) !== Math.sign(f.evaluate(b))) {
60-
roots.push(findRoot(f, derivative, a, b, epsilon));
61+
const b = rootsOfDerivative[i];
62+
const fb = f.evaluate(b);
63+
64+
if (Math.sign(fa) !== Math.sign(fb)) {
65+
const r = findRoot(f, derivative, a, b, fa, epsilon);
66+
if (f.degree === 3) {
67+
const deflated = deflate(f, r);
68+
return [r, ...findQuadraticRoots(deflated, startInterval, endInterval)];
69+
}
70+
71+
roots.push(r);
6172
}
73+
6274
a = b;
75+
fa = fb;
6376
}
6477

6578
return roots;
6679
}
6780

81+
function deflate(f: Polynomial, xr: number) {
82+
const c = f.coefficients[1];
83+
const b = f.coefficients[2];
84+
const a = f.coefficients[3];
85+
const ap = a;
86+
const bp = b + ap * xr;
87+
const cp = c + bp * xr;
88+
return new Polynomial([cp, bp, ap]);
89+
}
90+
6891
function findQuadraticRoots(f: Polynomial, startInterval: number, endInterval: number) {
6992
const c = f.coefficients[0];
7093
const b = f.coefficients[1];
@@ -100,7 +123,7 @@ function multSign(v: number, sign: number) {
100123
}
101124

102125
// Following http://www.cemyuksel.com/research/polynomials/polynomial_roots_hpg2022_supplemental.pdf
103-
function findRoot(f: Polynomial, deriv: Polynomial, x1: number, x2: number, epsilon: number): number {
126+
function findRoot(f: Polynomial, deriv: Polynomial, x1: number, x2: number, fx1: number, epsilon: number): number {
104127
let xr = (x1 + x2) / 2;
105128

106129
if (Math.abs(x2 - x1) <= 2 * epsilon) {
@@ -125,7 +148,7 @@ function findRoot(f: Polynomial, deriv: Polynomial, x1: number, x2: number, epsi
125148
}
126149
}
127150

128-
const y1 = f.evaluate(x1);
151+
const y1 = fx1;
129152

130153
let yr = f.evaluate(xr);
131154

0 commit comments

Comments
 (0)