def parallel_newton (func, x0, fprime=None, par_args=(), simple_args=(), tol=1.48e-8,
                     maxiter=50, parallel=True, **kwargs):
    """A parallelized version of :func:`scipy.optimize.newton`.

    Arguments:

    func
      The function to search for zeros, called as ``f(x, [*par_args...], [*simple_args...])``.
    x0
      The initial point for the zero search.
    fprime
      (Optional) The first derivative of *func*, called the same way.
    par_args
      Tuple of additional parallelized arguments.
    simple_args
      Tuple of additional arguments passed identically to every invocation.
    tol
      The allowable error of the zero value.
    maxiter
      Maximum number of iterations.
    parallel
      Controls parallelization; default uses all available cores. See
      :func:`pwkit.parallel.make_parallel_helper`.
    kwargs
      Passed to :func:`scipy.optimize.newton`.

    Returns: an array of locations of zeros.

    Finds zeros in parallel. The values *x0*, *tol*, *maxiter*, and the items
    of *par_args* should all be numeric, and may be N-dimensional Numpy
    arrays. They are all broadcast to a common shape, and one zero-finding run
    is performed for each element in the resulting array. The return value is
    an array of zero locations having the same shape as the common broadcast
    of the parameters named above.

    The *simple_args* are passed to each function identically for each
    integration. They do not need to be Pickle-able.

    Example::

       >>> parallel_newton (lambda x, a: x - 2 * a, 2,
                            par_args=(np.arange (6),))
       <<< array([  0.,   2.,   4.,   6.,   8.,  10.])
       >>> parallel_newton (lambda x: np.sin (x), np.arange (6))
       <<< array([  0.00000000e+00,   3.65526589e-26,   3.14159265e+00,
                    3.14159265e+00,   3.14159265e+00,   6.28318531e+00])

    """
    from scipy.optimize import newton

    from .parallel import make_parallel_helper
    phelp = make_parallel_helper (parallel)

    if not isinstance (par_args, tuple):
        raise ValueError ('par_args must be a tuple')

    if not isinstance (simple_args, tuple):
        raise ValueError ('simple_args must be a tuple')

    bc_raw = np.broadcast_arrays (x0, tol, maxiter, *par_args)
    bc_1d = tuple (np.atleast_1d (a) for a in bc_raw)

    def gen_var_args ():
        for i in range (bc_1d[0].size):
            yield tuple (x.flat[i] for x in bc_1d)

    def helper (i, _, var_args):
        x0, tol, maxiter = var_args[:3]
        args = var_args[3:] + simple_args
        return newton (func, x0, fprime=fprime, args=args, tol=tol,
                       maxiter=maxiter, **kwargs)

    with phelp.get_ppmap () as ppmap:
        result = np.asarray (ppmap (helper, None, gen_var_args ()))

    if bc_raw[0].ndim == 0:
        return np.asscalar (result)
    return result