quinta-feira, 6 de abril de 2017

Mandelbrot com expoentes complexos e animados

Depois de descobrir como gerar um mandelbrot com expoentes complexos, caí no dilema de como fazer uma animação disso. O dilema é qual intervalo escolher. Tendo duas dimensões, os números complexos podem variar tanto no eixo real quanto no eixo imaginário. Há infinitas variações possíveis.

Dado que a exponenciação com números imaginários tem o efeito de girar os números, decidi animar usando vetores de tamanho 2 em torno do zero. Então, eu começaria com 2 (o mandelbrot clássico). Girando o vetor por 90 graus no sentido anti-horário, eu chegaria a 0+2i. Mais 90 graus me levariam a -2, depois 0-2i e, finalmente, com 360 graus, eu estaria de volta a 2.

Então, a parte real varia com o cosseno e a parte imaginária varia com o seno. Basta dividir os 2*pi radianos do círculo pelo número de quadros que quisermos.

Testei primeiro com 4 passos e obtive as seguintes imagens:





Pareceu promisor. Então, decidi fazer uma animação com 1000 quadros.

use GD;
use strict;
use constant PI => 4 * atan2(1, 1);
use constant E => exp(1);
use constant MAXITER => 256;
use constant SIDE => 512;
use constant FRAMES => 1000;

sub mandel($$$$) {
  my $c=1;
  my ($x, $y, $rp, $ip)=@_;
  my $r=$x;
  my $i=$y;
  while($r**2+$i**2<4 && $c<MAXITER) {
    my $a=atan2($i,$r);
    my $b=$r**2+$i**2;
    my $p=($b**($rp/2))*(E**(-1*$ip*$a));
    my $nr=$p*cos($rp*$a+0.5*$ip*($b!=0?log($b):0));
    my $ni=$p*sin($rp*$a+0.5*$ip*($b!=0?log($b):0));
    $r=$nr+$x;
    $i=$ni+$y;
    $c++;
  }
  return $c;
}

my $img=new GD::Image(SIDE,SIDE); 
my @colours=();
for my $c (0..255) {
  $colours[$c]=$img->colorAllocate(($c*7)%255, ($c*13)%255, ($c*19)%255),
};

for my $frame (0..FRAMES-1) {
  my $theta=0;
  if($frame>0) {
    $theta=(2*PI)*($frame/FRAMES);
  }
  my $ip=2*sin($theta);
  my $rp=2*cos($theta);
  for my $a (0..SIDE) {
    for my $b (0..SIDE) {
      my $d=mandel(-2+(4*$a/SIDE), -2+(4*$b/SIDE), $rp, $ip);
      $img->setPixel($a, $b, $colours[$d]);
    }
  } 
  open(IMAGE, sprintf(">frame%04d.png", $frame));
  binmode IMAGE;
  print IMAGE $img->png;
  close IMAGE;
}

O código é simples e a biblioteca GD torna fácil a geração de PNGs (ou GIFs ou JPEGs...). Depois, usei ffmpeg para juntar tudo numa animação.

ffmpeg  -i frame%04d.png -vf curves=preset=lighter video.mp4

E o resultado final é este vídeo:


O vídeo está no youtube também.

segunda-feira, 3 de abril de 2017

Cobertura dos n primeiros primos II

Quando calculei a porcentagem de números inteiros que possuem fatores dentre os n primeiros primos, procurei calcular diretamente o valor. Deixei escapar uma maneira muito mais simples que é a de calcular quantos números não são cobertos pelos n primeiros primos.

Esse cálculo é mais simples. Calcular da maneira mais complexa me ensinou algumas coisas, mas a solução mais simples também vale a pena.

Os números que não são múltiplos de 2 são 1/2 do total. Os que não são múltiplos de 3 são 2/3 do total. Os que não são múltiplos de 5 são 4/5. E assim por diante.

Os que não são múltiplos de 2, 3, ou 5, são (1-1/2)(1-1/3)(1-1/5) do total. É muito mais simples a conta.

Em menos de 1s, um programa calculou que os primeiros 10.000 primos fazem parte de 95,14% de todos os inteiros. O programa original não consegue chegar a n=13 nesse tempo. Os primeiros 100.008 ajudam a quebrar a barreira dos 96%, mas o primeiro milhão não chega a 97% (96,6%).

O gráfico abaixo mostra o resultado dessas contas. Ele mostra, em escala, logarítmica, quantos inteiros não têm um fator dentre os primeiros n primos.



O código é tão insosso quanto o gráfico, então não vou publicá-lo.