În [1] am decorat pânza după resturile faţă de un anumit modul ale valorilor unei expresii de coordonate; acum încercăm o altă idee de colorare a punctelor: aplicăm fiecăruia un anumit procedeu de orbitare - care (depinzând de funcţia implicată) se finalizează într-o anumită mulţime restrânsă de puncte "atractoare" - şi colorăm originea fiecărei orbite după "limita" (şi lungimea) acesteia.

Ca "procedeu de orbitare" angajăm metoda lui Newton de aproximare a rădăcinilor unui polinom - obţinând (folosind R) "decoraţiuni" similare celor redate la Newton_fractal; considerăm numai polinoamele zn - 1, care conduc la rădăcinile de ordinul n ale unităţii (şi subliniem că privim lucrurile "cuminte", ca decoraţiuni; teoria din spatele fractalilor este mare şi ne depăşeşte).

Deocamdată, avem o modelare oarecum "defectuoasă" (dar instructivă); prin parametrul 'tip' specificăm modalitatea de colorare (după rădăcina atractoare, sau/şi după lungimea orbitei) şi prevedem (într-o instrucţiune "switch") câte o variantă de funcţie 'orbiting()' pentru fiecare caz (funcţie care primeşte ca argument un punct - şi nu un vector de puncte, cum s-ar cuveni în R - şi iterează după metoda lui Newton, returnând culoarea de asociat acelui punct):

# require(compiler)
# enableJIT(3)  # compilează (inclusiv ciclările) înainte de a executa programul 
eps <- 1E-10
attract <- function(n, tip = 1, x = seq(-1.3, 1.3, length=1000)) {
  roots <- exp(1i * 2*pi * (0:(n-1)/n))  # rădăcinile de ordinul n ale unităţii
  it <- 1  # contorizează iteraţiile orbitării spre o rădăcină, dintr-un punct dat
  switch(tip,
    { # 1 - colorează după rădăcina atractoare
      kol <- randomcoloR::distinctColorPalette(n)  # require("randomcoloR")
      orbiting <- function(z0) {
        while(TRUE) {
          for(r in roots) {  # testează dacă s-a ajuns la una dintre rădăcini
            if(abs(z0 - r) < eps) return(kol[match(r, roots)])
          }
          it <- it + 1
          if(it > 100) break
          z0 <- ((n - 1)*z0 + z0^(1-n)) / n  # tinde spre o rădăcină (metoda lui Newton)
        }
        return("NA")  # dacă nu atinge o rădăcină, nici după 100 de iteraţii
      }
    },
    { # 2 - colorează după lungimea orbitei
      kol <- randomcoloR::distinctColorPalette(100)
      orbiting <- function(z0) {
        while(TRUE) {
          for(r in roots) {
            if(abs(z0 - r) < eps) return(kol[it])
          }
          it <- it + 1
          if(it > 100) break
          z0 <- ((n - 1)*z0 + z0^(1-n)) / n
        }
        return("NA")
      }
    },
    { # 3 - colorează după atractor, gradând spre "white" după lungimea orbitei (modulo 10)
      kol <- randomcoloR::distinctColorPalette(n)
      ramp <- sapply(kol, function(k) colorRampPalette(c(k, "white", k))(10))
      orbiting <- function(z0) {
        while(TRUE) {
          for(r in roots) {
            if(abs(z0 - r) < eps) {
              id <- match(r, roots)  # identifică rădăcina atractoare
              return(ramp[1 + it %% 10, id])
            }
          }
          it <- it + 1
          if(it > 100) break
          z0 <- ((n - 1)*z0 + z0^(1-n)) / n
        }
        return("NA")
      }
    }
  )
  z <- c(outer(x, x, function(u, v) u + 1i*v))
  q <- sapply(z, orbiting)  # culorile punctelor, după valorile parametrului 'tip'
  plot(z, col=q, cex=0.2, pch=19)
}

# exemplu:
opar <- par(mar=rep(0.5, 4), bty="n", xaxt="n", yaxt="n")
prct <- system.time({attract(4, 3)})  # pentru rădăcinile de ordinul 4, {1, i, -1, -i}
print(prct)  # durata execuţiei
par(opar)  # enableJIT(-1)

Imaginile redate corespund polinomului z⁴-1; în prima, punctele colorate la fel orbitează spre o aceeaşi rădăcină (de exemplu, culoarea violet indică rădăcina 1); în a doua, punctele sunt colorate după lungimea orbitei (dar aceasta poate fi determinată numai investigând vectorul de culori folosit; în principiu, era mai potrivit să fi implicat pentru acest caz, un gradient de culoare precum cel returnat de gray.colors(100)); în a treia imagine avem patru culori "principale" (care indică rădăcina atractoare), nuanţate însă după lungimea orbitelor asociate punctelor (şi iarăşi, cu defectul că gradarea %10 a culorii nu reflectă ierarhia firească a acestor lungimi).

Am evidenţiat mai sus vreo patru defecte: am scris de trei ori (doar cu mici modificări) o aceeaşi funcţie "orbiting()"; această funcţie operează pe un singur punct (şi tocmai de aceea am adaptat-o după "tip", în cele trei variante) - ca urmare, a trebuit să o implicăm indirect, sub sapply() (nefiind "vectorizată", nu poate fi invocată direct în apelul outer()). Apoi (mai puţin important), gradarea culorilor este defectuoasă.

Dar s-a mai ivit un defect, deja stringent: durata execuţiei programului este mare. Tocmai de aceea am încercat compiler::enableJIT(), prin care fişierul-sursă a programului este în prealabil, compilat (cu anumite optimizări) - şi într-adevăr, timpul de execuţie se reduce astfel, măcar cu 15%; totuşi, pentru fiecare din imaginile redate mai sus, execuţia a durat cam 40 de secunde.

În plus, să zicem - investigând culorile atribuite punctelor (după reducerea valorii "length" din parametrul "x", la 100, la 200, etc., variind totodată numărul maxim prevăzut pentru iteraţii) am constatat că există puncte colorate "NA" (şi anume, cele aflate pe diagonalele x = ±y), însemnând că orbitele corespunzătoare nu ajung în nici una dintre rădăcini; s-ar fi cuvenit ca aceste puncte să fie reliefate cumva, pe imaginea produsă (mai ales dacă ne gândim şi la posibilele puncte "NA" pentru alte polinoame, decât cele corespunzătoare rădăcinilor unităţii).