def halley_newton ( fun            ,    ## the function 
                    x              ,    ## x 
                    deriv1         ,    ## the first derivative 
                    deriv2  = None ,    ## the second derivative
                    fx      = None ,    ## value of fun(x)
                    args    = ()   ) :  ## additional arguments for function calls 
    """Single step of Newton/Halley's algorithm

    Parameters
    -----------
    
    fun    : the function 
    x      : the initial guess for the root positon
    deriv  : the first derivative  (for Newton and Halley's method) 
    deriv2 : the second derivative (for Halley's method)
    fx     : (optional) the value of funnction at the guess point
    args   : (optional) parameters for function/derivative calls 

    Returns
    -------
    The next approximation to the root or `None`
    
    Example
    -------
    
    >>> fun   = lambda x : sin ( x )
    >>> deriv = lambda x : cos ( x )
    >>> x = 0.5
    >>> for i in  range ( 10 ) :
    ...      x = halley_newton ( fun , x , deriv )
    ...      print i, x

    References
    ----------
    - see https://en.wikipedia.org/wiki/Halley%27s_method
    - see https://en.wikipedia.org/wiki/Newton%27s_method     
   
    """
    
    ## newton corrections
    d1 = float ( deriv1 ( x , *args ) ) 
    fx = float ( fun    ( x , *args ) ) if fx is None else fx 
    
    if d1  : rn = fx / d1
    else   : return None                               ## error here! 
    
    ## make corrections: Halley's  steps
    if deriv2 :  
        d2 = float ( deriv2 ( x , *args ) )
        if d2 : rn /= ( 1.0  - 0.5 * rn * d2 / d1 ) ## Halley's correction 
        
    return x - rn