So I finally figured it out

I’d been trying for a while in my spare time to understand why my incredibly simple Perl Mandelbrot test, inspired by the Julia benchmarks, was returning wrong numbers. Yeah, they were wrong. As in incorrect values.

So I figured it out this morning. The punchline. There is a bug (which I haven’t quite yet found) in the Math::Complex library, specifically in the pathway for the abs(z) function.

How did I find this.

Start with the C code

int mandel(double complex z) {
    int maxiter = 80;
    double complex c = z;
    for (int n=0; n<maxiter ; ++n) {
        if (cabs(z) > 2.0) {
            return n;
        }
        z = z*z+c;
    }
    return maxiter;
}
</maxiter>

Very simple loop. Convert this to Perl

sub mandel {
    use Math::Complex;
    my $z = shift;
    my $c = $z;
    my $n;
    my $maxiter = 80;
    for($n=0; $n < $maxiter; $n++) {
	 if ( abs($z) > 2.0 ) { return $n  }
     $z = $z * $z + $c ;
    }
    return $maxiter;
}

Run both codes with their drivers.

C:

landman@lightning:~/work/benchmarking/Julia$ ./mandel.exe 
c,mandel,0.399828
sum: 14791


Perl:

landman@lightning:~/work/benchmarking/Julia$ ./m2.pl 
perl,mandel,187.34
 sum mandel = 14722

The sum is both a sanity check and an optimizer defeater. It should be the same. And it is, in every other language.

So what is going on here.

I finally had time to think this through, and follow good debugging practice. Remove the impossible from consideration, until the only thing that remains is the answer.

Instrumenting everything that made sense, comparing results along the iteration pathway led me to believe that there is an issue with the Math::Complex code. So I started out replacing key aspects of it.

Finally, I changed the abs($z) portion out for its equivalent, and simplified the expression

sub mandel {
    use Math::Complex;
    my $z = shift;
    my $c = $z;
    my $maxiter = 80;
    for(my $n=0; $n < $maxiter; $n++) {
        if ( (Re($z)**2+Im($z)**2) > 4.0 ) { return $n  }
        $z *= $z;
        $z += $c;
    }
    return $maxiter;
}

Basically I removed the call to abs($z)

landman@lightning:~/work/benchmarking/Julia$ ./m2.pl 
perl,mandel,364.876
 sum mandel = 14791

Now swapping the abs($z) code back into that one line …

sub mandel {
    use Math::Complex;
    my $z = shift;
    my $c = $z;
    my $maxiter = 80;
    for(my $n=0; $n < $maxiter; $n++) {
        #if ( (Re($z)**2+Im($z)**2) > 4.0 ) { return $n  }
        if ( abs($z) > 2.0 ) { return $n  }
        $z *= $z;
        $z += $c;
    }
    return $maxiter;
}

we get

landman@lightning:~/work/benchmarking/Julia$ ./m2.pl 
perl,mandel,192.367
 sum mandel = 14722

Note how the sums are wrong.

Ok. Looking at the abs function in Math::Complex, it basically overloads the original core library abs code, and decides between abs(real), abs(cartesian) and abs(polar)

sub abs {
        my ($z, $rho) = @_ ? @_ : $_;
        unless (ref $z) {
            if (@_ == 2) {
                $_[0] = $_[1];
            } else {
                return CORE::abs($z);
            }
        }
        if (defined $rho) {
            $z->{'polar'} = [ $rho, ${$z->_polar}[1] ];
            $z->{p_dirty} = 0;
            $z->{c_dirty} = 1;
            return $rho;
        } else {
            return ${$z->_polar}[0];
        }
}

Following through the code, it hits z->_polar, which is

sub _polar     {$_[0]->{p_dirty} ?
		   $_[0]->_update_polar : $_[0]->{'polar'}}

and we’ll hit _update_polar

sub _update_polar {
	my $self = shift;
	my ($x, $y) = @{$self->{'cartesian'}};
	$self->{p_dirty} = 0;
	return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
	return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y),
				   CORE::atan2($y, $x)];
}

where the actual computation is done.

Ok. Deep call stack for something that should be fast and called frequently.

This is bad design IMO. I understand the desire to handle r,θ type coordinates as well as cartesian. But this introduces many issues in the library. Unacceptable tradeoffs IMO.

Core library code should be simple, fast, and as much as possible bug free. Even more than this, there is a whole extra computation that is not needed here (note the atan2 call).

All I can say is … wow.

This one is bad enough that I might rewrite this library, stripping it down to the bare essentials for cartesian complex numbers.

But more frightening in this is that I am getting what appears to be a failure in a Core:: function. Which means I probably have to look through this as well.

Enough for now.

Viewed 41774 times by 5199 viewers

Facebooktwittergoogle_plusredditpinterestlinkedinmail