1#!/usr/bin/perl -w
2
3# Copyright 2010 Kevin Ryde
4
5# This file is part of Image-Base.
6#
7# Image-Base is free software; you can redistribute it and/or modify it
8# under the terms of the GNU Lesser General Public License as published by
9# the Free Software Foundation; either version 3, or (at your option) any
10# later version.
11#
12# Image-Base is distributed in the hope that it will be useful,
13# but WITHOUT ANY WARRANTY; without even the implied warranty of
14# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
15# Public License for more details.
16#
17# You should have received a copy of the GNU Lesser General Public License
18# along with Image-Base.  If not, see <http://www.gnu.org/licenses/>.
19
20use strict;
21use warnings;
22use List::Util qw(min max);
23
24# uncomment this to run the ### lines
25use Smart::Comments;
26
27
28my $a = 1000;
29my $b = 1000;
30my $aa = $a ** 2;
31my $bb = $b ** 2;
32
33my $x  = $a - int($a) ;  # 0 or 0.5
34my $y  = $b ;
35### initial: "start xy $x,$y"
36
37my $d = ($x ? 2.25*$bb : $bb) - ( $aa * $b ) + ( $aa / 4 ) ;
38
39my $max_d = 0;
40
41while( $y >= 1
42       && ( $aa * ( $y - 0.5 ) ) > ( $bb * ( $x + 1 ) ) ) {
43
44  ### assert: $d == ($x+1)**2 * $bb + ($y-.5)**2 * $aa - $aa * $bb
45  if( $d < 0 ) {
46    $d += ( $bb * ( ( 2 * $x ) + 3 ) ) ;
47    ++$x ;
48  }
49  else {
50    $d += ( ( $bb * ( (  2 * $x ) + 3 ) ) +
51            ( $aa * ( ( -2 * $y ) + 2 ) ) ) ;
52    ++$x ;
53    --$y ;
54  }
55  $max_d = max($d, $max_d);
56}
57
58# switch to d2 = E(x+1/2,y-1) by adding E(x+1/2,y-1) - E(x+1,y-1/2)
59$d += $aa*(.75-$y) - $bb*($x+.75);
60
61while( $y >= 1 ) {
62  if( $d < 0 ) {
63    $d += ( $bb * ( (  2 * $x ) + 2 ) ) +
64      ( $aa * ( ( -2 * $y ) + 3 ) ) ;
65    ++$x ;
66    --$y ;
67  }
68  else {
69    $d += ( $aa * ( ( -2 * $y ) + 3 ) ) ;
70    --$y ;
71  }
72  ### assert: $d == $bb*($x+0.5)**2 + $aa*($y-1)**2 - $aa*$bb
73
74  $max_d = max($d, $max_d);
75}
76
77### $max_d
78