Generuj fraktale Newtona

24

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:

Fraktal Newtona dla f (x) = x ^ 3-1 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.

FUZxxl
źródło
6
Nie jestem do końca pewien, czy nadaje się to do golfa kodowego. Moim zdaniem zadanie jest zbyt złożone. Może się jednak okazać, że się mylę.
Joey,
5
@Jey: Rzeczywiście. Chciałbym, żeby to było wyzwanie dla kodu .
Joey Adams,
2
... lub PPM w tym zakresie.
Joey,
1
@Joey: Moim zamiarem było stworzenie raczej trudnego zadania, ponieważ wiele osób nie lubi bardzo łatwych zadań.
FUZxxl
1
Można go łatwo podzielić na osobne zadania, a jeśli Twój język obsługuje natywnie złożone liczby zmiennoprzecinkowe, możesz zapisać dużą porcję. Mam nie w pełni golfową wersję z 1600 znakami, z czego 340 to klasa liczb złożonych. Nie rozpoznaje jeszcze korzeni i nie używa kolorów, ale staram się wyśledzić, co według mnie jest błędem w kodzie NR. (Znalezienie pierwiastka z x ^ 3-1 od -0,5 + 0,866i na pewno nie powinno się rozbierać!)
Peter Taylor

Odpowiedzi:

13

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.

wynik na przykład wielomian

Keith Randall
źródło
Dobry! Lubię to.
FUZxxl
14

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): Korzenie x ^ 5 + (- 2 + 7,5i) x ^ 4 + (23.0004-3.8i) x ^ 3 + 12i x ^ 2 + (5,1233 + 0,1i)

Peter Taylor
źródło
@FUZxxl, obraz pochodzi ze starej wersji. Prześlę jeden ze stopniem konwergencji później. Problem z zaznaczaniem pierwiastków polega na określeniu, który piksel należy oznaczyć. Jest to klasyczny problem polegający na tym, że przy zmiennoprzecinkowym nie można używać dokładnych testów równości, dlatego należy porównać je z epsilonem. W rezultacie nie mam wartości „kanonicznych” dla korzeni. Mógłbym oznaczyć piksele, które zbiegają się w jednym kroku, ale nie gwarantuje to niczego, i mogę zaznaczyć blok 4 pikseli dla jednego katalogu głównego.
Peter Taylor
@Peter Taylor: Jak widzisz, Keith Randall również znalazł rozwiązanie tego problemu. Dodałem ten wymóg jako dodatkową trudność. Jednym z podejść do tego byłoby obliczenie najbliższego piksela dla każdego pierwiastka, a następnie sprawdzenie, czy każdy piksel jest równy.
FUZxxl
@FUZxxl, nie zrozumiałeś mojego punktu. „Najbliższy piksel” pierwiastka nie jest dobrze zdefiniowany. Mogę jednak zhakować coś, co może działać we wszystkich testowanych przez ciebie przypadkach testowych i mam wrażenie, że to cię uszczęśliwi. Zamierzam go pokolorować na biało, a nie na czarno, bo to bardziej logiczne.
Peter Taylor
@Peter Taylor: OK.
FUZxxl
6
Moje zdjęcie profilowe powinno wkrótce się zmienić na x^6-9x^3+8, starannie zaprojektowane przez wybranie korzeni, a następnie użycie Wolfram Alpha w celu uproszczenia wielomianu. Ok, oszukiwałem zamieniając odcienie w GIMP.
Peter Taylor
3

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

Newton Fractal na log (z)

Serowy chleb
źródło
Możesz użyć krótszych nazw (1 znak) i usunąć białe znaki, łącząc wiele linii za pomocą ;. Usuń także wszystkie możliwe spacje.
mbomb007
Niektóre zwykłe gry w golfa zmniejszają to do zaledwie 353 bajtów ! Nie testowałem tego (nie matplotlibtutaj), więc nie ma gwarancji, że nadal działa.
Khuldraeseth na'Barya