#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Questo script genera il frattale di newton
# associato al polinomio x^alpha - 1; Può essere
# divertente cambiare alpha ed eventualemnte
# anche il numero di pixel dell'immagine.
#
 
def GetNewtonConvergenceSpeed(z, maxit = 255, eps = 10e-11, alpha = 3):
    """
    Ritorna, dato un numero complesso z ed un eps > 0, il numero
    di iterazioni necessarie (del metodo di Newton) perché
    |f(x_k)| < eps, oppure maxit se dopo maxit iterazioni
    non si è ancora raggiunta la precisione desiderata.
    """
 
    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):
    """
    Calcola la velocità di convergenza su una rete di punti
    size x size sul quadrato di [-2,2]^2 attorno all'origine
    del piano complesso. La salva poi in un'immagine dal nome
    newton.pnm
    """
 
    # Delta è l'ampiezza dell'intervallo fratto il numero di
    # divisioni - 1
    delta = 4 / float(size - 1)
 
    # Creiamo gli indici su cui iterare
    x_range = xrange(0 , size)
    y_range = x_range
 
 
    # 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_range:
        for y in y_range:
            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__":
 
    # Creiamo un'immagine 1024px x 1024px.
    Newton (1024)