Wszyscy znacie metodę Newtona do przybliżania pierwiastków funkcji, prawda? Moim celem w tym zadaniu jest wprowadzenie Cię w interesujący aspekt tego algorytmu.

Algorytm Newtona jest zbieżny tylko dla niektórych, ale przede wszystkim złożonych wartości wejściowych. Jeśli zobrazujesz zbieżność metody dla wszystkich wartości wejściowych na płaszczyźnie złożonej, zwykle otrzymasz piękny fraktal taki jak ten:

Obraz z Wikimedia commons

Dane techniczne

Celem tego zadania jest wygenerowanie takich fraktali. Oznacza to, że otrzymujesz wielomian jako dane wejściowe i musisz wydrukować odpowiedni fraktal jako obraz w wybranym formacie jako wynik.

Wkład

Dane wejściowe to oddzielona spacjami lista liczb zespolonych. Są zapisane w stylu <Real part><iImaginary part>, jak ten numer:5.32i3.05 . Możesz założyć, że liczba wejściowa ma nie więcej niż 4 miejsca po przecinku i jest mniejsza niż 1000. Pierwszy z nich nie może być zerowy. Może to być na przykład dane wejściowe do Twojego programu:

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

Liczby są interpretowane jako współczynniki wielomianu, zaczynając od najwyższej mocy. W pozostałej części tej specyfikacji wejściowy wielomian nosi nazwę P . Powyższe dane wejściowe są równe temu wielomianowi:

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

Dane wejściowe mogą pochodzić z wejścia standardowego, z argumentu przekazanego do programu lub z monitu wyświetlanego programowi. Możesz założyć, że dane wejściowe nie zawierają żadnych początkowych ani końcowych znaków spacji.

Wykonanie

Fraktal musisz wyrenderować w następujący sposób:

• Wybierz tyle kolorów, ile korzeni P. i dodatkowy kolor dla rozbieżności
• Dla każdej liczby w widocznej płaszczyźnie określ, czy metoda jest zbieżna, a jeśli tak, z którym pierwiastkiem. Pokoloruj punkt zgodnie z wynikiem.
• Nie drukuj linijek ani innych fantazyjnych rzeczy
• Wydrukuj czarny punkt w punktach, które są pierwiastkami wielomianów dla orientacji. Możesz wydrukować do czterech pikseli wokół każdego katalogu głównego.
• Znajdź sposób, aby wybrać widoczną płaszczyznę w taki sposób, aby wszystkie korzenie były rozróżnialne i, jeśli to możliwe, szeroko rozłożone na niej. Chociaż nie jest wymagane idealne umieszczenie ramki wyjściowej, zastrzegam sobie prawo do odmowy przyjęcia odpowiedzi, która wybiera ramkę w niedopuszczalny sposób, np. zawsze na tych samych współrzędnych, wszystkie pierwiastki są w jednym punkcie itp.
• Obraz wyjściowy powinien mieć rozmiar 1024 * 1024 pikseli.
• Czas renderowania wynosi maksymalnie 10 minut
• Wystarczy zastosować zmiennoprzecinkowe wartości pojedynczej precyzji

Wydajność

Wyjście powinno być obrazem rastrowym w wybranym formacie pliku, odczytywalnym przez standardowe oprogramowanie dla systemu operacyjnego marki X. Jeśli chcesz użyć rzadkiego formatu, rozważ dodanie linku do strony internetowej, z której można pobrać przeglądarkę.

Wyjście pliku na standardowe wyjście. Jeśli twój język nie obsługuje stdout lub jeśli ta opcja jest mniej wygodna, znajdź inny sposób. W jakikolwiek sposób musi być możliwe zapisanie wygenerowanego obrazu.

Ograniczenia

• Brak bibliotek przetwarzania obrazu
• Brak bibliotek generujących fraktale
• Najkrótszy kod wygrywa

Rozszerzenia

Jeśli podoba Ci się to zadanie, możesz spróbować pokolorować punkty zgodnie z prędkością zbieżności lub innymi kryteriami. Chciałbym zobaczyć ciekawe wyniki.

Python, 827 777 znaków

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

Znajduje granice wyświetlania (i pierwiastki), znajdując punkty zbieżności dla grupy losowych próbek. Następnie rysuje wykres, obliczając punkty zbieżności dla każdego punktu początkowego i używając funkcji skrótu, aby uzyskać losowe kolory dla każdego punktu zbieżności. Przyjrzyj się bardzo uważnie i zobaczysz zaznaczone korzenie.

Oto wynik przykładowego wielomianu.

Java, 1093 1058 1099 1077 znaków

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="";}}}

Dane wejściowe to argumenty wiersza polecenia - np java F 1 0 0 -1. Uruchom . Dane wyjściowe są wysyłane na standardowe wyjście w formacie PPM (piksel ASCII).

Skalę wybiera się za pomocą Fujiwary związanej z wartością bezwzględną złożonych pierwiastków wielomianu; Następnie mnożę to ograniczenie przez 1,5. Dostosowuję jasność za pomocą współczynnika konwergencji, więc korzenie będą w najjaśniejszych łatach. Dlatego logiczne jest używanie bieli zamiast czerni do zaznaczania przybliżonych lokalizacji korzeni (co kosztuje mnie 41 znaków za coś, czego nie można nawet zrobić „poprawnie”. Jeśli oznaczę wszystkie punkty, które zbiegają się w granicach 0,5 piksela od siebie następnie niektóre pierwiastki wychodzą nieoznakowane; jeśli oznaczę wszystkie punkty, które zbiegają się w obrębie 0,6 piksela od siebie, wówczas niektóre pierwiastki wyjdą oznakowane na więcej niż jednym pikselu; więc dla każdego pierwiastka oznaczam pierwszy napotkany punkt, który zbiegnie się w granicach 1 piksela od siebie ).

Obraz dla przykładowego wielomianu (przekonwertowanego do png za pomocą GIMP):

Python, 633 bajty

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()

Po przyspieszeniach i upiększaniu (756 bajtów)

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()

Poniższy wykres dotyczy funkcji Newtona Fraktala log (z).