Last active
September 23, 2025 09:06
-
-
Save kvedala/9713842d5e264c94e293a930b479f597 to your computer and use it in GitHub Desktop.
Durand-Kerner Roots - Python
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "metadata": { | |
| "toc": true | |
| }, | |
| "cell_type": "markdown", | |
| "source": "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Fast-Numerical-algorithm-to-find-all-the-roots-of-an-n-th-degree-polynomial\" data-toc-modified-id=\"Fast-Numerical-algorithm-to-find-all-the-roots-of-an-n-th-degree-polynomial-1\"><span class=\"toc-item-num\">1 </span>Fast-Numerical algorithm to find all the roots of an n-th degree polynomial</a></span></li><li><span><a href=\"#Using-Numba\" data-toc-modified-id=\"Using-Numba-2\"><span class=\"toc-item-num\">2 </span>Using Numba</a></span></li><li><span><a href=\"#Using-Cython\" data-toc-modified-id=\"Using-Cython-3\"><span class=\"toc-item-num\">3 </span>Using Cython</a></span></li><li><span><a href=\"#Testing-&-Benchmarking\" data-toc-modified-id=\"Testing-&-Benchmarking-4\"><span class=\"toc-item-num\">4 </span>Testing & Benchmarking</a></span><ul class=\"toc-item\"><li><span><a href=\"#Known-bugs-in-Durand–Kerner-algorithm\" data-toc-modified-id=\"Known-bugs-in-Durand–Kerner-algorithm-4.1\"><span class=\"toc-item-num\">4.1 </span>Known bugs in Durand–Kerner algorithm</a></span><ul class=\"toc-item\"><li><span><a href=\"#Fail-cases\" data-toc-modified-id=\"Fail-cases-4.1.1\"><span class=\"toc-item-num\">4.1.1 </span>Fail cases</a></span></li><li><span><a href=\"#Fail-cases-resolved\" data-toc-modified-id=\"Fail-cases-resolved-4.1.2\"><span class=\"toc-item-num\">4.1.2 </span>Fail cases resolved</a></span></li></ul></li></ul></li></ul></div>" | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "# Fast-Numerical algorithm to find all the roots of an n-th degree polynomial\nThe algorithm is an implementation of [Durand–Kerner method](https://en.wikipedia.org/wiki/Durand–Kerner_method). The algorithm is implemented using the [numba](https://numba.pydata.org) optimizer and also using [cython](http://cython.readthedocs.io/). The two implementations are then compared with the implementation of numpy and we see that both usually perform better than `numpy.solve`. " | |
| }, | |
| { | |
| "metadata": { | |
| "ExecuteTime": { | |
| "end_time": "2020-02-29T22:11:53.086454Z", | |
| "start_time": "2020-02-29T22:11:50.238096Z" | |
| }, | |
| "run_control": { | |
| "marked": false | |
| }, | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "import numpy as np\nimport numba\nfrom typing import Callable, List\nfrom math import sqrt\n%reload_ext cython", | |
| "execution_count": 1, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "# Using Numba" | |
| }, | |
| { | |
| "metadata": { | |
| "ExecuteTime": { | |
| "end_time": "2020-02-29T22:11:53.120963Z", | |
| "start_time": "2020-02-29T22:11:53.091450Z" | |
| }, | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "@numba.njit(nogil=True, boundscheck=False, cache=True)\ndef f(poly, x: complex) -> complex:\n out = 0 + 0j\n L = len(poly)\n for i in range(L):\n k = L - i - 1\n if poly[k] != 0:\n out += poly[k] * pow(x, i)\n\n# if poly[i] != 0:\n# out += poly[i] * (x ** i)\n return out\n\n\[email protected](parallel=True, \n nopython=True, boundscheck=False, nogil=True, cache=True)\ndef solve(poly: np.ndarray, max_iter: int = 50000):\n epsilon = 5e-16\n ii = 0\n scale = 500\n L = poly.shape[0] - 1\n roots = (np.random.random(L) + np.random.random(L) * 1j) - (0.5 + 0.5j)\n roots *= scale\n max_delta = 1\n\n while max_delta > epsilon and ii < max_iter:\n max_delta = 0\n ii += 1\n if L > 20:\n for i in numba.prange(L): # for ith root\n deno = 1 + 0j\n delta = f(poly, roots[i])\n for k in range(L):\n if k != i:\n delta /= (roots[i] - roots[k])\n roots[i] -= delta\n # print(ii, roots, delta)\n max_delta += ((delta.real * delta.real) +\n (delta.imag * delta.imag)) / L\n else:\n for i in range(L): # for ith root\n deno = 1 + 0j\n delta = f(poly, roots[i])\n for k in range(L):\n if k != i:\n delta /= (roots[i] - roots[k])\n roots[i] -= delta\n # print(ii, roots, delta)\n max_delta += ((delta.real * delta.real) +\n (delta.imag * delta.imag)) / L\n\n\n max_delta = np.sqrt(max_delta)\n return roots, max_delta\n\n\nsolve(np.array([1, 0, -1]))", | |
| "execution_count": 4, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "execution_count": 4, | |
| "data": { | |
| "text/plain": "(array([-1.+0.j, 1.+0.j]), 1.123017639172393e-22)" | |
| }, | |
| "metadata": {} | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "# Using Cython" | |
| }, | |
| { | |
| "metadata": { | |
| "ExecuteTime": { | |
| "end_time": "2020-02-29T22:12:11.454576Z", | |
| "start_time": "2020-02-29T22:11:53.126778Z" | |
| }, | |
| "scrolled": false, | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "%%cython -a -c=-Ofast -c=-flto=thin\n# cython: boundscheck=False, wraparound=False, cdivision=True\n\ncimport numpy as np\nimport numpy as np\n\ncimport cython\nfrom cython.parallel cimport prange\n\n# from libc.math cimport abs, sqrt, fmax\n# from libcpp.complex cimport pow as cpow, real, imag, abs\n\ncdef extern from \"complex.h\" nogil:\n np.complex128_t cpow(np.complex128_t, long)\n double cabs(np.complex128_t) \n double creal(np.complex128_t) \n double cimag(np.complex128_t) \n\ncdef extern from \"math.h\" nogil:\n long double sqrtl(long double)\n long double fabsl(long double)\n long double fmaxl(long double, long double)\n\n \ncdef np.complex128_t f_c(double[:] poly, np.complex128_t x) nogil:\n cdef:\n np.complex128_t out = 0 + 0j\n long i, k, L = poly.shape[0]\n \n for i in range(L):\n k = L - i - 1\n if poly[k] != 0:\n out += poly[k] * cpow(x, i)\n# if poly[i] != 0:\n# out += <np.complex128_t> poly[i] * (x ** i)\n return out \n\n\ndef solve_c(double[:] poly, long max_iter = 500000):\n cdef:\n double epsilon = 5e-10\n long ii = 0, i, k, L = poly.shape[0] - 1\n double scale = 500.0\n np.ndarray[np.complex128_t, ndim=1] roots = ((np.random.random(L) + np.random.random(L) * 1j) \\\n - (.5 + .5j)) * scale\n# np.complex128_t[:] roots = _roots\n np.complex128_t deno, delta\n long double max_delta = 1, d2, rd, id\n \n with nogil:\n while max_delta > epsilon and ii < max_iter:\n max_delta = 0\n ii += 1\n for i in range(L):\n# deno = 1 + 0j\n# delta = 0 + 0j\n delta = f_c(poly, roots[i])\n for k in range(L):\n if k != i:\n delta /= (roots[i] - roots[k])\n roots[i] -= delta\n rd = creal(delta)\n id = cimag(delta)\n max_delta += ((rd * rd) + (id * id)) / L\n# with gil:\n# print(ii, roots, delta)\n max_delta = sqrtl(max_delta)\n# print(\"max_delta = \", max_delta)\n return roots, max_delta", | |
| "execution_count": 17, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "execution_count": 17, | |
| "data": { | |
| "text/plain": "<IPython.core.display.HTML object>", | |
| "text/html": "<!DOCTYPE html>\n<!-- Generated by Cython 0.29.15 -->\n<html>\n<head>\n <meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n <title>Cython: _cython_magic_7d2d8e83eb391e7da013a05885fe9a7f.pyx</title>\n <style type=\"text/css\">\n \nbody.cython { font-family: courier; font-size: 12; }\n\n.cython.tag { }\n.cython.line { margin: 0em }\n.cython.code { font-size: 9; color: #444444; display: none; margin: 0px 0px 0px 8px; border-left: 8px none; }\n\n.cython.line .run { background-color: #B0FFB0; }\n.cython.line .mis { background-color: #FFB0B0; }\n.cython.code.run { border-left: 8px solid #B0FFB0; }\n.cython.code.mis { border-left: 8px solid #FFB0B0; }\n\n.cython.code .py_c_api { color: red; }\n.cython.code .py_macro_api { color: #FF7000; }\n.cython.code .pyx_c_api { color: #FF3000; }\n.cython.code .pyx_macro_api { color: #FF7000; }\n.cython.code .refnanny { color: #FFA000; }\n.cython.code .trace { color: #FFA000; }\n.cython.code .error_goto { color: #FFA000; }\n\n.cython.code .coerce { color: #008000; border: 1px dotted #008000 }\n.cython.code .py_attr { color: #FF0000; font-weight: bold; }\n.cython.code .c_attr { color: #0000FF; }\n.cython.code .py_call { color: #FF0000; font-weight: bold; }\n.cython.code .c_call { color: #0000FF; }\n\n.cython.score-0 {background-color: #FFFFff;}\n.cython.score-1 {background-color: #FFFFe7;}\n.cython.score-2 {background-color: #FFFFd4;}\n.cython.score-3 {background-color: #FFFFc4;}\n.cython.score-4 {background-color: #FFFFb6;}\n.cython.score-5 {background-color: #FFFFaa;}\n.cython.score-6 {background-color: #FFFF9f;}\n.cython.score-7 {background-color: #FFFF96;}\n.cython.score-8 {background-color: #FFFF8d;}\n.cython.score-9 {background-color: #FFFF86;}\n.cython.score-10 {background-color: #FFFF7f;}\n.cython.score-11 {background-color: #FFFF79;}\n.cython.score-12 {background-color: #FFFF73;}\n.cython.score-13 {background-color: #FFFF6e;}\n.cython.score-14 {background-color: #FFFF6a;}\n.cython.score-15 {background-color: #FFFF66;}\n.cython.score-16 {background-color: #FFFF62;}\n.cython.score-17 {background-color: #FFFF5e;}\n.cython.score-18 {background-color: #FFFF5b;}\n.cython.score-19 {background-color: #FFFF57;}\n.cython.score-20 {background-color: #FFFF55;}\n.cython.score-21 {background-color: #FFFF52;}\n.cython.score-22 {background-color: #FFFF4f;}\n.cython.score-23 {background-color: #FFFF4d;}\n.cython.score-24 {background-color: #FFFF4b;}\n.cython.score-25 {background-color: #FFFF48;}\n.cython.score-26 {background-color: #FFFF46;}\n.cython.score-27 {background-color: #FFFF44;}\n.cython.score-28 {background-color: #FFFF43;}\n.cython.score-29 {background-color: #FFFF41;}\n.cython.score-30 {background-color: #FFFF3f;}\n.cython.score-31 {background-color: #FFFF3e;}\n.cython.score-32 {background-color: #FFFF3c;}\n.cython.score-33 {background-color: #FFFF3b;}\n.cython.score-34 {background-color: #FFFF39;}\n.cython.score-35 {background-color: #FFFF38;}\n.cython.score-36 {background-color: #FFFF37;}\n.cython.score-37 {background-color: #FFFF36;}\n.cython.score-38 {background-color: #FFFF35;}\n.cython.score-39 {background-color: #FFFF34;}\n.cython.score-40 {background-color: #FFFF33;}\n.cython.score-41 {background-color: #FFFF32;}\n.cython.score-42 {background-color: #FFFF31;}\n.cython.score-43 {background-color: #FFFF30;}\n.cython.score-44 {background-color: #FFFF2f;}\n.cython.score-45 {background-color: #FFFF2e;}\n.cython.score-46 {background-color: #FFFF2d;}\n.cython.score-47 {background-color: #FFFF2c;}\n.cython.score-48 {background-color: #FFFF2b;}\n.cython.score-49 {background-color: #FFFF2b;}\n.cython.score-50 {background-color: #FFFF2a;}\n.cython.score-51 {background-color: #FFFF29;}\n.cython.score-52 {background-color: #FFFF29;}\n.cython.score-53 {background-color: #FFFF28;}\n.cython.score-54 {background-color: #FFFF27;}\n.cython.score-55 {background-color: #FFFF27;}\n.cython.score-56 {background-color: #FFFF26;}\n.cython.score-57 {background-color: #FFFF26;}\n.cython.score-58 {background-color: #FFFF25;}\n.cython.score-59 {background-color: #FFFF24;}\n.cython.score-60 {background-color: #FFFF24;}\n.cython.score-61 {background-color: #FFFF23;}\n.cython.score-62 {background-color: #FFFF23;}\n.cython.score-63 {background-color: #FFFF22;}\n.cython.score-64 {background-color: #FFFF22;}\n.cython.score-65 {background-color: #FFFF22;}\n.cython.score-66 {background-color: #FFFF21;}\n.cython.score-67 {background-color: #FFFF21;}\n.cython.score-68 {background-color: #FFFF20;}\n.cython.score-69 {background-color: #FFFF20;}\n.cython.score-70 {background-color: #FFFF1f;}\n.cython.score-71 {background-color: #FFFF1f;}\n.cython.score-72 {background-color: #FFFF1f;}\n.cython.score-73 {background-color: #FFFF1e;}\n.cython.score-74 {background-color: #FFFF1e;}\n.cython.score-75 {background-color: #FFFF1e;}\n.cython.score-76 {background-color: #FFFF1d;}\n.cython.score-77 {background-color: #FFFF1d;}\n.cython.score-78 {background-color: #FFFF1c;}\n.cython.score-79 {background-color: #FFFF1c;}\n.cython.score-80 {background-color: #FFFF1c;}\n.cython.score-81 {background-color: #FFFF1c;}\n.cython.score-82 {background-color: #FFFF1b;}\n.cython.score-83 {background-color: #FFFF1b;}\n.cython.score-84 {background-color: #FFFF1b;}\n.cython.score-85 {background-color: #FFFF1a;}\n.cython.score-86 {background-color: #FFFF1a;}\n.cython.score-87 {background-color: #FFFF1a;}\n.cython.score-88 {background-color: #FFFF1a;}\n.cython.score-89 {background-color: #FFFF19;}\n.cython.score-90 {background-color: #FFFF19;}\n.cython.score-91 {background-color: #FFFF19;}\n.cython.score-92 {background-color: #FFFF19;}\n.cython.score-93 {background-color: #FFFF18;}\n.cython.score-94 {background-color: #FFFF18;}\n.cython.score-95 {background-color: #FFFF18;}\n.cython.score-96 {background-color: #FFFF18;}\n.cython.score-97 {background-color: #FFFF17;}\n.cython.score-98 {background-color: #FFFF17;}\n.cython.score-99 {background-color: #FFFF17;}\n.cython.score-100 {background-color: #FFFF17;}\n.cython.score-101 {background-color: #FFFF16;}\n.cython.score-102 {background-color: #FFFF16;}\n.cython.score-103 {background-color: #FFFF16;}\n.cython.score-104 {background-color: #FFFF16;}\n.cython.score-105 {background-color: #FFFF16;}\n.cython.score-106 {background-color: #FFFF15;}\n.cython.score-107 {background-color: #FFFF15;}\n.cython.score-108 {background-color: #FFFF15;}\n.cython.score-109 {background-color: #FFFF15;}\n.cython.score-110 {background-color: #FFFF15;}\n.cython.score-111 {background-color: #FFFF15;}\n.cython.score-112 {background-color: #FFFF14;}\n.cython.score-113 {background-color: #FFFF14;}\n.cython.score-114 {background-color: #FFFF14;}\n.cython.score-115 {background-color: #FFFF14;}\n.cython.score-116 {background-color: #FFFF14;}\n.cython.score-117 {background-color: #FFFF14;}\n.cython.score-118 {background-color: #FFFF13;}\n.cython.score-119 {background-color: #FFFF13;}\n.cython.score-120 {background-color: #FFFF13;}\n.cython.score-121 {background-color: #FFFF13;}\n.cython.score-122 {background-color: #FFFF13;}\n.cython.score-123 {background-color: #FFFF13;}\n.cython.score-124 {background-color: #FFFF13;}\n.cython.score-125 {background-color: #FFFF12;}\n.cython.score-126 {background-color: #FFFF12;}\n.cython.score-127 {background-color: #FFFF12;}\n.cython.score-128 {background-color: #FFFF12;}\n.cython.score-129 {background-color: #FFFF12;}\n.cython.score-130 {background-color: #FFFF12;}\n.cython.score-131 {background-color: #FFFF12;}\n.cython.score-132 {background-color: #FFFF11;}\n.cython.score-133 {background-color: #FFFF11;}\n.cython.score-134 {background-color: #FFFF11;}\n.cython.score-135 {background-color: #FFFF11;}\n.cython.score-136 {background-color: #FFFF11;}\n.cython.score-137 {background-color: #FFFF11;}\n.cython.score-138 {background-color: #FFFF11;}\n.cython.score-139 {background-color: #FFFF11;}\n.cython.score-140 {background-color: #FFFF11;}\n.cython.score-141 {background-color: #FFFF10;}\n.cython.score-142 {background-color: #FFFF10;}\n.cython.score-143 {background-color: #FFFF10;}\n.cython.score-144 {background-color: #FFFF10;}\n.cython.score-145 {background-color: #FFFF10;}\n.cython.score-146 {background-color: #FFFF10;}\n.cython.score-147 {background-color: #FFFF10;}\n.cython.score-148 {background-color: #FFFF10;}\n.cython.score-149 {background-color: #FFFF10;}\n.cython.score-150 {background-color: #FFFF0f;}\n.cython.score-151 {background-color: #FFFF0f;}\n.cython.score-152 {background-color: #FFFF0f;}\n.cython.score-153 {background-color: #FFFF0f;}\n.cython.score-154 {background-color: #FFFF0f;}\n.cython.score-155 {background-color: #FFFF0f;}\n.cython.score-156 {background-color: #FFFF0f;}\n.cython.score-157 {background-color: #FFFF0f;}\n.cython.score-158 {background-color: #FFFF0f;}\n.cython.score-159 {background-color: #FFFF0f;}\n.cython.score-160 {background-color: #FFFF0f;}\n.cython.score-161 {background-color: #FFFF0e;}\n.cython.score-162 {background-color: #FFFF0e;}\n.cython.score-163 {background-color: #FFFF0e;}\n.cython.score-164 {background-color: #FFFF0e;}\n.cython.score-165 {background-color: #FFFF0e;}\n.cython.score-166 {background-color: #FFFF0e;}\n.cython.score-167 {background-color: #FFFF0e;}\n.cython.score-168 {background-color: #FFFF0e;}\n.cython.score-169 {background-color: #FFFF0e;}\n.cython.score-170 {background-color: #FFFF0e;}\n.cython.score-171 {background-color: #FFFF0e;}\n.cython.score-172 {background-color: #FFFF0e;}\n.cython.score-173 {background-color: #FFFF0d;}\n.cython.score-174 {background-color: #FFFF0d;}\n.cython.score-175 {background-color: #FFFF0d;}\n.cython.score-176 {background-color: #FFFF0d;}\n.cython.score-177 {background-color: #FFFF0d;}\n.cython.score-178 {background-color: #FFFF0d;}\n.cython.score-179 {background-color: #FFFF0d;}\n.cython.score-180 {background-color: #FFFF0d;}\n.cython.score-181 {background-color: #FFFF0d;}\n.cython.score-182 {background-color: #FFFF0d;}\n.cython.score-183 {background-color: #FFFF0d;}\n.cython.score-184 {background-color: #FFFF0d;}\n.cython.score-185 {background-color: #FFFF0d;}\n.cython.score-186 {background-color: #FFFF0d;}\n.cython.score-187 {background-color: #FFFF0c;}\n.cython.score-188 {background-color: #FFFF0c;}\n.cython.score-189 {background-color: #FFFF0c;}\n.cython.score-190 {background-color: #FFFF0c;}\n.cython.score-191 {background-color: #FFFF0c;}\n.cython.score-192 {background-color: #FFFF0c;}\n.cython.score-193 {background-color: #FFFF0c;}\n.cython.score-194 {background-color: #FFFF0c;}\n.cython.score-195 {background-color: #FFFF0c;}\n.cython.score-196 {background-color: #FFFF0c;}\n.cython.score-197 {background-color: #FFFF0c;}\n.cython.score-198 {background-color: #FFFF0c;}\n.cython.score-199 {background-color: #FFFF0c;}\n.cython.score-200 {background-color: #FFFF0c;}\n.cython.score-201 {background-color: #FFFF0c;}\n.cython.score-202 {background-color: #FFFF0c;}\n.cython.score-203 {background-color: #FFFF0b;}\n.cython.score-204 {background-color: #FFFF0b;}\n.cython.score-205 {background-color: #FFFF0b;}\n.cython.score-206 {background-color: #FFFF0b;}\n.cython.score-207 {background-color: #FFFF0b;}\n.cython.score-208 {background-color: #FFFF0b;}\n.cython.score-209 {background-color: #FFFF0b;}\n.cython.score-210 {background-color: #FFFF0b;}\n.cython.score-211 {background-color: #FFFF0b;}\n.cython.score-212 {background-color: #FFFF0b;}\n.cython.score-213 {background-color: #FFFF0b;}\n.cython.score-214 {background-color: #FFFF0b;}\n.cython.score-215 {background-color: #FFFF0b;}\n.cython.score-216 {background-color: #FFFF0b;}\n.cython.score-217 {background-color: #FFFF0b;}\n.cython.score-218 {background-color: #FFFF0b;}\n.cython.score-219 {background-color: #FFFF0b;}\n.cython.score-220 {background-color: #FFFF0b;}\n.cython.score-221 {background-color: #FFFF0b;}\n.cython.score-222 {background-color: #FFFF0a;}\n.cython.score-223 {background-color: #FFFF0a;}\n.cython.score-224 {background-color: #FFFF0a;}\n.cython.score-225 {background-color: #FFFF0a;}\n.cython.score-226 {background-color: #FFFF0a;}\n.cython.score-227 {background-color: #FFFF0a;}\n.cython.score-228 {background-color: #FFFF0a;}\n.cython.score-229 {background-color: #FFFF0a;}\n.cython.score-230 {background-color: #FFFF0a;}\n.cython.score-231 {background-color: #FFFF0a;}\n.cython.score-232 {background-color: #FFFF0a;}\n.cython.score-233 {background-color: #FFFF0a;}\n.cython.score-234 {background-color: #FFFF0a;}\n.cython.score-235 {background-color: #FFFF0a;}\n.cython.score-236 {background-color: #FFFF0a;}\n.cython.score-237 {background-color: #FFFF0a;}\n.cython.score-238 {background-color: #FFFF0a;}\n.cython.score-239 {background-color: #FFFF0a;}\n.cython.score-240 {background-color: #FFFF0a;}\n.cython.score-241 {background-color: #FFFF0a;}\n.cython.score-242 {background-color: #FFFF0a;}\n.cython.score-243 {background-color: #FFFF0a;}\n.cython.score-244 {background-color: #FFFF0a;}\n.cython.score-245 {background-color: #FFFF0a;}\n.cython.score-246 {background-color: #FFFF09;}\n.cython.score-247 {background-color: #FFFF09;}\n.cython.score-248 {background-color: #FFFF09;}\n.cython.score-249 {background-color: #FFFF09;}\n.cython.score-250 {background-color: #FFFF09;}\n.cython.score-251 {background-color: #FFFF09;}\n.cython.score-252 {background-color: #FFFF09;}\n.cython.score-253 {background-color: #FFFF09;}\n.cython.score-254 {background-color: #FFFF09;}\n.cython .hll { background-color: #ffffcc }\n.cython { background: #f8f8f8; }\n.cython .c { color: #408080; font-style: italic } /* Comment */\n.cython .err { border: 1px solid #FF0000 } /* Error */\n.cython .k { color: #008000; font-weight: bold } /* Keyword */\n.cython .o { color: #666666 } /* Operator */\n.cython .ch { color: #408080; font-style: italic } /* Comment.Hashbang */\n.cython .cm { color: #408080; font-style: italic } /* Comment.Multiline */\n.cython .cp { color: #BC7A00 } /* Comment.Preproc */\n.cython .cpf { color: #408080; font-style: italic } /* Comment.PreprocFile */\n.cython .c1 { color: #408080; font-style: italic } /* Comment.Single */\n.cython .cs { color: #408080; font-style: italic } /* Comment.Special */\n.cython .gd { color: #A00000 } /* Generic.Deleted */\n.cython .ge { font-style: italic } /* Generic.Emph */\n.cython .gr { color: #FF0000 } /* Generic.Error */\n.cython .gh { color: #000080; font-weight: bold } /* Generic.Heading */\n.cython .gi { color: #00A000 } /* Generic.Inserted */\n.cython .go { color: #888888 } /* Generic.Output */\n.cython .gp { color: #000080; font-weight: bold } /* Generic.Prompt */\n.cython .gs { font-weight: bold } /* Generic.Strong */\n.cython .gu { color: #800080; font-weight: bold } /* Generic.Subheading */\n.cython .gt { color: #0044DD } /* Generic.Traceback */\n.cython .kc { color: #008000; font-weight: bold } /* Keyword.Constant */\n.cython .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */\n.cython .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */\n.cython .kp { color: #008000 } /* Keyword.Pseudo */\n.cython .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */\n.cython .kt { color: #B00040 } /* Keyword.Type */\n.cython .m { color: #666666 } /* Literal.Number */\n.cython .s { color: #BA2121 } /* Literal.String */\n.cython .na { color: #7D9029 } /* Name.Attribute */\n.cython .nb { color: #008000 } /* Name.Builtin */\n.cython .nc { color: #0000FF; font-weight: bold } /* Name.Class */\n.cython .no { color: #880000 } /* Name.Constant */\n.cython .nd { color: #AA22FF } /* Name.Decorator */\n.cython .ni { color: #999999; font-weight: bold } /* Name.Entity */\n.cython .ne { color: #D2413A; font-weight: bold } /* Name.Exception */\n.cython .nf { color: #0000FF } /* Name.Function */\n.cython .nl { color: #A0A000 } /* Name.Label */\n.cython .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */\n.cython .nt { color: #008000; font-weight: bold } /* Name.Tag */\n.cython .nv { color: #19177C } /* Name.Variable */\n.cython .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */\n.cython .w { color: #bbbbbb } /* Text.Whitespace */\n.cython .mb { color: #666666 } /* Literal.Number.Bin */\n.cython .mf { color: #666666 } /* Literal.Number.Float */\n.cython .mh { color: #666666 } /* Literal.Number.Hex */\n.cython .mi { color: #666666 } /* Literal.Number.Integer */\n.cython .mo { color: #666666 } /* Literal.Number.Oct */\n.cython .sa { color: #BA2121 } /* Literal.String.Affix */\n.cython .sb { color: #BA2121 } /* Literal.String.Backtick */\n.cython .sc { color: #BA2121 } /* Literal.String.Char */\n.cython .dl { color: #BA2121 } /* Literal.String.Delimiter */\n.cython .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */\n.cython .s2 { color: #BA2121 } /* Literal.String.Double */\n.cython .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */\n.cython .sh { color: #BA2121 } /* Literal.String.Heredoc */\n.cython .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */\n.cython .sx { color: #008000 } /* Literal.String.Other */\n.cython .sr { color: #BB6688 } /* Literal.String.Regex */\n.cython .s1 { color: #BA2121 } /* Literal.String.Single */\n.cython .ss { color: #19177C } /* Literal.String.Symbol */\n.cython .bp { color: #008000 } /* Name.Builtin.Pseudo */\n.cython .fm { color: #0000FF } /* Name.Function.Magic */\n.cython .vc { color: #19177C } /* Name.Variable.Class */\n.cython .vg { color: #19177C } /* Name.Variable.Global */\n.cython .vi { color: #19177C } /* Name.Variable.Instance */\n.cython .vm { color: #19177C } /* Name.Variable.Magic */\n.cython .il { color: #666666 } /* Literal.Number.Integer.Long */\n </style>\n</head>\n<body class=\"cython\">\n<p><span style=\"border-bottom: solid 1px grey;\">Generated by Cython 0.29.15</span></p>\n<p>\n <span style=\"background-color: #FFFF00\">Yellow lines</span> hint at Python interaction.<br />\n Click on a line that starts with a \"<code>+</code>\" to see the C code that Cython generated for it.\n</p>\n<div class=\"cython\"><pre class=\"cython line score-8\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">01</span>: <span class=\"c\"># cython: boundscheck=False, wraparound=False, cdivision=True</span></pre>\n<pre class='cython code score-8 '> __pyx_t_1 = <span class='pyx_c_api'>__Pyx_PyDict_NewPresized</span>(0);<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 1, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);\n if (<span class='py_c_api'>PyDict_SetItem</span>(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) <span class='error_goto'>__PYX_ERR(0, 1, __pyx_L1_error)</span>\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;\n</pre><pre class=\"cython line score-0\"> <span class=\"\">02</span>: </pre>\n<pre class=\"cython line score-0\"> <span class=\"\">03</span>: <span class=\"k\">cimport</span> <span class=\"nn\">numpy</span> <span class=\"k\">as</span> <span class=\"nn\">np</span></pre>\n<pre class=\"cython line score-8\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">04</span>: <span class=\"k\">import</span> <span class=\"nn\">numpy</span> <span class=\"k\">as</span> <span class=\"nn\">np</span></pre>\n<pre class='cython code score-8 '> __pyx_t_1 = <span class='pyx_c_api'>__Pyx_Import</span>(__pyx_n_s_numpy, 0, 0);<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 4, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);\n if (<span class='py_c_api'>PyDict_SetItem</span>(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) <span class='error_goto'>__PYX_ERR(0, 4, __pyx_L1_error)</span>\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;\n</pre><pre class=\"cython line score-0\"> <span class=\"\">05</span>: </pre>\n<pre class=\"cython line score-0\"> <span class=\"\">06</span>: <span class=\"k\">cimport</span> <span class=\"nn\">cython</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">07</span>: <span class=\"k\">from</span> <span class=\"nn\">cython.parallel</span> <span class=\"k\">cimport</span> <span class=\"n\">prange</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">08</span>: </pre>\n<pre class=\"cython line score-0\"> <span class=\"\">09</span>: <span class=\"c\"># from libc.math cimport abs, sqrt, fmax</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">10</span>: <span class=\"c\"># from libcpp.complex cimport pow as cpow, real, imag, abs</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">11</span>: </pre>\n<pre class=\"cython line score-0\"> <span class=\"\">12</span>: <span class=\"k\">cdef</span> <span class=\"kr\">extern</span> <span class=\"k\">from</span> <span class=\"s\">"complex.h"</span> <span class=\"k\">nogil</span><span class=\"p\">:</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">13</span>: <span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">complex128_t</span> <span class=\"n\">cpow</span><span class=\"p\">(</span><span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">complex128_t</span><span class=\"p\">,</span> <span class=\"nb\">long</span><span class=\"p\">)</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">14</span>: <span class=\"n\">double</span> <span class=\"n\">cabs</span><span class=\"p\">(</span><span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">complex128_t</span><span class=\"p\">)</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">15</span>: <span class=\"n\">double</span> <span class=\"n\">creal</span><span class=\"p\">(</span><span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">complex128_t</span><span class=\"p\">)</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">16</span>: <span class=\"n\">double</span> <span class=\"n\">cimag</span><span class=\"p\">(</span><span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">complex128_t</span><span class=\"p\">)</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">17</span>: </pre>\n<pre class=\"cython line score-0\"> <span class=\"\">18</span>: <span class=\"k\">cdef</span> <span class=\"kr\">extern</span> <span class=\"k\">from</span> <span class=\"s\">"math.h"</span> <span class=\"k\">nogil</span><span class=\"p\">:</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">19</span>: <span class=\"nb\">long</span> <span class=\"n\">double</span> <span class=\"n\">sqrtl</span><span class=\"p\">(</span><span class=\"nb\">long</span> <span class=\"n\">double</span><span class=\"p\">)</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">20</span>: <span class=\"nb\">long</span> <span class=\"n\">double</span> <span class=\"n\">fabsl</span><span class=\"p\">(</span><span class=\"nb\">long</span> <span class=\"n\">double</span><span class=\"p\">)</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">21</span>: <span class=\"nb\">long</span> <span class=\"n\">double</span> <span class=\"n\">fmaxl</span><span class=\"p\">(</span><span class=\"nb\">long</span> <span class=\"n\">double</span><span class=\"p\">,</span> <span class=\"nb\">long</span> <span class=\"n\">double</span><span class=\"p\">)</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">22</span>: </pre>\n<pre class=\"cython line score-0\"> <span class=\"\">23</span>: </pre>\n<pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">24</span>: <span class=\"k\">cdef</span> <span class=\"kt\">np</span>.<span class=\"kt\">complex128_t</span> <span class=\"nf\">f_c</span><span class=\"p\">(</span><span class=\"n\">double</span><span class=\"p\">[:]</span> <span class=\"n\">poly</span><span class=\"p\">,</span> <span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">complex128_t</span> <span class=\"n\">x</span><span class=\"p\">)</span> <span class=\"k\">nogil</span><span class=\"p\">:</span></pre>\n<pre class='cython code score-0 '>static __pyx_t_double_complex __pyx_f_46_cython_magic_7d2d8e83eb391e7da013a05885fe9a7f_f_c(__Pyx_memviewslice __pyx_v_poly, __pyx_t_double_complex __pyx_v_x) {\n __pyx_t_double_complex __pyx_v_out;\n long __pyx_v_i;\n long __pyx_v_k;\n long __pyx_v_L;\n __pyx_t_double_complex __pyx_r;\n/* … */\n /* function exit code */\n __pyx_L0:;\n return __pyx_r;\n}\n</pre><pre class=\"cython line score-0\"> <span class=\"\">25</span>: <span class=\"k\">cdef</span><span class=\"p\">:</span></pre>\n<pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">26</span>: <span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">complex128_t</span> <span class=\"n\">out</span> <span class=\"o\">=</span> <span class=\"mf\">0</span> <span class=\"o\">+</span> <span class=\"mf\">0</span><span class=\"n\">j</span></pre>\n<pre class='cython code score-0 '> __pyx_v_out = __Pyx_c_sum_double(__pyx_t_double_complex_from_parts(0, 0), __pyx_t_double_complex_from_parts(0, 0.0));\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">27</span>: <span class=\"nb\">long</span> <span class=\"n\">i</span><span class=\"p\">,</span> <span class=\"n\">k</span><span class=\"p\">,</span> <span class=\"n\">L</span> <span class=\"o\">=</span> <span class=\"n\">poly</span><span class=\"o\">.</span><span class=\"n\">shape</span><span class=\"p\">[</span><span class=\"mf\">0</span><span class=\"p\">]</span></pre>\n<pre class='cython code score-0 '> __pyx_v_L = (__pyx_v_poly.shape[0]);\n</pre><pre class=\"cython line score-0\"> <span class=\"\">28</span>: </pre>\n<pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">29</span>: <span class=\"k\">for</span> <span class=\"n\">i</span> <span class=\"ow\">in</span> <span class=\"nb\">range</span><span class=\"p\">(</span><span class=\"n\">L</span><span class=\"p\">):</span></pre>\n<pre class='cython code score-0 '> __pyx_t_1 = __pyx_v_L;\n __pyx_t_2 = __pyx_t_1;\n for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) {\n __pyx_v_i = __pyx_t_3;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">30</span>: <span class=\"n\">k</span> <span class=\"o\">=</span> <span class=\"n\">L</span> <span class=\"o\">-</span> <span class=\"n\">i</span> <span class=\"o\">-</span> <span class=\"mf\">1</span></pre>\n<pre class='cython code score-0 '> __pyx_v_k = ((__pyx_v_L - __pyx_v_i) - 1);\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">31</span>: <span class=\"k\">if</span> <span class=\"n\">poly</span><span class=\"p\">[</span><span class=\"n\">k</span><span class=\"p\">]</span> <span class=\"o\">!=</span> <span class=\"mf\">0</span><span class=\"p\">:</span></pre>\n<pre class='cython code score-0 '> __pyx_t_4 = __pyx_v_k;\n __pyx_t_5 = (((*((double *) ( /* dim=0 */ (__pyx_v_poly.data + __pyx_t_4 * __pyx_v_poly.strides[0]) ))) != 0.0) != 0);\n if (__pyx_t_5) {\n/* … */\n }\n }\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">32</span>: <span class=\"n\">out</span> <span class=\"o\">+=</span> <span class=\"n\">poly</span><span class=\"p\">[</span><span class=\"n\">k</span><span class=\"p\">]</span> <span class=\"o\">*</span> <span class=\"n\">cpow</span><span class=\"p\">(</span><span class=\"n\">x</span><span class=\"p\">,</span> <span class=\"n\">i</span><span class=\"p\">)</span></pre>\n<pre class='cython code score-0 '> __pyx_t_6 = __pyx_v_k;\n __pyx_v_out = __Pyx_c_sum_double(__pyx_v_out, __Pyx_c_prod_double(__pyx_t_double_complex_from_parts((*((double *) ( /* dim=0 */ (__pyx_v_poly.data + __pyx_t_6 * __pyx_v_poly.strides[0]) ))), 0), cpow(__pyx_v_x, __pyx_v_i)));\n</pre><pre class=\"cython line score-0\"> <span class=\"\">33</span>: <span class=\"c\"># if poly[i] != 0:</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">34</span>: <span class=\"c\"># out += <np.complex128_t> poly[i] * (x ** i)</span></pre>\n<pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">35</span>: <span class=\"k\">return</span> <span class=\"n\">out</span></pre>\n<pre class='cython code score-0 '> __pyx_r = __pyx_v_out;\n goto __pyx_L0;\n</pre><pre class=\"cython line score-0\"> <span class=\"\">36</span>: </pre>\n<pre class=\"cython line score-0\"> <span class=\"\">37</span>: </pre>\n<pre class=\"cython line score-60\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">38</span>: <span class=\"k\">def</span> <span class=\"nf\">solve_c</span><span class=\"p\">(</span><span class=\"n\">double</span><span class=\"p\">[:]</span> <span class=\"n\">poly</span><span class=\"p\">,</span> <span class=\"nb\">long</span> <span class=\"n\">max_iter</span> <span class=\"o\">=</span> <span class=\"mf\">500000</span><span class=\"p\">):</span></pre>\n<pre class='cython code score-60 '>/* Python wrapper */\nstatic PyObject *__pyx_pw_46_cython_magic_7d2d8e83eb391e7da013a05885fe9a7f_1solve_c(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/\nstatic PyMethodDef __pyx_mdef_46_cython_magic_7d2d8e83eb391e7da013a05885fe9a7f_1solve_c = {\"solve_c\", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_46_cython_magic_7d2d8e83eb391e7da013a05885fe9a7f_1solve_c, METH_VARARGS|METH_KEYWORDS, 0};\nstatic PyObject *__pyx_pw_46_cython_magic_7d2d8e83eb391e7da013a05885fe9a7f_1solve_c(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {\n __Pyx_memviewslice __pyx_v_poly = { 0, 0, { 0 }, { 0 }, { 0 } };\n long __pyx_v_max_iter;\n PyObject *__pyx_r = 0;\n <span class='refnanny'>__Pyx_RefNannyDeclarations</span>\n <span class='refnanny'>__Pyx_RefNannySetupContext</span>(\"solve_c (wrapper)\", 0);\n {\n static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_poly,&__pyx_n_s_max_iter,0};\n PyObject* values[2] = {0,0};\n if (unlikely(__pyx_kwds)) {\n Py_ssize_t kw_args;\n const Py_ssize_t pos_args = <span class='py_macro_api'>PyTuple_GET_SIZE</span>(__pyx_args);\n switch (pos_args) {\n case 2: values[1] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 1);\n CYTHON_FALLTHROUGH;\n case 1: values[0] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 0);\n CYTHON_FALLTHROUGH;\n case 0: break;\n default: goto __pyx_L5_argtuple_error;\n }\n kw_args = <span class='py_c_api'>PyDict_Size</span>(__pyx_kwds);\n switch (pos_args) {\n case 0:\n if (likely((values[0] = <span class='pyx_c_api'>__Pyx_PyDict_GetItemStr</span>(__pyx_kwds, __pyx_n_s_poly)) != 0)) kw_args--;\n else goto __pyx_L5_argtuple_error;\n CYTHON_FALLTHROUGH;\n case 1:\n if (kw_args > 0) {\n PyObject* value = <span class='pyx_c_api'>__Pyx_PyDict_GetItemStr</span>(__pyx_kwds, __pyx_n_s_max_iter);\n if (value) { values[1] = value; kw_args--; }\n }\n }\n if (unlikely(kw_args > 0)) {\n if (unlikely(<span class='pyx_c_api'>__Pyx_ParseOptionalKeywords</span>(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, \"solve_c\") < 0)) <span class='error_goto'>__PYX_ERR(0, 38, __pyx_L3_error)</span>\n }\n } else {\n switch (<span class='py_macro_api'>PyTuple_GET_SIZE</span>(__pyx_args)) {\n case 2: values[1] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 1);\n CYTHON_FALLTHROUGH;\n case 1: values[0] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 0);\n break;\n default: goto __pyx_L5_argtuple_error;\n }\n }\n __pyx_v_poly = <span class='pyx_c_api'>__Pyx_PyObject_to_MemoryviewSlice_ds_double</span>(values[0], PyBUF_WRITABLE);<span class='error_goto'> if (unlikely(!__pyx_v_poly.memview)) __PYX_ERR(0, 38, __pyx_L3_error)</span>\n if (values[1]) {\n __pyx_v_max_iter = <span class='pyx_c_api'>__Pyx_PyInt_As_long</span>(values[1]); if (unlikely((__pyx_v_max_iter == (long)-1) && <span class='py_c_api'>PyErr_Occurred</span>())) <span class='error_goto'>__PYX_ERR(0, 38, __pyx_L3_error)</span>\n } else {\n __pyx_v_max_iter = ((long)0x7A120);\n }\n }\n goto __pyx_L4_argument_unpacking_done;\n __pyx_L5_argtuple_error:;\n <span class='pyx_c_api'>__Pyx_RaiseArgtupleInvalid</span>(\"solve_c\", 0, 1, 2, <span class='py_macro_api'>PyTuple_GET_SIZE</span>(__pyx_args)); <span class='error_goto'>__PYX_ERR(0, 38, __pyx_L3_error)</span>\n __pyx_L3_error:;\n <span class='pyx_c_api'>__Pyx_AddTraceback</span>(\"_cython_magic_7d2d8e83eb391e7da013a05885fe9a7f.solve_c\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n <span class='refnanny'>__Pyx_RefNannyFinishContext</span>();\n return NULL;\n __pyx_L4_argument_unpacking_done:;\n __pyx_r = __pyx_pf_46_cython_magic_7d2d8e83eb391e7da013a05885fe9a7f_solve_c(__pyx_self, __pyx_v_poly, __pyx_v_max_iter);\n\n /* function exit code */\n <span class='refnanny'>__Pyx_RefNannyFinishContext</span>();\n return __pyx_r;\n}\n\nstatic PyObject *__pyx_pf_46_cython_magic_7d2d8e83eb391e7da013a05885fe9a7f_solve_c(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_poly, long __pyx_v_max_iter) {\n double __pyx_v_epsilon;\n long __pyx_v_ii;\n long __pyx_v_i;\n long __pyx_v_k;\n long __pyx_v_L;\n double __pyx_v_scale;\n PyArrayObject *__pyx_v_roots = 0;\n __pyx_t_double_complex __pyx_v_delta;\n long double __pyx_v_max_delta;\n long double __pyx_v_rd;\n long double __pyx_v_id;\n __Pyx_LocalBuf_ND __pyx_pybuffernd_roots;\n __Pyx_Buffer __pyx_pybuffer_roots;\n PyObject *__pyx_r = NULL;\n <span class='refnanny'>__Pyx_RefNannyDeclarations</span>\n <span class='refnanny'>__Pyx_RefNannySetupContext</span>(\"solve_c\", 0);\n __pyx_pybuffer_roots.pybuffer.buf = NULL;\n __pyx_pybuffer_roots.refcount = 0;\n __pyx_pybuffernd_roots.data = NULL;\n __pyx_pybuffernd_roots.rcbuffer = &__pyx_pybuffer_roots;\n/* … */\n /* function exit code */\n __pyx_L1_error:;\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_1);\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_2);\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_3);\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_4);\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_5);\n { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;\n __Pyx_PyThreadState_declare\n __Pyx_PyThreadState_assign\n <span class='pyx_c_api'>__Pyx_ErrFetch</span>(&__pyx_type, &__pyx_value, &__pyx_tb);\n <span class='pyx_c_api'>__Pyx_SafeReleaseBuffer</span>(&__pyx_pybuffernd_roots.rcbuffer->pybuffer);\n <span class='pyx_c_api'>__Pyx_ErrRestore</span>(__pyx_type, __pyx_value, __pyx_tb);}\n <span class='pyx_c_api'>__Pyx_AddTraceback</span>(\"_cython_magic_7d2d8e83eb391e7da013a05885fe9a7f.solve_c\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n __pyx_r = NULL;\n goto __pyx_L2;\n __pyx_L0:;\n <span class='pyx_c_api'>__Pyx_SafeReleaseBuffer</span>(&__pyx_pybuffernd_roots.rcbuffer->pybuffer);\n __pyx_L2:;\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>((PyObject *)__pyx_v_roots);\n __PYX_XDEC_MEMVIEW(&__pyx_v_poly, 1);\n <span class='refnanny'>__Pyx_XGIVEREF</span>(__pyx_r);\n <span class='refnanny'>__Pyx_RefNannyFinishContext</span>();\n return __pyx_r;\n}\n/* … */\n __pyx_tuple__26 = <span class='py_c_api'>PyTuple_Pack</span>(15, __pyx_n_s_poly, __pyx_n_s_max_iter, __pyx_n_s_epsilon, __pyx_n_s_ii, __pyx_n_s_i, __pyx_n_s_k, __pyx_n_s_L, __pyx_n_s_scale, __pyx_n_s_roots, __pyx_n_s_deno, __pyx_n_s_delta, __pyx_n_s_max_delta, __pyx_n_s_d2, __pyx_n_s_rd, __pyx_n_s_id);<span class='error_goto'> if (unlikely(!__pyx_tuple__26)) __PYX_ERR(0, 38, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_tuple__26);\n <span class='refnanny'>__Pyx_GIVEREF</span>(__pyx_tuple__26);\n/* … */\n __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_7d2d8e83eb391e7da013a05885fe9a7f_1solve_c, NULL, __pyx_n_s_cython_magic_7d2d8e83eb391e7da0);<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 38, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);\n if (<span class='py_c_api'>PyDict_SetItem</span>(__pyx_d, __pyx_n_s_solve_c, __pyx_t_1) < 0) <span class='error_goto'>__PYX_ERR(0, 38, __pyx_L1_error)</span>\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;\n __pyx_codeobj__27 = (PyObject*)<span class='pyx_c_api'>__Pyx_PyCode_New</span>(2, 0, 15, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__26, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_kvedala_ipython_cython__c, __pyx_n_s_solve_c, 38, __pyx_empty_bytes);<span class='error_goto'> if (unlikely(!__pyx_codeobj__27)) __PYX_ERR(0, 38, __pyx_L1_error)</span>\n</pre><pre class=\"cython line score-0\"> <span class=\"\">39</span>: <span class=\"k\">cdef</span><span class=\"p\">:</span></pre>\n<pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">40</span>: <span class=\"n\">double</span> <span class=\"n\">epsilon</span> <span class=\"o\">=</span> <span class=\"mf\">5e-10</span></pre>\n<pre class='cython code score-0 '> __pyx_v_epsilon = 5e-10;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">41</span>: <span class=\"nb\">long</span> <span class=\"n\">ii</span> <span class=\"o\">=</span> <span class=\"mf\">0</span><span class=\"p\">,</span> <span class=\"n\">i</span><span class=\"p\">,</span> <span class=\"n\">k</span><span class=\"p\">,</span> <span class=\"n\">L</span> <span class=\"o\">=</span> <span class=\"n\">poly</span><span class=\"o\">.</span><span class=\"n\">shape</span><span class=\"p\">[</span><span class=\"mf\">0</span><span class=\"p\">]</span> <span class=\"o\">-</span> <span class=\"mf\">1</span></pre>\n<pre class='cython code score-0 '> __pyx_v_ii = 0;\n __pyx_v_L = ((__pyx_v_poly.shape[0]) - 1);\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">42</span>: <span class=\"n\">double</span> <span class=\"n\">scale</span> <span class=\"o\">=</span> <span class=\"mf\">500.0</span></pre>\n<pre class='cython code score-0 '> __pyx_v_scale = 500.0;\n</pre><pre class=\"cython line score-69\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">43</span>: <span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">ndarray</span><span class=\"p\">[</span><span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">complex128_t</span><span class=\"p\">,</span> <span class=\"n\">ndim</span><span class=\"o\">=</span><span class=\"mf\">1</span><span class=\"p\">]</span> <span class=\"n\">roots</span> <span class=\"o\">=</span> <span class=\"p\">((</span><span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">random</span><span class=\"o\">.</span><span class=\"n\">random</span><span class=\"p\">(</span><span class=\"n\">L</span><span class=\"p\">)</span> <span class=\"o\">+</span> <span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">random</span><span class=\"o\">.</span><span class=\"n\">random</span><span class=\"p\">(</span><span class=\"n\">L</span><span class=\"p\">)</span> <span class=\"o\">*</span> <span class=\"mf\">1</span><span class=\"n\">j</span><span class=\"p\">)</span> \\</pre>\n<pre class='cython code score-69 '> <span class='pyx_c_api'>__Pyx_GetModuleGlobalName</span>(__pyx_t_2, __pyx_n_s_np);<span class='error_goto'> if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_2);\n __pyx_t_3 = <span class='pyx_c_api'>__Pyx_PyObject_GetAttrStr</span>(__pyx_t_2, __pyx_n_s_random);<span class='error_goto'> if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_3);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_2); __pyx_t_2 = 0;\n __pyx_t_2 = <span class='pyx_c_api'>__Pyx_PyObject_GetAttrStr</span>(__pyx_t_3, __pyx_n_s_random);<span class='error_goto'> if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_2);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_3); __pyx_t_3 = 0;\n __pyx_t_3 = <span class='pyx_c_api'>__Pyx_PyInt_From_long</span>(__pyx_v_L);<span class='error_goto'> if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_3);\n __pyx_t_4 = NULL;\n if (CYTHON_UNPACK_METHODS && likely(<span class='py_c_api'>PyMethod_Check</span>(__pyx_t_2))) {\n __pyx_t_4 = <span class='py_macro_api'>PyMethod_GET_SELF</span>(__pyx_t_2);\n if (likely(__pyx_t_4)) {\n PyObject* function = <span class='py_macro_api'>PyMethod_GET_FUNCTION</span>(__pyx_t_2);\n <span class='pyx_macro_api'>__Pyx_INCREF</span>(__pyx_t_4);\n <span class='pyx_macro_api'>__Pyx_INCREF</span>(function);\n <span class='pyx_macro_api'>__Pyx_DECREF_SET</span>(__pyx_t_2, function);\n }\n }\n __pyx_t_1 = (__pyx_t_4) ? __Pyx_PyObject_Call2Args(__pyx_t_2, __pyx_t_4, __pyx_t_3) : <span class='pyx_c_api'>__Pyx_PyObject_CallOneArg</span>(__pyx_t_2, __pyx_t_3);\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_4); __pyx_t_4 = 0;\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_3); __pyx_t_3 = 0;\n if (unlikely(!__pyx_t_1)) <span class='error_goto'>__PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_2); __pyx_t_2 = 0;\n <span class='pyx_c_api'>__Pyx_GetModuleGlobalName</span>(__pyx_t_3, __pyx_n_s_np);<span class='error_goto'> if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_3);\n __pyx_t_4 = <span class='pyx_c_api'>__Pyx_PyObject_GetAttrStr</span>(__pyx_t_3, __pyx_n_s_random);<span class='error_goto'> if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_4);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_3); __pyx_t_3 = 0;\n __pyx_t_3 = <span class='pyx_c_api'>__Pyx_PyObject_GetAttrStr</span>(__pyx_t_4, __pyx_n_s_random);<span class='error_goto'> if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_3);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_4); __pyx_t_4 = 0;\n __pyx_t_4 = <span class='pyx_c_api'>__Pyx_PyInt_From_long</span>(__pyx_v_L);<span class='error_goto'> if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_4);\n __pyx_t_5 = NULL;\n if (CYTHON_UNPACK_METHODS && likely(<span class='py_c_api'>PyMethod_Check</span>(__pyx_t_3))) {\n __pyx_t_5 = <span class='py_macro_api'>PyMethod_GET_SELF</span>(__pyx_t_3);\n if (likely(__pyx_t_5)) {\n PyObject* function = <span class='py_macro_api'>PyMethod_GET_FUNCTION</span>(__pyx_t_3);\n <span class='pyx_macro_api'>__Pyx_INCREF</span>(__pyx_t_5);\n <span class='pyx_macro_api'>__Pyx_INCREF</span>(function);\n <span class='pyx_macro_api'>__Pyx_DECREF_SET</span>(__pyx_t_3, function);\n }\n }\n __pyx_t_2 = (__pyx_t_5) ? __Pyx_PyObject_Call2Args(__pyx_t_3, __pyx_t_5, __pyx_t_4) : <span class='pyx_c_api'>__Pyx_PyObject_CallOneArg</span>(__pyx_t_3, __pyx_t_4);\n <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_5); __pyx_t_5 = 0;\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_4); __pyx_t_4 = 0;\n if (unlikely(!__pyx_t_2)) <span class='error_goto'>__PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_2);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_3); __pyx_t_3 = 0;\n __pyx_t_3 = <span class='py_c_api'>PyComplex_FromDoubles</span>(0.0, 1.0);<span class='error_goto'> if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_3);\n __pyx_t_4 = <span class='py_c_api'>PyNumber_Multiply</span>(__pyx_t_2, __pyx_t_3);<span class='error_goto'> if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_4);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_2); __pyx_t_2 = 0;\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_3); __pyx_t_3 = 0;\n __pyx_t_3 = <span class='py_c_api'>PyNumber_Add</span>(__pyx_t_1, __pyx_t_4);<span class='error_goto'> if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 43, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_3);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_4); __pyx_t_4 = 0;\n</pre><pre class=\"cython line score-29\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">44</span>: <span class=\"o\">-</span> <span class=\"p\">(</span><span class=\"o\">.</span><span class=\"mf\">5</span> <span class=\"o\">+</span> <span class=\"o\">.</span><span class=\"mf\">5</span><span class=\"n\">j</span><span class=\"p\">))</span> <span class=\"o\">*</span> <span class=\"n\">scale</span></pre>\n<pre class='cython code score-29 '> __pyx_t_6 = __Pyx_c_sum_double(__pyx_t_double_complex_from_parts(.5, 0), __pyx_t_double_complex_from_parts(0, 0.5));\n __pyx_t_4 = __pyx_<span class='py_c_api'>PyComplex_FromComplex</span>(__pyx_t_6);<span class='error_goto'> if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 44, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_4);\n __pyx_t_1 = <span class='py_c_api'>PyNumber_Subtract</span>(__pyx_t_3, __pyx_t_4);<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 44, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_3); __pyx_t_3 = 0;\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_4); __pyx_t_4 = 0;\n __pyx_t_4 = <span class='py_c_api'>PyFloat_FromDouble</span>(__pyx_v_scale);<span class='error_goto'> if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 44, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_4);\n __pyx_t_3 = <span class='py_c_api'>PyNumber_Multiply</span>(__pyx_t_1, __pyx_t_4);<span class='error_goto'> if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 44, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_3);\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;\n <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_4); __pyx_t_4 = 0;\n if (!(likely(((__pyx_t_3) == Py_None) || likely(<span class='pyx_c_api'>__Pyx_TypeTest</span>(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) <span class='error_goto'>__PYX_ERR(0, 44, __pyx_L1_error)</span>\n __pyx_t_7 = ((PyArrayObject *)__pyx_t_3);\n {\n __Pyx_BufFmt_StackElem __pyx_stack[1];\n if (unlikely(<span class='pyx_c_api'>__Pyx_GetBufferAndValidate</span>(&__pyx_pybuffernd_roots.rcbuffer->pybuffer, (PyObject*)__pyx_t_7, &__Pyx_TypeInfo___pyx_t_double_complex, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 1, 0, __pyx_stack) == -1)) {\n __pyx_v_roots = ((PyArrayObject *)Py_None); <span class='pyx_macro_api'>__Pyx_INCREF</span>(Py_None); __pyx_pybuffernd_roots.rcbuffer->pybuffer.buf = NULL;\n <span class='error_goto'>__PYX_ERR(0, 43, __pyx_L1_error)</span>\n } else {__pyx_pybuffernd_roots.diminfo[0].strides = __pyx_pybuffernd_roots.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_roots.diminfo[0].shape = __pyx_pybuffernd_roots.rcbuffer->pybuffer.shape[0];\n }\n }\n __pyx_t_7 = 0;\n __pyx_v_roots = ((PyArrayObject *)__pyx_t_3);\n __pyx_t_3 = 0;\n</pre><pre class=\"cython line score-0\"> <span class=\"\">45</span>: <span class=\"c\"># np.complex128_t[:] roots = _roots</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">46</span>: <span class=\"n\">np</span><span class=\"o\">.</span><span class=\"n\">complex128_t</span> <span class=\"n\">deno</span><span class=\"p\">,</span> <span class=\"n\">delta</span></pre>\n<pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">47</span>: <span class=\"nb\">long</span> <span class=\"n\">double</span> <span class=\"n\">max_delta</span> <span class=\"o\">=</span> <span class=\"mf\">1</span><span class=\"p\">,</span> <span class=\"n\">d2</span><span class=\"p\">,</span> <span class=\"n\">rd</span><span class=\"p\">,</span> <span class=\"nb\">id</span></pre>\n<pre class='cython code score-0 '> __pyx_v_max_delta = 1.0;\n</pre><pre class=\"cython line score-0\"> <span class=\"\">48</span>: </pre>\n<pre class=\"cython line score-4\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">49</span>: <span class=\"k\">with</span> <span class=\"k\">nogil</span><span class=\"p\">:</span></pre>\n<pre class='cython code score-4 '> {\n #ifdef WITH_THREAD\n PyThreadState *_save;\n Py_UNBLOCK_THREADS\n <span class='pyx_c_api'>__Pyx_FastGIL_Remember</span>();\n #endif\n /*try:*/ {\n/* … */\n /*finally:*/ {\n /*normal exit:*/{\n #ifdef WITH_THREAD\n <span class='pyx_c_api'>__Pyx_FastGIL_Forget</span>();\n Py_BLOCK_THREADS\n #endif\n goto __pyx_L5;\n }\n __pyx_L5:;\n }\n }\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">50</span>: <span class=\"k\">while</span> <span class=\"n\">max_delta</span> <span class=\"o\">></span> <span class=\"n\">epsilon</span> <span class=\"ow\">and</span> <span class=\"n\">ii</span> <span class=\"o\"><</span> <span class=\"n\">max_iter</span><span class=\"p\">:</span></pre>\n<pre class='cython code score-0 '> while (1) {\n __pyx_t_9 = ((__pyx_v_max_delta > __pyx_v_epsilon) != 0);\n if (__pyx_t_9) {\n } else {\n __pyx_t_8 = __pyx_t_9;\n goto __pyx_L8_bool_binop_done;\n }\n __pyx_t_9 = ((__pyx_v_ii < __pyx_v_max_iter) != 0);\n __pyx_t_8 = __pyx_t_9;\n __pyx_L8_bool_binop_done:;\n if (!__pyx_t_8) break;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">51</span>: <span class=\"n\">max_delta</span> <span class=\"o\">=</span> <span class=\"mf\">0</span></pre>\n<pre class='cython code score-0 '> __pyx_v_max_delta = 0.0;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">52</span>: <span class=\"n\">ii</span> <span class=\"o\">+=</span> <span class=\"mf\">1</span></pre>\n<pre class='cython code score-0 '> __pyx_v_ii = (__pyx_v_ii + 1);\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">53</span>: <span class=\"k\">for</span> <span class=\"n\">i</span> <span class=\"ow\">in</span> <span class=\"nb\">range</span><span class=\"p\">(</span><span class=\"n\">L</span><span class=\"p\">):</span></pre>\n<pre class='cython code score-0 '> __pyx_t_10 = __pyx_v_L;\n __pyx_t_11 = __pyx_t_10;\n for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {\n __pyx_v_i = __pyx_t_12;\n</pre><pre class=\"cython line score-0\"> <span class=\"\">54</span>: <span class=\"c\"># deno = 1 + 0j</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">55</span>: <span class=\"c\"># delta = 0 + 0j</span></pre>\n<pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">56</span>: <span class=\"n\">delta</span> <span class=\"o\">=</span> <span class=\"n\">f_c</span><span class=\"p\">(</span><span class=\"n\">poly</span><span class=\"p\">,</span> <span class=\"n\">roots</span><span class=\"p\">[</span><span class=\"n\">i</span><span class=\"p\">])</span></pre>\n<pre class='cython code score-0 '> __pyx_t_13 = __pyx_v_i;\n __pyx_v_delta = __pyx_f_46_cython_magic_7d2d8e83eb391e7da013a05885fe9a7f_f_c(__pyx_v_poly, (*__Pyx_BufPtrStrided1d(__pyx_t_double_complex *, __pyx_pybuffernd_roots.rcbuffer->pybuffer.buf, __pyx_t_13, __pyx_pybuffernd_roots.diminfo[0].strides)));\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">57</span>: <span class=\"k\">for</span> <span class=\"n\">k</span> <span class=\"ow\">in</span> <span class=\"nb\">range</span><span class=\"p\">(</span><span class=\"n\">L</span><span class=\"p\">):</span></pre>\n<pre class='cython code score-0 '> __pyx_t_14 = __pyx_v_L;\n __pyx_t_15 = __pyx_t_14;\n for (__pyx_t_16 = 0; __pyx_t_16 < __pyx_t_15; __pyx_t_16+=1) {\n __pyx_v_k = __pyx_t_16;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">58</span>: <span class=\"k\">if</span> <span class=\"n\">k</span> <span class=\"o\">!=</span> <span class=\"n\">i</span><span class=\"p\">:</span></pre>\n<pre class='cython code score-0 '> __pyx_t_8 = ((__pyx_v_k != __pyx_v_i) != 0);\n if (__pyx_t_8) {\n/* … */\n }\n }\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">59</span>: <span class=\"n\">delta</span> <span class=\"o\">/=</span> <span class=\"p\">(</span><span class=\"n\">roots</span><span class=\"p\">[</span><span class=\"n\">i</span><span class=\"p\">]</span> <span class=\"o\">-</span> <span class=\"n\">roots</span><span class=\"p\">[</span><span class=\"n\">k</span><span class=\"p\">])</span></pre>\n<pre class='cython code score-0 '> __pyx_t_17 = __pyx_v_i;\n __pyx_t_18 = __pyx_v_k;\n __pyx_v_delta = __Pyx_c_quot_double(__pyx_v_delta, __Pyx_c_diff_double((*__Pyx_BufPtrStrided1d(__pyx_t_double_complex *, __pyx_pybuffernd_roots.rcbuffer->pybuffer.buf, __pyx_t_17, __pyx_pybuffernd_roots.diminfo[0].strides)), (*__Pyx_BufPtrStrided1d(__pyx_t_double_complex *, __pyx_pybuffernd_roots.rcbuffer->pybuffer.buf, __pyx_t_18, __pyx_pybuffernd_roots.diminfo[0].strides))));\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">60</span>: <span class=\"n\">roots</span><span class=\"p\">[</span><span class=\"n\">i</span><span class=\"p\">]</span> <span class=\"o\">-=</span> <span class=\"n\">delta</span></pre>\n<pre class='cython code score-0 '> __pyx_t_19 = __pyx_v_i;\n *__Pyx_BufPtrStrided1d(__pyx_t_double_complex *, __pyx_pybuffernd_roots.rcbuffer->pybuffer.buf, __pyx_t_19, __pyx_pybuffernd_roots.diminfo[0].strides) -= __pyx_v_delta;\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">61</span>: <span class=\"n\">rd</span> <span class=\"o\">=</span> <span class=\"n\">creal</span><span class=\"p\">(</span><span class=\"n\">delta</span><span class=\"p\">)</span></pre>\n<pre class='cython code score-0 '> __pyx_v_rd = creal(__pyx_v_delta);\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">62</span>: <span class=\"nb\">id</span> <span class=\"o\">=</span> <span class=\"n\">cimag</span><span class=\"p\">(</span><span class=\"n\">delta</span><span class=\"p\">)</span></pre>\n<pre class='cython code score-0 '> __pyx_v_id = cimag(__pyx_v_delta);\n</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">63</span>: <span class=\"n\">max_delta</span> <span class=\"o\">+=</span> <span class=\"p\">((</span><span class=\"n\">rd</span> <span class=\"o\">*</span> <span class=\"n\">rd</span><span class=\"p\">)</span> <span class=\"o\">+</span> <span class=\"p\">(</span><span class=\"nb\">id</span> <span class=\"o\">*</span> <span class=\"nb\">id</span><span class=\"p\">))</span> <span class=\"o\">/</span> <span class=\"n\">L</span></pre>\n<pre class='cython code score-0 '> __pyx_v_max_delta = (__pyx_v_max_delta + (((__pyx_v_rd * __pyx_v_rd) + (__pyx_v_id * __pyx_v_id)) / ((long double)__pyx_v_L)));\n }\n</pre><pre class=\"cython line score-0\"> <span class=\"\">64</span>: <span class=\"c\"># with gil:</span></pre>\n<pre class=\"cython line score-0\"> <span class=\"\">65</span>: <span class=\"c\"># print(ii, roots, delta)</span></pre>\n<pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">66</span>: <span class=\"n\">max_delta</span> <span class=\"o\">=</span> <span class=\"n\">sqrtl</span><span class=\"p\">(</span><span class=\"n\">max_delta</span><span class=\"p\">)</span></pre>\n<pre class='cython code score-0 '> __pyx_v_max_delta = sqrtl(__pyx_v_max_delta);\n }\n }\n</pre><pre class=\"cython line score-0\"> <span class=\"\">67</span>: <span class=\"c\"># print("max_delta = ", max_delta)</span></pre>\n<pre class=\"cython line score-14\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">68</span>: <span class=\"k\">return</span> <span class=\"n\">roots</span><span class=\"p\">,</span> <span class=\"n\">max_delta</span></pre>\n<pre class='cython code score-14 '> <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_r);\n __pyx_t_3 = <span class='py_c_api'>PyFloat_FromDouble</span>(__pyx_v_max_delta);<span class='error_goto'> if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 68, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_3);\n __pyx_t_4 = <span class='py_c_api'>PyTuple_New</span>(2);<span class='error_goto'> if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 68, __pyx_L1_error)</span>\n <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_4);\n <span class='pyx_macro_api'>__Pyx_INCREF</span>(((PyObject *)__pyx_v_roots));\n <span class='refnanny'>__Pyx_GIVEREF</span>(((PyObject *)__pyx_v_roots));\n <span class='py_macro_api'>PyTuple_SET_ITEM</span>(__pyx_t_4, 0, ((PyObject *)__pyx_v_roots));\n <span class='refnanny'>__Pyx_GIVEREF</span>(__pyx_t_3);\n <span class='py_macro_api'>PyTuple_SET_ITEM</span>(__pyx_t_4, 1, __pyx_t_3);\n __pyx_t_3 = 0;\n __pyx_r = __pyx_t_4;\n __pyx_t_4 = 0;\n goto __pyx_L0;\n</pre></div></body></html>" | |
| }, | |
| "metadata": {} | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "# Testing & Benchmarking" | |
| }, | |
| { | |
| "metadata": { | |
| "ExecuteTime": { | |
| "end_time": "2020-02-29T22:12:11.494565Z", | |
| "start_time": "2020-02-29T22:12:11.466458Z" | |
| }, | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "def benchmark(polynom):\n L = len(polynom)\n msg = ''\n for i in range(L):\n if polynom[i] == 0:\n continue\n if i > 0:\n msg += '+ '\n msg += '({})x^{} '.format(polynom[i], L-i-1)\n print(msg)\n print('--------------')\n \n print(\"Using numpy.solve\")\n %timeit np.roots(polynom)\n roots1 = np.roots(polynom)\n print(roots1)\n print('--------------')\n print(\"Using numba implementation\")\n %timeit roots = solve(polynom)\n roots = solve(polynom)\n print(roots)\n print('--------------')\n print(\"Using cython implementation\")\n %timeit roots2, _ = solve_c(polynom)\n roots2, _ = solve_c(polynom)\n print(roots2)\n print('--------------')", | |
| "execution_count": 10, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "## Known bugs in Durand–Kerner algorithm\n### Fail cases" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# 1 - fail example:\npolynom = np.array([4, 0, -1], dtype=float)\nbenchmark(polynom)\nprint('================')\n# 2 - fail example\npolynom = np.array([2, 13, -7], dtype=float)\nbenchmark(polynom)", | |
| "execution_count": 14, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "text": "(4.0)x^2 + (-1.0)x^0 \n--------------\nUsing numpy.solve\n76.2 µs ± 5.7 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n[ 0.5 -0.5]\n--------------\nUsing numba implementation\n122 µs ± 39.5 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n(array([-inf+infj, nan+nanj]), nan)\n--------------\nUsing cython implementation\n65.7 µs ± 488 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n[-90.81607955-24.37936125j inf +nanj]\n--------------\n(2.0)x^2 + (13.0)x^1 + (-7.0)x^0 \n--------------\nUsing numpy.solve\n74.1 µs ± 7.7 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n[-7. 0.5]\n--------------\nUsing numba implementation\n11.1 ms ± 167 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n(array([ -10.15437377 -6.61299036j, 4040.40420859+4516.35012194j]), 8559.543666648191)\n--------------\nUsing cython implementation\n190 ms ± 1.86 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n[-6.11547323e+04-7.82040356e+04j 2.73825269e+01+3.89510547e+01j]\n--------------\n", | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### Fail cases resolved" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# 1 - fail example - resolved:\npolynom = np.array([4, 0, -1])\npolynom = polynom / polynom[0]\nbenchmark(polynom)\n# 2 - fail example - resolved\npolynom = np.array([2, 13, -7])\npolynom = polynom / polynom[0]\nbenchmark(polynom)", | |
| "execution_count": 18, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "text": "(1.0)x^2 + (-0.25)x^0 \n--------------\nUsing numpy.solve\n98.8 µs ± 24.2 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n[ 0.5 -0.5]\n--------------\nUsing numba implementation\n71.7 µs ± 5.13 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n(array([ 0.5+0.j, -0.5+0.j]), 3.1526380599683343e-22)\n--------------\nUsing cython implementation\n13.9 µs ± 89.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n[ 0.5+0.00000000e+00j -0.5+6.34935309e-18j]\n--------------\n(1.0)x^2 + (6.5)x^1 + (-3.5)x^0 \n--------------\nUsing numpy.solve\n71.4 µs ± 3.71 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n[-7. 0.5]\n--------------\nUsing numba implementation\n63.7 µs ± 4.66 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n(array([ 0.5+0.00000000e+00j, -7. +4.89163707e-16j]), 1.3882670517149822e-16)\n--------------\nUsing cython implementation\n15.4 µs ± 1.7 µs per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n[-7. -4.62060783e-16j 0.5-3.85185989e-33j]\n--------------\n", | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": { | |
| "ExecuteTime": { | |
| "end_time": "2020-02-29T22:12:22.905023Z", | |
| "start_time": "2020-02-29T22:12:11.526177Z" | |
| }, | |
| "scrolled": false, | |
| "trusted": false | |
| }, | |
| "cell_type": "code", | |
| "source": "polynom = np.array([1, 0, 0, 10, 0, 0, -27])\nbenchmark(polynom)", | |
| "execution_count": 5, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "text": "(1)x^6 + (10)x^3 + (-27)x^0 \n--------------\nUsing numpy.solve\n113 µs ± 27.8 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n[-2.30277564+0.j 1.15138782+1.9942622j 1.15138782-1.9942622j\n -0.65138782+1.1282368j -0.65138782-1.1282368j 1.30277564+0.j ]\n--------------\nUsing numba implementation\nThe slowest run took 577.42 times longer than the fastest. This could mean that an intermediate result is being cached.\n20.7 ms ± 12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n(array([ 1.15138782-1.99426220e+00j, -0.65138782-1.12823680e+00j,\n -2.30277564-1.27822331e-16j, 1.30277564+1.02494282e-32j,\n -0.65138782+1.12823680e+00j, 1.15138782+1.99426220e+00j]), 2.2283351559332415e-16)\n--------------\nUsing cython implementation\n246 ms ± 83.8 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n[-0.65138782-1.12823680e+00j 1.30277564+5.02605867e-32j\n -0.65138782+1.12823680e+00j 1.15138782-1.99426220e+00j\n 1.15138782+1.99426220e+00j -2.30277564-7.30758745e-16j]\n--------------\n", | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": { | |
| "ExecuteTime": { | |
| "start_time": "2020-02-29T22:11:50.565Z" | |
| }, | |
| "trusted": false | |
| }, | |
| "cell_type": "code", | |
| "source": "polynom = np.array([1,0,0,0,0,0,0,0,0,-1])\nbenchmark(polynom)", | |
| "execution_count": 6, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "text": "(1)x^9 + (-1)x^0 \n--------------\nUsing numpy.solve\n225 µs ± 89.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n[-0.93969262+0.34202014j -0.93969262-0.34202014j -0.5 +0.8660254j\n -0.5 -0.8660254j 0.17364818+0.98480775j 0.17364818-0.98480775j\n 1. +0.j 0.76604444+0.64278761j 0.76604444-0.64278761j]\n--------------\nUsing numba implementation\n91.3 µs ± 2.22 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n(array([ 1. +1.50463277e-36j, -0.5 -8.66025404e-01j,\n 0.17364818+9.84807753e-01j, 0.76604444-6.42787610e-01j,\n 0.76604444+6.42787610e-01j, -0.93969262-3.42020143e-01j,\n -0.93969262+3.42020143e-01j, -0.5 +8.66025404e-01j,\n 0.17364818-9.84807753e-01j]), 1.6244362503173417e-16)\n--------------\nUsing cython implementation\n88.8 µs ± 7.2 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n[-0.93969262+3.42020143e-01j -0.93969262-3.42020143e-01j\n 0.17364818-9.84807753e-01j 0.76604444-6.42787610e-01j\n 1. -2.61012179e-54j 0.17364818+9.84807753e-01j\n 0.76604444+6.42787610e-01j -0.5 -8.66025404e-01j\n -0.5 +8.66025404e-01j]\n--------------\n", | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": { | |
| "ExecuteTime": { | |
| "start_time": "2020-02-29T22:11:50.593Z" | |
| }, | |
| "trusted": false | |
| }, | |
| "cell_type": "code", | |
| "source": "polynom = np.array([1, 0, 0, 0, -1, -1])\nbenchmark(polynom)", | |
| "execution_count": 7, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "text": "(1)x^5 + (-1)x^1 + (-1)x^0 \n--------------\nUsing numpy.solve\n79.6 µs ± 6.29 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n[ 1.16730398+0.j 0.18123244+1.0839541j 0.18123244-1.0839541j\n -0.76488443+0.35247155j -0.76488443-0.35247155j]\n--------------\nUsing numba implementation\n44.7 µs ± 896 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n(array([-0.76488443+3.52471546e-01j, 1.16730398-7.24892101e-33j,\n -0.76488443-3.52471546e-01j, 0.18123244+1.08395410e+00j,\n 0.18123244-1.08395410e+00j]), 1.56743240643403e-16)\n--------------\nUsing cython implementation\n40.7 µs ± 286 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n[ 1.16730398-3.65926689e-33j -0.76488443+3.52471546e-01j\n -0.76488443-3.52471546e-01j 0.18123244-1.08395410e+00j\n 0.18123244+1.08395410e+00j]\n--------------\n", | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": { | |
| "ExecuteTime": { | |
| "start_time": "2020-02-29T22:11:50.665Z" | |
| }, | |
| "trusted": false | |
| }, | |
| "cell_type": "code", | |
| "source": "polynom = np.array([1,-3,3,-5])\nbenchmark(polynom)", | |
| "execution_count": 8, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "text": "(1)x^3 + (-3)x^2 + (3)x^1 + (-5)x^0 \n--------------\nUsing numpy.solve\n74 µs ± 4.73 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n[2.58740105+0.j 0.20629947+1.37472964j 0.20629947-1.37472964j]\n--------------\nUsing numba implementation\n20 µs ± 1.84 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n(array([0.20629947-1.37472964e+00j, 0.20629947+1.37472964e+00j,\n 2.58740105+5.03794111e-32j]), 3.581396832712688e-16)\n--------------\nUsing cython implementation\n27 µs ± 1.66 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n[2.58740105+1.24235219e-32j 0.20629947-1.37472964e+00j\n 0.20629947+1.37472964e+00j]\n--------------\n", | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": { | |
| "scrolled": false, | |
| "trusted": false | |
| }, | |
| "cell_type": "code", | |
| "source": "N = 61\npolynom = np.zeros(N, dtype=int)\npolynom[0] = 1\npolynom[N-1] = -1\nbenchmark(polynom)", | |
| "execution_count": 9, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "text": "(1)x^60 + (-1)x^0 \n--------------\nUsing numpy.solve\n951 µs ± 95.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n[-1.00000000e+00+0.j -9.94521895e-01+0.10452846j\n -9.94521895e-01-0.10452846j -9.78147601e-01+0.20791169j\n -9.78147601e-01-0.20791169j -9.51056516e-01+0.30901699j\n -9.51056516e-01-0.30901699j -9.13545458e-01+0.40673664j\n -9.13545458e-01-0.40673664j -8.66025404e-01+0.5j\n -8.66025404e-01-0.5j -8.09016994e-01+0.58778525j\n -8.09016994e-01-0.58778525j -7.43144825e-01+0.66913061j\n -7.43144825e-01-0.66913061j -6.69130606e-01+0.74314483j\n -6.69130606e-01-0.74314483j -5.87785252e-01+0.80901699j\n -5.87785252e-01-0.80901699j -5.00000000e-01+0.8660254j\n -5.00000000e-01-0.8660254j -4.06736643e-01+0.91354546j\n -4.06736643e-01-0.91354546j -3.09016994e-01+0.95105652j\n -3.09016994e-01-0.95105652j -2.07911691e-01+0.9781476j\n -2.07911691e-01-0.9781476j -1.04528463e-01+0.9945219j\n -1.04528463e-01-0.9945219j -5.68989300e-16+1.j\n -5.68989300e-16-1.j 1.04528463e-01+0.9945219j\n 1.04528463e-01-0.9945219j 2.07911691e-01+0.9781476j\n 2.07911691e-01-0.9781476j 3.09016994e-01+0.95105652j\n 3.09016994e-01-0.95105652j 4.06736643e-01+0.91354546j\n 4.06736643e-01-0.91354546j 5.00000000e-01+0.8660254j\n 5.00000000e-01-0.8660254j 5.87785252e-01+0.80901699j\n 5.87785252e-01-0.80901699j 6.69130606e-01+0.74314483j\n 6.69130606e-01-0.74314483j 7.43144825e-01+0.66913061j\n 7.43144825e-01-0.66913061j 8.09016994e-01+0.58778525j\n 8.09016994e-01-0.58778525j 1.00000000e+00+0.j\n 9.94521895e-01+0.10452846j 9.94521895e-01-0.10452846j\n 9.78147601e-01+0.20791169j 9.78147601e-01-0.20791169j\n 9.51056516e-01+0.30901699j 9.51056516e-01-0.30901699j\n 9.13545458e-01+0.40673664j 9.13545458e-01-0.40673664j\n 8.66025404e-01+0.5j 8.66025404e-01-0.5j ]\n--------------\nUsing numba implementation\n4.64 ms ± 330 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n(array([ 4.06736643e-01+9.13545458e-01j, 9.51056516e-01+3.09016994e-01j,\n 9.78147601e-01-2.07911691e-01j, -5.87785252e-01-8.09016994e-01j,\n 6.69130606e-01+7.43144825e-01j, 8.09016994e-01+5.87785252e-01j,\n 5.87785252e-01-8.09016994e-01j, -3.09016994e-01-9.51056516e-01j,\n -1.00000000e+00+1.19426323e-16j, 8.66025404e-01-5.00000000e-01j,\n -2.07911691e-01-9.78147601e-01j, 1.00000000e+00+1.06910588e-50j,\n 9.94521895e-01+1.04528463e-01j, -2.77756903e-16-1.00000000e+00j,\n -9.94521895e-01-1.04528463e-01j, -5.00000000e-01-8.66025404e-01j,\n -6.78056107e-17+1.00000000e+00j, 9.94521895e-01-1.04528463e-01j,\n -1.04528463e-01+9.94521895e-01j, 8.66025404e-01+5.00000000e-01j,\n 1.04528463e-01-9.94521895e-01j, 7.43144825e-01+6.69130606e-01j,\n -7.43144825e-01-6.69130606e-01j, 6.69130606e-01-7.43144825e-01j,\n -9.94521895e-01+1.04528463e-01j, 5.87785252e-01+8.09016994e-01j,\n -3.09016994e-01+9.51056516e-01j, -9.51056516e-01+3.09016994e-01j,\n 9.13545458e-01+4.06736643e-01j, 2.07911691e-01-9.78147601e-01j,\n -9.78147601e-01+2.07911691e-01j, 1.04528463e-01+9.94521895e-01j,\n -9.78147601e-01-2.07911691e-01j, 9.13545458e-01-4.06736643e-01j,\n -8.66025404e-01-5.00000000e-01j, -8.09016994e-01-5.87785252e-01j,\n -6.69130606e-01+7.43144825e-01j, -4.06736643e-01-9.13545458e-01j,\n -8.66025404e-01+5.00000000e-01j, -9.13545458e-01-4.06736643e-01j,\n 9.51056516e-01-3.09016994e-01j, -5.00000000e-01+8.66025404e-01j,\n -9.13545458e-01+4.06736643e-01j, -9.51056516e-01-3.09016994e-01j,\n 8.09016994e-01-5.87785252e-01j, 3.09016994e-01+9.51056516e-01j,\n 5.00000000e-01-8.66025404e-01j, -4.06736643e-01+9.13545458e-01j,\n 7.43144825e-01-6.69130606e-01j, -8.09016994e-01+5.87785252e-01j,\n -6.69130606e-01-7.43144825e-01j, 4.06736643e-01-9.13545458e-01j,\n 2.07911691e-01+9.78147601e-01j, 9.78147601e-01+2.07911691e-01j,\n -2.07911691e-01+9.78147601e-01j, -7.43144825e-01+6.69130606e-01j,\n -5.87785252e-01+8.09016994e-01j, -1.04528463e-01-9.94521895e-01j,\n 5.00000000e-01+8.66025404e-01j, 3.09016994e-01-9.51056516e-01j]), 1.5530898495486913e-16)\n--------------\nUsing cython implementation\n6.3 ms ± 352 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n[ 5.87785252e-01-8.09016994e-01j -5.00000000e-01+8.66025404e-01j\n -5.87785252e-01+8.09016994e-01j 7.43144825e-01-6.69130606e-01j\n 3.09016994e-01-9.51056516e-01j -3.09016994e-01+9.51056516e-01j\n -8.66025404e-01+5.00000000e-01j -2.07911691e-01-9.78147601e-01j\n 6.69130606e-01-7.43144825e-01j -5.00000000e-01-8.66025404e-01j\n -1.04528463e-01-9.94521895e-01j -8.66025404e-01-5.00000000e-01j\n -4.06736643e-01+9.13545458e-01j -1.04528463e-01+9.94521895e-01j\n 5.87785252e-01+8.09016994e-01j -8.09016994e-01+5.87785252e-01j\n 1.04528463e-01-9.94521895e-01j -1.26407118e-16+1.00000000e+00j\n -9.51056516e-01+3.09016994e-01j -8.09016994e-01-5.87785252e-01j\n 9.13545458e-01-4.06736643e-01j 1.00000000e+00-2.23032868e-57j\n 8.66025404e-01-5.00000000e-01j -6.69130606e-01+7.43144825e-01j\n 6.69130606e-01+7.43144825e-01j 1.04528463e-01+9.94521895e-01j\n 9.51056516e-01-3.09016994e-01j 9.78147601e-01+2.07911691e-01j\n 7.43144825e-01+6.69130606e-01j -6.69130606e-01-7.43144825e-01j\n 5.00000000e-01-8.66025404e-01j 3.09016994e-01+9.51056516e-01j\n -9.78147601e-01-2.07911691e-01j 8.66025404e-01+5.00000000e-01j\n -3.09016994e-01-9.51056516e-01j 4.06736643e-01-9.13545458e-01j\n -9.51263465e-17-1.00000000e+00j 2.07911691e-01-9.78147601e-01j\n -5.87785252e-01-8.09016994e-01j -4.06736643e-01-9.13545458e-01j\n -9.51056516e-01-3.09016994e-01j -9.78147601e-01+2.07911691e-01j\n 9.13545458e-01+4.06736643e-01j -9.13545458e-01+4.06736643e-01j\n -7.43144825e-01+6.69130606e-01j 9.94521895e-01+1.04528463e-01j\n 5.00000000e-01+8.66025404e-01j 8.09016994e-01+5.87785252e-01j\n -9.94521895e-01+1.04528463e-01j -1.00000000e+00-1.46665411e-16j\n 4.06736643e-01+9.13545458e-01j 9.94521895e-01-1.04528463e-01j\n 9.78147601e-01-2.07911691e-01j 8.09016994e-01-5.87785252e-01j\n -9.13545458e-01-4.06736643e-01j -9.94521895e-01-1.04528463e-01j\n -2.07911691e-01+9.78147601e-01j 2.07911691e-01+9.78147601e-01j\n -7.43144825e-01-6.69130606e-01j 9.51056516e-01+3.09016994e-01j]\n--------------\n", | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": false | |
| }, | |
| "cell_type": "code", | |
| "source": "", | |
| "execution_count": null, | |
| "outputs": [] | |
| } | |
| ], | |
| "metadata": { | |
| "_draft": { | |
| "nbviewer_url": "https://gist.github.com/9713842d5e264c94e293a930b479f597" | |
| }, | |
| "gist": { | |
| "id": "9713842d5e264c94e293a930b479f597", | |
| "data": { | |
| "description": "make the cofficient of x^n equal to 1. This improves the stability of the algorithm", | |
| "public": true | |
| } | |
| }, | |
| "hide_input": false, | |
| "kernelspec": { | |
| "name": "python3", | |
| "display_name": "Python 3", | |
| "language": "python" | |
| }, | |
| "nbTranslate": { | |
| "hotkey": "alt-t", | |
| "sourceLang": "en", | |
| "targetLang": "fr", | |
| "displayLangs": [ | |
| "*" | |
| ], | |
| "langInMainMenu": true, | |
| "useGoogleTranslate": true | |
| }, | |
| "toc": { | |
| "nav_menu": {}, | |
| "number_sections": true, | |
| "sideBar": true, | |
| "skip_h1_title": false, | |
| "base_numbering": 1, | |
| "title_cell": "Table of Contents", | |
| "title_sidebar": "Contents", | |
| "toc_cell": true, | |
| "toc_position": {}, | |
| "toc_section_display": true, | |
| "toc_window_display": false | |
| }, | |
| "varInspector": { | |
| "window_display": true, | |
| "cols": { | |
| "lenName": 16, | |
| "lenType": 16, | |
| "lenVar": 40 | |
| }, | |
| "kernels_config": { | |
| "python": { | |
| "library": "var_list.py", | |
| "delete_cmd_prefix": "del ", | |
| "delete_cmd_postfix": "", | |
| "varRefreshCmd": "print(var_dic_list())" | |
| }, | |
| "r": { | |
| "library": "var_list.r", | |
| "delete_cmd_prefix": "rm(", | |
| "delete_cmd_postfix": ") ", | |
| "varRefreshCmd": "cat(var_dic_list()) " | |
| } | |
| }, | |
| "types_to_exclude": [ | |
| "module", | |
| "function", | |
| "builtin_function_or_method", | |
| "instance", | |
| "_Feature" | |
| ], | |
| "position": { | |
| "left": "1228px", | |
| "top": "197px", | |
| "width": "312px", | |
| "height": "384px", | |
| "right": "20px" | |
| } | |
| }, | |
| "language_info": { | |
| "name": "python", | |
| "version": "3.7.7", | |
| "mimetype": "text/x-python", | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 3 | |
| }, | |
| "pygments_lexer": "ipython3", | |
| "nbconvert_exporter": "python", | |
| "file_extension": ".py" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 4 | |
| } |
Author
Thank you so much for doing this! I was a bit lost after reading the Wikipedia article, but this write-up made it so much clearer. I would also like to propose an alternative implementation that, even if not 100% equivalent, makes the computation a lot more SIMD-friendly (the actual SIMD implementation is missing, though):
import numpy as np
def durand_kerner(
p: np.ndarray,
max_iter: int = 1000,
tol: float = 1e-8
) -> np.ndarray:
# Polynomial scaling (roots remain unchanged)
p = p / p[0]
# Deterministic init using the complex unit circle
deg = len(p) - 1
roots = np.exp(np.arange(deg) * 2j * np.pi / deg)
f0, cnt = np.polyval(p, roots), 0
avg_delta = tol + 1 # do-while enforcing
while avg_delta > tol and cnt < max_iter:
cnt += 1
# Simultaneous distance calculation using advanced indexing
dens = np.prod(roots[:, None] - roots[None, :] + np.eye(deg), axis=1)
roots = roots - f0 / dens
f1 = np.polyval(p, roots)
# 2-norm of the change in polyval
avg_delta = np.sqrt(np.sum((f1 - f0)**2) / deg)
f0 = f1
return roots
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Performance improvement
Performing
i.e., by ensuing the coefficient of the highest power of
xto be1significantly improves the numerical stability.For example, the above algorithm fails for:
but the problem gets resolved when executing as such: