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