findComplexRoots = { var out = new Array(2); var params = {}; function cmod (a, b) { var c var aa = Math.abs(a) var ab = Math.abs(b) if (aa === 0 && ab === 0) return 0 if (aa >= ab) { c = b / a return aa * Math.sqrt(1 + c * c) } c = a / b return ab * Math.sqrt(1 + c * c) } function cdiv (out, a, b, c, d ) { var e, f if (Math.abs(c) >= Math.abs(d)) { e = d / c f = 1 / (c + d * e) out[0] = (a + b * e) * f out[1] = (b - a * e) * f } else { e = c / d f = 1 / (c * e + d) out[0] = (a * e + b) * f out[1] = (b * e - a) * f } return out } var dzrsub = [ 0, 0.81, 0.5727564927611035, 0, -0.5727564927611035, -0.81, -0.5727564927611035, 0, 0.5727564927611035 ]; var dzisub = [ 0, 0, 0.5727564927611035, 0.81, 0.5727564927611035, 0, -0.5727564927611035, -0.81, -0.5727564927611035 ]; var rsub = [ 0.5, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667 ]; function s0Integrand (theta) { var c, s, a, b, rtmp, itmp; c = Math.cos(theta); s = Math.sin(theta); a = params.z0r + c * params.r; b = params.z0i + s * params.r; // tmp <- f(z) params.f(out, a, b); rtmp = out[0]; itmp = out[1]; // out <- f'(z) out[0] = out[2]; out[1] = out[3]; // out /= f(z) cdiv(out, out[0], out[1], rtmp, itmp); // re(tmp * dz / (2 pi i)) return out[0] * c - out[1] * s; } function smIntegrand (out, theta) { var c, s, a, b, rtmp, itmp; c = Math.cos(theta); s = Math.sin(theta); // out <- f(z) params.f(out, params.z0r + c * params.r, params.z0i + s * params.r); // out <- f'(z) / f(z) cdiv(out, out[2], out[3], out[0], out[1]); // out *= dz (/ (2 pi) but leave this part until outside the loop) rtmp = out[0]; out[0] = out[0] * c - out[1] * s; out[1] = rtmp * s + out[1] * c; // out *= z rtmp = out[0]; out[0] = (out[0] * c - out[1] * s); out[1] = (rtmp * s + out[1] * c); // Recursively multiply by z to compute successive moments z^N: // out[i] = out[i - 1] * z for (var i = 2; i < 2 * params.M; i += 2) { out[i] = out[i - 2] * c - out[i - 1] * s; out[i + 1] = out[i - 2] * s + out[i - 1] * c; } return out } function locateZeros (qr, qi, f, zr, zi, r, z0r, z0i, r0, tol, maxIntegrationDepth, s0IntegrationTolerance, s0IntegrationDepth, curDepth, maxRootsPerCircle, searchCircleOutput) { var p, s0, er, ei, i, j, qri, qii, qmod, qC, M, rr; var qNr, qNi, tmp, qNew, success; // Exit if we've reached the recursion limit: if (curDepth <= 0) { return false; } // Exit early if this circle is completely outside the search region if (cmod(zr - z0r, zi - z0i) - r > r0) { return true } params.z0r = zr; params.z0i = zi; params.r = r; params.f = f; s0 = adaptiveSimpson(s0Integrand, 0, Math.PI * 2, s0IntegrationTolerance, s0IntegrationDepth) * r * (0.5 / Math.PI); M = Math.round(s0); qC = 0; for (i = qr.length - 1; i >= 0; i--) { // Count this zero if it's inside the circle: if (cmod(qr[i] - zr, qi[i] - zi) < r) qC++; } if (M - qC <= 0.9 || s0 < 0.9) { } else if (Math.abs(s0 - M) > 1e-1 || M - qC > maxRootsPerCircle) { // If this happens, then we had a pretty inaccurate integration: for (i = 0, success = true; i < 9; i++) { rr = r * rsub[i]; success = locateZeros(qr, qi, f, zr + r * dzrsub[i], zi + r * dzisub[i], rr, z0r, z0i, r0, tol, maxIntegrationDepth, s0IntegrationTolerance, s0IntegrationDepth, curDepth - 1, maxRootsPerCircle, searchCircleOutput ) && success; } } else { // Compute the zeros: params.M = M - qC; var p = vectorAdaptiveSimpson([], smIntegrand, 0, Math.PI * 2, tol, maxIntegrationDepth, 2 * params.M); // post-multiply the integral by this constant factor var f = r * 0.5 / Math.PI; for (var i = 2 * params.M - 1; i >= 0; i--) { p[i] *= f; } // Deflate the result by known zeros inside the current circle: for (i = qr.length - 1; i >= 0; i--) { // Scale to the unit circle in which we're working qri = (qr[i] - zr) / r; qii = (qi[i] - zi) / r; if (cmod(qri, qii) < 1) { // Initialize the powers of qN var qNr = qri; var qNi = qii; for (j = 0, qNr = 1, qNi = 0; j < p.length; j += 2) { tmp = qNr; qNr = tmp * qri - qNi * qii; qNi = tmp * qii + qNi * qri; // Subtract off the power of the zero: p[j] -= qNr; p[j + 1] -= qNi; } } } // Apply Newton's Identities, constructing a polynomial sharing the roots // of the analytic function: er = []; ei = []; newtonsIdentities(er, ei, params.M, p); // Compute the roots of the constructed polynomial: qNew = findRoots(er, ei, 100 * params.M * params.M, tol); // Append and filter to zeros inside the original contour: for (i = qNew[0].length - 1; i >= 0; i--) { var qrRescaled = zr + qNew[0][i] * r; var qiRescaled = zi + qNew[1][i] * r; //if (cmod(qrRescaled - z0r, qiRescaled - z0i) < r0) { qr.push(qrRescaled); qi.push(qiRescaled); //} } } if (searchCircleOutput) { searchCircleOutput.push({z: [zr, zi], r: r, text: 's₀ = ' + s0.toFixed(5).toString().replace(/(\.[0-9]*[1-9]|\.0)0*$/,'$1').replace(/\.0$/,'').replace(/^-0$/, '0')}) } return [qr, qi]; } return function delvesLyness (f, z0, r, opts) { opts = opts || {} var s0IntegrationDepth = opts.s0IntegrationDepth === undefined ? 13 : opts.s0IntegrationDepth var s0IntegrationTolerance = opts.s0IntegrationTolerance === undefined ? 1e-6 : opts.s0IntegrationTolerance var maxIntegrationDepth = opts.maxIntegrationDepth === undefined ? 12 : opts.maxIntegrationDepth var integrationTolerance = opts.integrationTolerance === undefined ? 1e-8 : opts.integrationTolerance var maxSearchDepth = opts.maxSearchDepth === undefined ? 4 : opts.maxSearchDepth var maxRootsPerCircle = opts.maxRootsPerCircle === undefined ? 4 : opts.maxRootsPerCircle if (opts.searchOutput) opts.searchOutput.length = 0; var foundRoots = locateZeros([], [], f, z0[0], z0[1], r, z0[0], z0[1], r, integrationTolerance, maxIntegrationDepth, s0IntegrationTolerance, s0IntegrationDepth, maxSearchDepth, maxRootsPerCircle, opts.searchOutput) var output = []; for (var i = 0; i < foundRoots[0].length; i++) { if (cmod(foundRoots[0][i] - z0[0], foundRoots[1][i] - z0[1]) <= r) { output.push([foundRoots[0][i], foundRoots[1][i]]) } } return output; } }