#!/usr/bin/enb python
# -*- coding: utf-8 -*-
#
# Sample script to generate a PNM image
# of the Newton's fractal associated with
# x^\alpha - 1 polynomial.
#
 
def GetNewtonConvergenceSpeed(z, maxit = 255, eps = 10e-11, alpha = 3):
    """
    This function returns, given a complex number z,
    an integer between 0 and maxit representing the
    number of iteration to approximate a root of
    the polynomial x^\alpha - 1 = f(x) so that
    |f(x)| < eps.
    """
 
    fz = pow(z, alpha) - 1
    iterations = 0
 
    while (abs(fz) > eps and iterations < maxit):
 
        # Usiamo una var temporanea per non dover
        # ricalcolare le parti di comuni di funzione
        # e derivata.
        t = pow(z, alpha - 1)
        fz = t * z - 1
        fpz = alpha * t
 
 
        # Questo per prevenire la divisione per zero
        # (tanto la nostra derivata si annulla solo lì)
        if (fpz == 0):
            return maxit
 
        z = z - fz / fpz
        iterations += 1
 
    return iterations
 
def Newton(size = 200):
    """Compute the newton's method on a net of points
    in the set { x + iy \in C | |x| < 2 and |y| < 2 }
    and fill the image matrix with the value returned
    by the GetNewtonConvergenceSpeed ()
    """
 
    # Delta è l'ampiezza dell'intervallo fratto il numero di
    # divisioni - 1
    delta = 4 / float(size - 1)
 
    # Creiamo gli indici su cui iterare
    x_values = xrange(0 , size)
    y_values = x_values
 
 
    # Apriamo il file dove salveremo il nostro lavoro
    # e scriviamo l'intestazione del file PNM.
    f = open("newton.pnm", 'w')
    f.write("P3\n")
    f.write(str(size) + " " + str(size) + "\n")
    f.write("255\n")
 
    # Scriviamo la matrice su file
    for x in x_values:
        for y in y_values:
            value = GetNewtonConvergenceSpeed(complex(-2 + x * delta,
                                                      -2 + y * delta))
            f.write( str(3 * value) + " " +
                     str(5 * value) + " " +
                     str(9 * value) + "\n")
    f.close ()
 
 
if __name__ == "__main__":
 
    Newton (500)