#!/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.
#
 
from numpy import linspace
 
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 ()
    """
 
 
    # Creiamo gli indici su cui iterare
    x_values = linspace(-2,2,size,True)
    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. Salviamola prima però
    # in una stringa. Sì, occupa un sacco di RAM ma tanto
    # quella in genere non manca, e le scritture ripetute
    # sono esose di risorse.
    fileContent = ""
    for x in x_values:
        for y in y_values:
            value = GetNewtonConvergenceSpeed(complex(x,y))
            fileContent += " ".join([str(3*value), str(5*value), str(9*value)]) + "\n"
    f.write(fileContent)
    f.close ()
 
 
if __name__ == "__main__":
 
    Newton (500)