Vous connaissez tous la méthode de Newton pour approximer les racines d'une fonction, n'est-ce pas? Mon objectif dans cette tâche est de vous présenter un aspect intéressant de cet algorithme.

L'algorithme de Newton ne converge que pour certaines, mais surtout pour les valeurs d'entrée complexes. Si vous imaginez la convergence de la méthode pour toutes les valeurs d'entrée sur le plan complexe, vous obtenez généralement une belle fractale comme celle-ci:

Image de Wikimedia commons

Caractéristiques

Le but de cette tâche est de générer de telles fractales. Cela signifie que vous obtenez un polynôme en entrée et que vous devez imprimer la fractale correspondante en tant qu'image dans un format de votre choix en tant que sortie.

Contribution

L'entrée est une liste de nombres complexes séparés par des espaces. Ils sont écrits dans le style <Real part><iImaginary part>, comme ce numéro: 5.32i3.05. Vous pouvez supposer que le nombre d'entrée n'a pas plus de 4 décimales et est inférieur à 1000. Le premier d'entre eux ne doit pas être nul. Par exemple, cela pourrait être une entrée pour votre programme:

1 -2i7,5 23,0004i-3,8 i12 0 5,1233i0,1

Les nombres sont interprétés comme les coefficients d'un polynôme, en commençant par la puissance la plus élevée. Dans le reste de cette description, le polynôme d'entrée est appelé P . L'entrée ci-dessus est égale à ce polynôme:

f (x) = x 5 + (-2 + 7,5 i ) x 4 + (23,0004 - 3,8 i ) x 3 + 12 i x 2 + 5,1233 + 0,1 i

L'entrée peut vous venir soit de stdin, soit d'un argument passé au programme, soit d'une invite affichée à votre programme. Vous pouvez supposer que l'entrée ne contient aucun espace blanc de début ou de fin.

Le rendu

Vous devez rendre la fractale de la manière suivante:

• Choisissez autant de couleurs que de racines de P plus une couleur supplémentaire pour la divergence
• Pour chaque nombre dans le plan visible, déterminez si la méthode converge et si oui vers quelle racine. Colorez le point en fonction du résultat.
• N'imprimez pas de règles ou d'autres objets de fantaisie
• Imprimez un point noir aux points, qui sont les racines des polynômes d'orientation. Vous pouvez imprimer jusqu'à quatre pixels autour de chaque racine.
• Trouvez un moyen de choisir le plan visible de manière à ce que toutes les racines se distinguent et se propagent largement à travers celui-ci, si possible. Bien qu'un placement parfait du cadre de sortie ne soit pas requis, je me réserve le droit de refuser d'accepter une réponse qui choisit le cadre d'une manière inacceptable, par exemple. toujours aux mêmes coordonnées, toutes les racines sont en un seul point, etc.
• L'image de sortie doit avoir une taille de 1024 * 1024 pixels.
• Le temps de rendu est de 10 minutes maximum
• L'utilisation de valeurs à virgule flottante simple précision suffit

Sortie

La sortie doit être une image graphique raster dans un format de fichier de votre choix, lisible par un logiciel standard pour un système d'exploitation de marque X. Si vous souhaitez utiliser un format rare, pensez à ajouter un lien vers un site Web où l'on peut télécharger une visionneuse pour celui-ci.

Sortez le fichier sur stdout. Si votre langue ne prend pas en charge la mise de quelque chose sur stdout ou si vous trouvez cette option moins pratique, trouvez un autre moyen. De toute façon, il doit être possible d'enregistrer l'image générée.

Restrictions

• Aucune bibliothèque de traitement d'image
• Aucune bibliothèque génératrice de fractales
• Le code le plus court gagne

Extensions

Si vous aimez cette tâche, vous pouvez essayer de colorer les points en fonction de la vitesse de convergence ou d'autres critères. J'aimerais voir des résultats intéressants.

Python, 827 777 caractères

import re,random
N=1024
M=N*N
R=range
P=map(lambda x:eval(re.sub('i','+',x)+'j'if 'i'in x else x),raw_input().split())[::-1]
Q=[i*P[i]for i in R(len(P))][1:]
E=lambda p,x:sum(x**k*p[k]for k in R(len(p)))
def Z(x):
 for j in R(99):
  f=E(P,x);g=E(Q,x)
  if abs(f)<1e-9:return x,1
  if abs(x)>1e5or g==0:break
  x-=f/g
 return x,0
T=[]
a=9e9
b=-a
for i in R(999):
 x,f=Z((random.randrange(-9999,9999)+1j*random.randrange(-9999,9999))/99)
 if f:a=min(a,x.real,x.imag);b=max(b,x.real,x.imag);T+=[x]
s=b-a
a,b=a-s/2,b+s/2
s=b-a
C=[[255]*3]*M
H=lambda x,k:int(x.real*k)+87*int(x.imag*k)&255
for i in R(M):
 x,f=Z(a+i%N*s/N+(a+i/N*s/N)*1j)
 if f:C[i]=H(x,99),H(x,57),H(x,76)
for r in T:C[N*int(N*(r.imag-a)/s)+int(N*(r.real-a)/s)]=0,0,0
print'P3',N,N,255
for c in C:print'%d %d %d'%c

Les trouvailles affichent les limites (et les racines) en trouvant des points de convergence pour un groupe d'échantillons aléatoires. Il trace ensuite le graphique en calculant les points de convergence pour chaque point de départ et en utilisant une fonction de hachage pour obtenir des couleurs aléatoires pour chaque point de convergence. Regardez de très près et vous pouvez voir les racines marquées.

Voici le résultat de l'exemple de polynôme.

Java, 1093 1058 1099 1077 caractères

public class F{double r,i,a,b;F(double R,double I){r=R;i=I;}F a(F c){return
new F(r+c.r,i+c.i);}F m(F c){return new F(r*c.r-i*c.i,r*c.i+i*c.r);}F
r(){a=r*r+i*i;return new F(-r/a,i/a);}double l(F c){a=r-c.r;b=i-c.i;return
Math.sqrt(a*a+b*b);}public static void main(String[]a){int
n=a.length,i=0,j,x,K=1024,r[]=new int[n];String o="P3\n"+K+" "+K+"\n255 ",s[];F z=new
F(0,0),P[]=new F[n],R[]=new F[n],c,d,e,p,q;for(;i<n;)P[i]=new
F((s=a[i++].split("i"))[0].isEmpty()?0:Float.parseFloat(s[0]),s.length==1?0:Float.parseFloat(s[1]));double
B=Math.pow(P[n-1].m(P[0].r()).l(z)/2,1./n),b,S;for(i=1;i<n;){b=Math.pow(P[i].m(P[i-1].r()).l(z),1./i++);B=b>B?b:B;}S=6*B/K;for(x=0;x<K*K;){e=d=c=new
F(x%K*S-3*B,x++/K*S-3*B);for(j=51;j-->1;){p=P[0];q=p.m(new
F(n-1,0));for(i=1;i<n;){if(i<n-1)q=q.m(c).a(P[i].m(new
F(n-1-i,0)));p=p.m(c).a(P[i++]);}c=c.a(d=q.r().m(p));if(d.l(z)<S/2)break;}i=j>0?0:n;for(;i<n;i++){if(R[i]==null)R[i]=c;if(R[i].l(c)<S)break;}i=java.awt.Color.HSBtoRGB(i*1f/n,j<1||e.l(c)<S&&r[i]++<1?0:1,j*.02f);for(j=0;j++<3;){o+=(i&255)+" ";i>>=8;}System.out.println(o);o="";}}}

L'entrée est des arguments de ligne de commande - par exemple, run java F 1 0 0 -1. La sortie est à stdout au format PPM (pixmap ASCII).

L'échelle est choisie en utilisant la borne de Fujiwara sur la valeur absolue des racines complexes d'un polynôme; Je multiplie ensuite celui lié par 1,5. J'ajuste la luminosité par le taux de convergence, donc les racines seront dans les patchs les plus brillants. Par conséquent, il est logique d'utiliser du blanc plutôt que du noir pour marquer les emplacements approximatifs des racines (ce qui me coûte 41 caractères pour quelque chose qui ne peut même pas être fait "correctement". Si j'étiquette tous les points qui convergent à moins de 0,5 pixel d'eux-mêmes) alors certaines racines sortent sans étiquette; si j'étiquette tous les points qui convergent à moins de 0,6 pixels d'eux-mêmes, alors certaines racines sortent étiquetées sur plus d'un pixel; donc pour chaque racine j'étiquette le premier point rencontré pour converger à moins de 1 pixel de lui-même ).

Image pour l'exemple de polynôme (converti en png avec le GIMP):

Python, 633 octets

import numpy as np
import matplotlib.pyplot as plt
from colorsys import hls_to_rgb
def f(z):
    return (z**4 - 1)
def df(z):
    return (4*z**3) 
def cz(z):
    r = np.abs(z)
    arg = np.angle(z)   
    h = (arg + np.pi)  / (3 * np.pi)
    l = 1.0 - 1.0/(1.0 + r**0.1)
    s = 0.8 
    c = np.vectorize(hls_to_rgb) (h,l,s)
    c = np.array(c)
    c = c.swapaxes(0,2) 
    return c    
x, y = np.ogrid[-1.5:1.5:2001j, -1.5:1.5:2001j]
z = x + 1j*y    
for i in range(10):
    z -= (f(z) / df(z))
zz = z
zz[np.isnan(zz)]=0
zz=cz(zz)
plt.figure()
plt.imshow(zz, interpolation='nearest')
plt.axis('off')
plt.savefig('plots/nf.svg')
plt.close()

Après les accélérations et l'embellissement (756 octets)

import numpy as np
from numba import jit
import matplotlib.pyplot as plt
from colorsys import hls_to_rgb 

@jit(nopython=True, parallel=True, nogil=True)
def f(z):
    return (z**4 - 1)   

@jit(nopython=True, parallel=True, nogil=True)
def df(z):
    return (4*z**3) 

def cz(z):
    r = np.abs(z)
    arg = np.angle(z)   

    h = (arg + np.pi)  / (3 * np.pi)
    l = 1.0 - 1.0/(1.0 + r**0.1)
    s = 0.8 

    c = np.vectorize(hls_to_rgb) (h,l,s)
    c = np.array(c)
    c = c.swapaxes(0,2) 
    return c    

x, y = np.ogrid[-1.5:1.5:2001j, -1.5:1.5:2001j]
z = x + 1j*y    

for i in range(10):
    z -= (f(z) / df(z))

zz = z
zz[np.isnan(zz)]=0
zz=cz(zz)
plt.figure()
plt.imshow(zz, interpolation='nearest')
plt.axis('off')
plt.savefig('plots/nf.svg')
plt.close()

Le graphique ci-dessous concerne la fractale de Newton de la fonction log (z).