Math-PlanePath
view release on metacpan or search on metacpan
lib/Math/PlanePath/DiagonalsOctant.pm view on Meta::CPAN
}
my $direction = ($self->{'direction'} ||= 'down');
if (! ($direction eq 'up' || $direction eq 'down')) {
croak "Unrecognised direction option: ", $direction;
}
return $self;
}
# start from integers
# [ 1, 2, 3, 4, 5 ],
# [ 1, 3, 7, 13, 21 ]
# N = (d^2 - d + 1)
# = ($d**2 - $d + 1)
# = (($d - 1)*$d + 1)
#
# starting 0.5 back from the odd Y, numbering d=1 for the Y=1
# [ 1, 2, 3, 4, 5 ],
# [ 2-.5, 5-.5, 10-.5, 17-.5, 26-.5 ] # squares +0.5
# N = (d^2 + 1/2)
# = ($d**2 + 1/2)
# = ($d**2 + 1/2)
# d = 0 + sqrt(1 * $n + -1/2)
# = sqrt($n -1/2)
# = sqrt(4*$n-2)/2
#
# 9 | 26
# 8 | 21
# 7 | 17 22
# 6 | 13 18 23
# 5 | 10 14 19 24
# 4 | 7 11 15 20 25
# 3 | 5 8 12 16
# 2 | 3 6 9
# 1 | 2 4
# Y=0 | 1
# +-------------
#
sub n_to_xy {
my ($self, $n) = @_;
### DiagonalsOctant n_to_xy(): "$n ".(ref $n || '')
# adjust to N=1 at origin X=0,Y=0
$n = $n - $self->{'n_start'} + 1;
my $d = int(4*$n)-2; # for sqrt
if ($d < 0) {
### nothing at N < 0.5 ...
return;
}
$d = int( _sqrtint($d)/2 );
### $d
# remainder positive or negative relative to the start of the following
# diagonal
#
$n -= $d*($d+1) + 1;
### remainder: $n
# $n first in formulas to preserve n=BigFloat when d=integer is BigInt
#
if ($self->{'direction'} eq 'up') {
if (2*$n >= -1) {
return (-$n + $d,
$n + $d);
} else {
return (-$n - 1,
$n + 2*$d);
}
} else {
if (2*$n >= -1) {
# stripe of d+1 many points starting at Y even, eg. N=13
return ($n,
-$n + 2*$d);
} else {
# stripe of d many points starting at Y odd, eg. N=10
return ($n + $d,
-$n + $d - 1);
}
}
}
sub xy_to_n {
my ($self, $x, $y) = @_;
### xy_to_n(): $x, $y
$x = round_nearest ($x);
if ($self->{'direction'} eq 'up') {
$y = round_nearest ($y);
} else {
$y = - round_nearest (- $y);
}
### rounded
### $x
### $y
if ($y < 0 || $x < 0 || $x > $y) {
### outside upper octant ...
return undef;
}
if ($self->{'direction'} eq 'up') {
my $d = $x + $y + 2;
### $d
return ($d*$d - ($d % 2))/4 - $x + $self->{'n_start'} - 1;
} else {
my $d = $x + $y + 1;
### $d
return ($d*$d - ($d % 2))/4 + $x + $self->{'n_start'};
}
}
# exact
sub rect_to_n_range {
my ($self, $x1,$y1, $x2,$y2) = @_;
$x1 = round_nearest ($x1);
$x2 = round_nearest ($x2);
if ($self->{'direction'} eq 'up') {
lib/Math/PlanePath/DiagonalsOctant.pm view on Meta::CPAN
=head1 FUNCTIONS
See L<Math::PlanePath/FUNCTIONS> for behaviour common to all path classes.
=over 4
=item C<$path = Math::PlanePath::DiagonalsOctant-E<gt>new ()>
=item C<$path = Math::PlanePath::DiagonalsOctant-E<gt>new (direction =E<gt> $str, n_start =E<gt> $n)>
Create and return a new path object.
=item C<($x,$y) = $path-E<gt>n_to_xy ($n)>
Return the X,Y coordinates of point number C<$n> on the path.
For C<$n E<lt> 0.5> the return is an empty list, it being considered the
path begins at 1.
=item C<$n = $path-E<gt>xy_to_n ($x,$y)>
Return the point number for coordinates C<$x,$y>. C<$x> and C<$y> are
each rounded to the nearest integer, which has the effect of treating each
point C<$n> as a square of side 1.
=item C<($n_lo, $n_hi) = $path-E<gt>rect_to_n_range ($x1,$y1, $x2,$y2)>
The returned range is exact, meaning C<$n_lo> and C<$n_hi> are the smallest
and biggest in the rectangle.
=back
=head1 FORMULAS
=head2 N to X,Y
To break N into X,Y it's convenient to take two diagonals at a time, since
the length then changes by 1 each pair making a quadratic. Starting at each
X=0,Y=odd just after perfect square N allows just a sqrt.
Nstart = d*d+1
where d numbers diagonal pairs, eg. d=3 for X=0,Y=5 going down. This is
easily reversed as
d = floor sqrt(N-1)
The code reckons the start of the diagonal as 0.5 further back, so that
N=9.5 is at X=-.5,Y=5.5. To do that d is formed as
d = floor sqrt(N-0.5)
= int( sqrt(int(4*$n)-2)/2 )
Taking /2 out of the sqrt helps with C<Math::BigInt> which circa Perl 5.14
doesn't inter-operate with flonums very well.
In any case N-Nstart is an offset into two diagonals, the first of length d
many points and the second d+1. For example d=3 starting Y=5 for points
N=10,11,12 followed by Y=6 N=13,14,15,16.
The formulas are simplified by calculating a remainder relative to the
second diagonal, so it's negative for the first and positive for the second,
Nrem = N - (d*(d+1)+1)
d*(d+1)+1 is 1 past the pronic numbers when end each first diagonal, as
described above. In any case for example d=3 is relative to N=13 making
Nrem=-3,-2,-1 or Nrem=0,1,2,3.
To include the preceding 0.5 in the second diagonal simply means reckoning
NremE<gt>=-0.5 as belonging to the second. In that base
if Nrem >= -0.5
X = Nrem # direction="down"
Y = 2*d - Nrem
else
X = Nrem + d
Y = d - Nrem - 1
For example N=15 Nrem=1 is the first case, X=1, Y=2*3-1=5. Or N=11 Nrem=-2
the second X=-2+3=1, Y=3-(-2)-1=4.
For "up" direction the Nrem and d are the same, but the coordinate
directions reverse.
if Nrem >= -0.5
X = d - Nrem # direction="up"
Y = d + Nrem
else
X = -Nrem - 1
Y = 2d + Nrem
Another way is to reckon Nstart from the X=0,Y=even diagonals, which is then
two diagonals of the same length and d formed by a sqrt inverting a pronic
Nstart=d*(d+1).
=head2 Rectangle to N Range
Within each row increasing X is increasing N, and in each column increasing
Y is increasing N. This is so in both "down" and "up" arrangements. On
that basis in a rectangle the lower left corner is the minimum N and the
upper right is the maximum N.
If the rectangle is partly outside the covered octant then the corners must
be shifted to put them in range, ie. trim off any rows or columns entirely
outside the rectangle. For the lower left this means,
| | /
| | /
+-------- if x1 < 0 then x1 = 0
x1 | / increase x1 to within octant
|/
+
| |/
| | if y1 < x1 then y1 = x1
| /| increase y1 to bottom-left within octant
|/ +----y1
+ x1
And for the top right,
( run in 0.531 second using v1.01-cache-2.11-cpan-df04353d9ac )