From b80ca5edcb014e9bbd2c8c8ac822f7846a31144a Mon Sep 17 00:00:00 2001
From: Leonardo Robol <leo@robol.it>
Date: Tue, 23 Feb 2010 08:56:33 +0100
Subject: [PATCH] =?UTF-8?q?Versione=20che=20usa=20numpy=20(solo=20pi=C3=B9?=
 =?UTF-8?q?=20pulito=20il=20codice,=20non=20c'=C3=A8=20nesssun=20vantaggio?=
 =?UTF-8?q?=20in=20termini=20di=20prestazioni=20per=20ora).?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 NewtonFractalNumpy.py | 78 +++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 78 insertions(+)
 create mode 100644 NewtonFractalNumpy.py

diff --git a/NewtonFractalNumpy.py b/NewtonFractalNumpy.py
new file mode 100644
index 0000000..65a45be
--- /dev/null
+++ b/NewtonFractalNumpy.py
@@ -0,0 +1,78 @@
+#!/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)
-- 
2.1.4