title: Vectorize a Newton method in Python/Numpy
body: <p>I am trying to figure out if Python/Numpy is a viable alternative to develop my numerical software which is already available in C++. In order to get performance in Python/Numpy, one need to "vectorize" the code. But it turns out that as soon as I move away from very simple examples, I struggle to vectorize the code (I am not talking about SIMD instructions but "efficient Numpy code" without loops). Here is an algorithm that I want to get efficiently in Python/Numpy.</p>
  
  <ol>
  <li>Create an numpy array containing: 1.0, 1.0 + 1/n, 1.0 + 2/n, ..., 2.0</li>
  <li>For every u in the array, compute the root of x^2 - u, using a Newton method, stopping when |dx| &lt;= 1.0e-7. Store the result in an array result.</li>
  <li>Sum all the elements of the result array</li>
  </ol>
  
  <p>Here is the algorithm in Python I want to speed up</p>
  
  <pre><code>import numpy as np
  
  n = 1000000
  data = np.arange(1.0, 2.0, 1.0 / n)
  
  def newton(u):
    x = 2.0
    while True:
      f = x**2 - u
      df_dx = 2 * x
      dx = f / df_dx
      if (abs(dx) &lt;= 1.0e-7):
        break
      x -= dx
    return x
  
    result = map(newton, data)
  
    print result[n - 1]
  </code></pre>
  
  <p>Here is a version of the algorithm in C++11</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;vector&gt;
  #include &lt;cmath&gt;
  
  int main (int argc, char const *argv[]) {
    auto n = std::size_t{100000000};
  
    auto v = std::vector&lt;double&gt;(n + 1);
    for(size_t k = 0; k &lt; v.size(); ++k) {
      v[k] = 1.0 + static_cast&lt;double&gt;(k) / n;
    }
  
    auto result = std::vector&lt;double&gt;(n + 1);
    for(size_t k = 0; k &lt; v.size(); ++k) {
      auto x = double{2.0};
      while(true) {
        auto f = double{x * x - v[k]};
        auto df_dx = double{2 * x};
        auto dx = double{f / df_dx};
        if (std::abs(dx) &lt;= 1.0e-7) {
          break;
        }
        x -= dx;
      }
      result[k] = x;
    }
  
    auto somme = double{0.0};
    for(size_t k = 0; k &lt; result.size(); ++k) {
      somme += result[k];
    }
  
    std::cout &lt;&lt; somme &lt;&lt; std::endl;
    return 0;
  }
  </code></pre>
  
  <p>It takes 2.9 seconds to run on my machine. Is there a way to make a fast Python/Numpy algorithm that does the same thing (I am willing to get something that is less than 5 times slower).</p>
  
  <p>Thanks.</p>