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 divir 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:

video

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.

segunda-feira, 30 de janeiro de 2017

Conectando Postgresql a Oracle com Perl

Postgresql oferece a possibilidade de mapear tabelas e vistas via FDWs (Foreign Data Wrappers). Isso, por si só, já é muito útil. Entretanto, melhor ainda seria poder invocar rotinas. Como o Postgresql permite escrever rotinas em Python, TCL, e Perl, a saída é óbvia.


create or replace 
function call_remote_function(varchar, varchar[]) returns varchar as $$
  use DBI;
  use DBD::Oracle;
  
  my $proc=shift;
  $proc=~s/[^A-Za-z0-9_.]//g;
  my $vars=shift;
  my $num_vars=scalar(@{$vars});
  my $placeholders=join(',',('?')x$num_vars);
    
  my $host=$ENV{'REMOTE_HOST'};
  my $schema=$ENV{'REMOTE_SCHEMA'};
  my $sid=$ENV{'REMOTE_SID'};
  my $pass=$ENV{'REMOTE_PASS'};
  
  my $db = DBI->connect("dbi:Oracle:host=$host;sid=$sid;server=POOLED", 
                        "$schema", "$pass", {ora_drcp=>1})
    || die($DBI::errstr."\n");
  $db->{AutoCommit}=0;
  $db->{RaiseError}=1;

  my $result;
  my $cmd="BEGIN ?:=$proc($placeholders); END;";
  my $stmt=$db->prepare($cmd);
  my $i=1;
  $stmt->bind_param_inout($i, \$result, {ora_type=>SQLT_CHR});
  for my $var (@{$vars}) {
    $stmt->bind_param(++$i, $var);
  };
  $stmt->execute();
  $db->disconnect();
  return $result;
$$ language plperlu;

Esta é uma versão que vale para qualquer número de parâmetros. Escrevi também uma versão para nenhum parâmetro, um parâmetro, e dois parâmetros.

Como essa função faz algumas coisas perigosas, ela só pode ser criada por uma conta administrativa. Além disso, a linguagem escolhida é a plperlu, quando rotinas normais podem ser declaradas com plperl.

Os parâmetros são inseridos como variáveis de ligação (bind variables), então não me preocupei muito com a segurança delas. O nome da rotina, entretanto é inserida diretamente no comando. Então, adicionei uma linha para limpar caracteres indesejados:

  $proc=~s/[^A-Za-z0-9_.]//g;

Essa linha retira todos os caracteres que não forem alfabéticos, numéricos, traço-baixo, ou ponto.

Os dados de conexão são escondidos como variáveis de ambiente do usuário que roda o banco. Além disso, no lado do Oracle, defini um pool de conexões (usando DRCP - Database Resident Connection Pooling). O tempo de execução das chamadas variava entre 300ms e 500ms. Com o DRCP, o teto passou a ser 300ms e normalmente elas são invocadas entre 100ms e 200ms.

No banco Oracle, adicionei uma entrada ao tnsnames.ora:

DB_POOLED =
  (DESCRIPTION =
    (ADDRESS = (PROTOCOL = TCP)(HOST = host.com)(PORT = 1521))
    (CONNECT_DATA =
      (SERVER = POOLED)
      (SID = DB)
    )
  )

E depois iniciei a pool com um comando:

execute dbms_connection_pool.start_pool();


Posso invocar uma função remota assim:

select call_remote_function('substr', array['postgresql', '2', '3']);

Mas prefiro criar funções locais para esconder a pilantragem:

CREATE OR REPLACE FUNCTION 
oracle_substr(p_str varchar, p_start varchar, p_length varchar) RETURNS varchar AS $$
begin
  return call_remote_function('substr', array[p_str, p_start, p_length]);
end;
$$ LANGUAGE plpgsql;

quinta-feira, 5 de janeiro de 2017

2017

Este ano é especial porque, além de ser primo, eu farei um aniversário primo. Então, resolvi descobrir quantos aniversários primos eu já completei.

Escrevi, para esta tarefa um pouco de Perl:

#!/usr/bin/perl
my @primes=qw(
   2 3 5 7 11 13 17 19 23 29
31 37 41 43 47 53 59 61 67 71
73 79 83 89 97 101 103 107 109 113
127 131 137 139 149 151 157 163 167 173
179 181 191 193 197 199 211 223 227 229
233 239 241 251 257 263 269 271 277 281
283 293 307 311 313 317 331 337 347 349
353 359 367 373 379 383 389 397 401 409
419 421 431 433 439 443 449 457 461 463
467 479 487 491 499 503 509 521 523 541
547 557 563 569 571 577 587 593 599 601
607 613 617 619 631 641 643 647 653 659
661 673 677 683 691 701 709 719 727 733
739 743 751 757 761 769 773 787 797 809
811 821 823 827 829 839 853 857 859 863
877 881 883 887 907 911 919 929 937 941
947 953 967 971 977 983 991 997 1009 1013
1019 1021 1031 1033 1039 1049 1051 1061 1063 1069
1087 1091 1093 1097 1103 1109 1117 1123 1129 1151
1153 1163 1171 1181 1187 1193 1201 1213 1217 1223
1229 1231 1237 1249 1259 1277 1279 1283 1289 1291
1297 1301 1303 1307 1319 1321 1327 1361 1367 1373
1381 1399 1409 1423 1427 1429 1433 1439 1447 1451
1453 1459 1471 1481 1483 1487 1489 1493 1499 1511
1523 1531 1543 1549 1553 1559 1567 1571 1579 1583
1597 1601 1607 1609 1613 1619 1621 1627 1637 1657
1663 1667 1669 1693 1697 1699 1709 1721 1723 1733
1741 1747 1753 1759 1777 1783 1787 1789 1801 1811
1823 1831 1847 1861 1867 1871 1873 1877 1879 1889
1901 1907 1913 1931 1933 1949 1951 1973 1979 1987
1993 1997 1999 2003 2011 2017 2027 2029 2039 2053
2063 2069 2081 2083 2087 2089 2099 2111 2113 2129
2131 2137 2141 2143 2153 2161 2179 2203 2207 2213);

my %primes=map { $_ => 1 } @primes;

my $birth=1966;

for my $i (1..100) {
  if (exists $primes{$i} && exists $primes{$birth+$i}) {
    print "$i ".($birth+$i)."\n";
  }
}

Eu estou no meu sétimo aniversário primo em ano primo. Para uma pessoa nascida em 1966, o resultado seria este:

7 1973
13 1979
31 1997
37 2003
61 2027
73 2039
97 2063


Curiosamente, muitas pessoas nunca experimentarão esse evento.

Então, como uma pergunta sempre leva a outra, pensei: quantos anos há em que as pessoas que neles nasceram nunca experimentarão um aniversário primo num ano primo e qual o ano que levaria ao maior número desses eventos? Testei do ano 1 ao ano 2000.

Do ano 1 ao ano 2000, as pessoas nascidas em 1302 diferentes anos tiveram essa sorte. Se seu reduzisse para 70 anos a idade máxima, seriam 1301 anos. Então, cerca de um terço das pessoas nunca passará por isso.

Os que nasceram nos anos 6AD e 30AD podiam ter 13 aniversários primos em anos primos, se vivessem até os 70 anos. Se eu aumentar para 100 anos a expectativa de vida, os nascidos no ano 30AD teriam 18 desses eventos.

No segundo milênio, os melhores anos para nascer foram 1410 e 1470 (14 aniversários primos em anos primos para quem viver até os 100 anos). No terceiro milênio, os anos mais interessantes são 2640 e 2670 (15 eventos para quem viver até os 100). Nos últimos 100 anos (de 1917 a 2017), os anos de 1980 e 1956 resultam em 12 eventos.

Os anos que nunca fazem isso são especiais também. Alguns, como 1943, 1975, e 1979 não passam por isso nem em 100 mil anos, o que eu poderia ter deduzido, já que sempre fazem aniversário par em ano ímpar. Nascidos em 1977 terão o prazer de fazer 2 em 1979 e só. Então, quem nasceu num ano ímpar, só pode fazer um aniversário primo num ano primo no seu segundo aniversário. Isso reduz bastante as possibilidades e torna mais interessante o fato de nascer dois anos antes de um ano primo.