use strict;
use warnings;

use PDL;
use PDL::Graphics::Simple;

use Inline Pdlpp => <<'EOPP';
pp_def('pp_mandel',
    # Pars (signature) specs threadable arguments.
    Pars => 'coord(n=2); [o]o(); complex [t]c()',
    # OtherPars specs scalar arguments.
    OtherPars => 'int max_it',
    Code => <<'EOC',
    /* All this code gets wrapped automagically in a thread loop.  */
    /* The $GENERIC(c) macro is the type of ndarray "c".        */
    // using the Wikipedia formula: z_{n+1} = {z_n}^2 + c
    int i;                        // iterator.
    $c() = $coord(n=>0) + I*$coord(n=>1); // Copy the initial value
    $GENERIC(c) z = $c(); // Copy for the initial iteration
    // the OtherPars are in the $COMP macro.
    for(i=$COMP(max_it); cabs(z) < 2 && i; i--)
        z = z*z + $c();
    $o()= i; // Assign the iterator to the output value
EOC
);
EOPP

my $cen = pdl(-0.74897,0.05708);
my $coords = $cen + (ndcoords(501,501)/250 - 1) * 0.001;
my $w = pgswin();

$w->image( $coords->using(0,1), $coords->pp_mandel(2500),
  {title=>"Mandelbrot"} );

<STDIN>