#!/usr/bin/perl -w
# Copyright 2010 Kevin Ryde
# This file is part of Image-Base.
#
# Image-Base is free software; you can redistribute it and/or modify it
# under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation; either version 3, or (at your option) any
# later version.
#
# Image-Base is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with Image-Base. If not, see <http://www.gnu.org/licenses/>.
use strict;
use warnings;
use List::Util qw(min max);
# uncomment this to run the ### lines
use Smart::Comments;
my $a = 1000;
my $b = 1000;
my $aa = $a ** 2;
my $bb = $b ** 2;
my $x = $a - int($a) ; # 0 or 0.5
my $y = $b ;
### initial: "start xy $x,$y"
my $d = ($x ? 2.25*$bb : $bb) - ( $aa * $b ) + ( $aa / 4 ) ;
my $max_d = 0;
while( $y >= 1
&& ( $aa * ( $y - 0.5 ) ) > ( $bb * ( $x + 1 ) ) ) {
### assert: $d == ($x+1)**2 * $bb + ($y-.5)**2 * $aa - $aa * $bb
if( $d < 0 ) {
$d += ( $bb * ( ( 2 * $x ) + 3 ) ) ;
++$x ;
}
else {
$d += ( ( $bb * ( ( 2 * $x ) + 3 ) ) +
( $aa * ( ( -2 * $y ) + 2 ) ) ) ;
++$x ;
--$y ;
}
$max_d = max($d, $max_d);
}
# switch to d2 = E(x+1/2,y-1) by adding E(x+1/2,y-1) - E(x+1,y-1/2)
$d += $aa*(.75-$y) - $bb*($x+.75);
while( $y >= 1 ) {
if( $d < 0 ) {
$d += ( $bb * ( ( 2 * $x ) + 2 ) ) +
( $aa * ( ( -2 * $y ) + 3 ) ) ;
++$x ;
--$y ;
}
else {
$d += ( $aa * ( ( -2 * $y ) + 3 ) ) ;
--$y ;
}
### assert: $d == $bb*($x+0.5)**2 + $aa*($y-1)**2 - $aa*$bb
$max_d = max($d, $max_d);
}
### $max_d