File Coverage

blib/lib/Math/PlanePath/GcdRationals.pm
Criterion Covered Total %
statement 137 159 86.1
branch 35 48 72.9
condition 12 23 52.1
subroutine 27 30 90.0
pod 7 7 100.0
total 218 267 81.6


line stmt bran cond sub pod time code
1             # Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 Kevin Ryde
2              
3             # This file is part of Math-PlanePath.
4             #
5             # Math-PlanePath is free software; you can redistribute it and/or modify
6             # it under the terms of the GNU General Public License as published by the
7             # Free Software Foundation; either version 3, or (at your option) any later
8             # version.
9             #
10             # Math-PlanePath is distributed in the hope that it will be useful, but
11             # WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12             # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13             # for more details.
14             #
15             # You should have received a copy of the GNU General Public License along
16             # with Math-PlanePath. If not, see .
17              
18              
19             # A003989 diagonals from (1,1)
20             # A109004 0,1,1,2,1,2,3,1,1,3,4,1,2,1,4,5,1,1,1,1
21             # gcd by diagonals (0,0)=0
22             # (1,0)=1 (0,1)=1
23             # (2,0)=2 (1,1)=1 (0,2)=2
24             # A050873 gcd rows n>=1, k=1..n
25             # 1,1,2,1,1,3,1,2,1,4,1,1,1,1,5,1,2,3,2,1,6,1,1,1,
26             # add 0,1,0,1,1,0,1,1,1,0,1,1,1,1,0,1,1,1,1,1,0 A023532 0 at m(m+3)/2
27             # IntXY 1,0,2,0,0,3,0,1,0,4,0,0,0,0,5,0,1,2,1,0,6,
28             # IntXY+1 2,1,3,1,1,4,1,2,1,5,1,1,1,1,6,1,2,3,2,1,7
29             # diff 1,0,1,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,0,0,1 A023531
30             # A178340 1,2,1,3,1,1,4,1,2,1,5,1,1,1,1,6,1,2,3,2,1,7,1,1 Bernoulli
31             # T(n,m) = A003989(n-m+1,m) m>=1, except when factor cancels
32              
33             # diagonals_down even/odd in wedges, and other modulo
34              
35             # math-image --path=GcdRationals --expression='i<30*31/2?i:0' --text --size=40
36             # math-image --path=GcdRationals --output=numbers --expression='i<100?i:0'
37             # math-image --path=GcdRationals --all --output=numbers
38              
39             # Y = v = j/g
40             # X = (g-1)*v + u
41             # = (g-1)*j/g + i/g
42             # = ((g-1)*j + i)/g
43              
44             # j=5 11 ...
45             # j=4 7 8 9 10
46             # j=3 4 5 6
47             # j=2 2 3
48             # j=1 1
49             #
50             # N = (1/2 d^2 - 1/2 d + 1)
51             # = (1/2*$d**2 - 1/2*$d + 1)
52             # = ((1/2*$d - 1/2)*$d + 1)
53             # j = 1/2 + sqrt(2 * $n + -7/4)
54             # = [ 1 + 2*sqrt(2 * $n + -7/4) ] /2
55             # = [ 1 + sqrt(8*$n -7) ] /2
56             #
57              
58             # Primes
59             # i=3*a,j=3*b
60             # N=3*a*(3*b-1)/2
61              
62              
63             package Math::PlanePath::GcdRationals;
64 5     5   6863 use 5.004;
  5         20  
65 5     5   28 use strict;
  5         11  
  5         186  
66 5     5   40 use Carp 'croak';
  5         11  
  5         424  
67             #use List::Util 'min','max';
68             *min = \&Math::PlanePath::_min;
69             *max = \&Math::PlanePath::_max;
70              
71 5     5   33 use vars '$VERSION', '@ISA';
  5         7  
  5         350  
72             $VERSION = 129;
73 5     5   898 use Math::PlanePath;
  5         10  
  5         272  
74             *_sqrtint = \&Math::PlanePath::_sqrtint;
75             @ISA = ('Math::PlanePath');
76              
77             use Math::PlanePath::Base::Generic
78 5         299 'is_infinite',
79 5     5   50 'round_nearest';
  5         41  
80             *_divrem = \&Math::PlanePath::_divrem;
81              
82 5     5   1126 use Math::PlanePath::CoprimeColumns;
  5         10  
  5         275  
83             *_coprime = \&Math::PlanePath::CoprimeColumns::_coprime;
84              
85              
86             # uncomment this to run the ### lines
87             #use Smart::Comments;
88              
89 5     5   38 use constant class_x_negative => 0;
  5         68  
  5         271  
90 5     5   28 use constant class_y_negative => 0;
  5         9  
  5         245  
91 5     5   39 use constant x_minimum => 1;
  5         11  
  5         261  
92 5     5   32 use constant y_minimum => 1;
  5         8  
  5         248  
93 5     5   34 use constant gcdxy_maximum => 1; # no common factor
  5         12  
  5         387  
94              
95 5         1299 use constant parameter_info_array =>
96             [ { name => 'pairs_order',
97             display => 'Pairs Order',
98             type => 'enum',
99             default => 'rows',
100             choices => ['rows','rows_reverse','diagonals_down','diagonals_up'],
101             choices_display => ['Rows',
102             'Rows Reverse',
103             'Diagonals Down',
104             'Diagonals Up'],
105             description => 'Order in the i,j pairs.',
106 5     5   34 } ];
  5         9  
107              
108             sub absdy_minimum {
109 0     0 1 0 my ($self) = @_;
110 0 0       0 return ($self->{'pairs_order'} eq 'diagonals_down'
111             ? 1
112             : 0);
113             }
114              
115             {
116             my %dir_minimum_dxdy
117             = (rows => [1,0], # N=4 to N=5 horiz
118             rows_reverse => [1,0], # N=1 to N=2 horiz
119             diagonals_down => [0,1], # N=1 to N=2 vertical, nothing less
120             diagonals_up => [1,0], # N=4 to N=5 horiz
121             );
122             sub dir_minimum_dxdy {
123 0     0 1 0 my ($self) = @_;
124 0         0 return @{$dir_minimum_dxdy{$self->{'pairs_order'}}};
  0         0  
125             }
126             }
127             {
128             my %dir_maximum_dxdy
129             = (rows => [1,-1], # N=2 to N=3 SE diagonal
130             rows_reverse => [2,-1], # N=3 to N=4 dX=2,dY=-1
131             diagonals_down => [1,-1], # N=5 to N=6 SE diagonal
132             diagonals_up => [2,-1], # N=9 to N=10 dX=2,dY=-1
133             );
134             sub dir_maximum_dxdy {
135 0     0 1 0 my ($self) = @_;
136 0         0 return @{$dir_maximum_dxdy{$self->{'pairs_order'}}};
  0         0  
137             }
138             }
139              
140             #------------------------------------------------------------------------------
141              
142             # all rationals X,Y >= 1 no common factor
143 5     5   2935 use Math::PlanePath::DiagonalRationals;
  5         14  
  5         7073  
144             *xy_is_visited = Math::PlanePath::DiagonalRationals->can('xy_is_visited');
145              
146             sub new {
147 10     10 1 1184 my $self = shift->SUPER::new(@_);
148              
149 10   100     68 my $pairs_order = ($self->{'pairs_order'} ||= 'rows');
150             (($self->{'pairs_order_n_to_xy'}
151             = $self->can("_pairs_order__${pairs_order}__n_to_xy"))
152 10 50 33     109 && ($self->{'pairs_order_xygr_to_n'}
153             = $self->can("_pairs_order__${pairs_order}__xygr_to_n")))
154             or croak "Unrecognised pairs_order: ",$pairs_order;
155              
156 10         34 return $self;
157             }
158              
159             sub n_to_xy {
160 237     237 1 34251 my ($self, $n) = @_;
161             ### GcdRationals n_to_xy(): "$n"
162              
163 237 50       582 if ($n < 1) { return; }
  0         0  
164 237 50       694 if (is_infinite($n)) { return ($n,$n); }
  0         0  
165              
166             # what to do for fractional $n?
167             {
168 237         820 my $int = int($n);
  237         372  
169 237 50       453 if ($n != $int) {
170             ### frac ...
171 0         0 my $frac = $n - $int; # inherit possible BigFloat/BigRat
172 0         0 my ($x1,$y1) = $self->n_to_xy($int);
173 0         0 my ($x2,$y2) = $self->n_to_xy($int+1);
174 0         0 my $dx = $x2-$x1;
175 0         0 my $dy = $y2-$y1;
176 0         0 return ($frac*$dx + $x1, $frac*$dy + $y1);
177             }
178 237         404 $n = $int;
179             }
180              
181 237         546 my ($x,$y) = $self->{'pairs_order_n_to_xy'}->($n);
182              
183             # if ($self->{'pairs_order'} eq 'rows'
184             # || $self->{'pairs_order'} eq 'rows_reverse') {
185             # $y = int((sqrt(8*$n-7) + 1) / 2);
186             # $x = $n - ($y - 1)*$y/2;
187             #
188             # if ($self->{'pairs_order'} eq 'rows_reverse') {
189             # $x = $y - ($x-1);
190             # }
191             #
192             # # require Math::PlanePath::PyramidRows;
193             # # my ($x,$y) = Math::PlanePath::PyramidRows->new(step=>1)->n_to_xy($n);
194             # # $x+=1;
195             # # $y+=1;
196             #
197             # } else {
198             # require Math::PlanePath::DiagonalsOctant;
199             # ($x,$y) = Math::PlanePath::DiagonalsOctant->new->n_to_xy($n);
200             # if ($self->{'pairs_order'} eq 'diagonals_up') {
201             # my $d = $x+$y; # top 0,d measure diag down by x
202             # my $e = int($d/2); # end e,d-e
203             # ($x,$y) = ($e-$x, $d - ($e-$x));
204             # }
205             # $x+=1;
206             # $y+=1;
207             # }
208             ### triangle: "$x,$y"
209              
210 237         1367 my $gcd = _gcd($x,$y);
211 237         589 $x /= $gcd;
212 237         474 $y /= $gcd;
213              
214             ### $gcd
215             ### reduced: "$x,$y"
216             ### push out to x: $x + ($gcd-1)*$y
217              
218 237         797 return ($x + ($gcd-1)*$y, $y);
219             }
220              
221             sub _pairs_order__rows__n_to_xy {
222 102     102   1787 my ($n) = @_;
223 102         276 my $y = int( (_sqrtint(8*$n-7) + 1) / 2 );
224 102         3369 return ($n - ($y-1)*$y/2,
225             $y);
226             }
227             sub _pairs_order__rows_reverse__n_to_xy {
228 65     65   1549 my ($n) = @_;
229 65         171 my $y = int( (_sqrtint(8*$n-7) + 1) / 2 );
230 65         166 return ($y*($y+1)/2 + 1 - $n,
231             $y);
232             }
233             sub _pairs_order__diagonals_down__n_to_xy {
234 80     80   2997 my ($n) = @_;
235 80         195 my $d = _sqrtint($n-1); # eg. N=10 d=3
236 80         143 $n -= $d*($d+1); # eg. d=3 subtract 12
237 80 100       145 if ($n > 0) {
238 43         108 return ($n,
239             2 - $n + 2*$d);
240             } else {
241 37         86 return ($n + $d,
242             1 - $n + $d);
243             }
244             }
245             sub _pairs_order__diagonals_up__n_to_xy {
246 80     80   3028 my ($n) = @_;
247 80         203 my $d = _sqrtint($n-1);
248 80         136 $n -= $d*($d+1);
249 80 100       153 if ($n > 0) {
250 43         111 return (-$n + $d + 2,
251             $n + $d);
252             } else {
253 37         86 return (1 - $n,
254             $n + 2*$d);
255             }
256             }
257              
258              
259             # X=(g-1)*v+u
260             # Y=v
261             # u = x % y
262             # i = u*g
263             # = (x % y)*g
264             # = (x % y)*(floor(x/y)+1)
265             #
266             # Better:
267             # g-1 = floor(x/y)
268             # Y = j/g
269             # X = ((g-1)*j + i)/g
270             # j = Y*g
271             # (g-1)*j + i = X*g
272             # i = X*g - (g-1)*j
273             # = X*g - (g-1)*Y*g
274             # N = i + j*(j-1)/2
275             # = X*g - (g-1)*Y*g + Y*g*(Y*g-1)/2
276             # = X*g + Y*g * (-(g-1) + (Y*g-1)/2) # but Y*g-1 may be odd
277             # = X*g + Y*g * (Y*g-1 - (2g-2))/2
278             # = X*g + Y*g * (Y*g-1 - 2g + 2))/2
279             # = X*g + Y*g * (Y*g - 2g + 1))/2
280             # = X*g + Y*g * ((Y-2)*g + 1) / 2
281             # = g * [ X + Y*((Y-2)*g + 1) / 2 ]
282             #
283             # N = X*g - (g-1)*Y*g + Y*g*(Y*g-1)/2
284             # = [ 2*X*g - 2*(g-1)*Y*g + Y*g*(Y*g-1) ] / 2
285             # = [ 2*X - 2*(g-1)*Y + Y*(Y*g-1) ] * g / 2
286             # = [ 2*X + Y*(- 2*(g-1) + (Y*g-1)) ] * g / 2
287             # = [ 2*X + Y*(-2g + 2 + Y*g - 1) ] * g / 2
288             # = [ 2*X + Y*((Y-2)*g + 1) ] * g / 2
289             # = X*g + [(Y-2)*g + 1]*Y*g/2
290             #
291             # if Y and g both odd then (Y-2)*g+1 is odd+1 so even
292              
293             # q=int(x/y)
294             # x = qy+r qy=x-r
295             # r = x % y
296             # g-1 = q
297             # g = q+1
298             # g*y = (q+1)*y
299             # = q*y + y
300             # = x-r + y
301             #
302             # N = X*g + Y*g * ((Y-2)*g + 1) / 2
303             # = X*g + (X+Y-r) * ((Y-2)*g + 1) / 2
304             # = X*g + (X+Y-r) * ((g*Y-2*g + 1) / 2
305             # = X*g + (X+Y-r) * (((X+Y-r) - 2*g + 1) / 2
306             # ... not much better
307              
308             sub xy_to_n {
309 896     896 1 16088 my ($self, $x, $y) = @_;
310 896         1742 $x = round_nearest ($x);
311 896         1716 $y = round_nearest ($y);
312             ### GcdRationals xy_to_n(): "$x,$y"
313              
314 896 50       1686 if (is_infinite($x)) { return $x; }
  0         0  
315 896 50       9812 if (is_infinite($y)) { return $y; }
  0         0  
316 896 100 33     10818 if ($x < 1 || $y < 1 || ! _coprime($x,$y)) {
      66        
317 259         573 return undef;
318             }
319              
320 637         3051 my ($p,$r) = _divrem ($x,$y);
321             ### $x
322             ### $y
323             ### $p
324             ### $r
325 637         1453 return $self->{'pairs_order_xygr_to_n'}->($x,$y,$p+1,$r);
326              
327              
328             # my $g = int($x/$y) + 1;
329             # ### g: "$g"
330             # ### halve: ''.$y*(($y-2)*$g + 1)
331             # return $self->{'pairs_order_xygr_to_n'}->($x,$y,$g);
332             }
333              
334             sub _pairs_order__rows__xygr_to_n {
335 590     590   5183 my ($x,$y,$g,$r) = @_;
336             ### j: $x+$y-$r
337             ### i: $g*$r
338 590         828 $x += $y;
339 590         2571 $x -= $r; # j=X+Y-r
340 590         5028 return $x*($x-1)/2 + $g*$r; # i=g*r
341             }
342              
343             # i = X*g - (g-1)*g*Y
344             # = X*g - (g-1)*(X+Y-r)
345             # = X*g - g*(X+Y-r) + *(X+Y-r)
346             # = X*g - g*X - g*Y + g*r + (X+Y-r)
347             # = X*g - g*X - (X+Y-r) + g*r + (X+Y-r)
348             # = g*r
349             #
350             # N = j-i+1 + j*(j-1)/2
351             # = [2j-2i + 2 + $j*($j-1)] / 2
352             # = [-2i + 2 + 2j+ j*(j-1)] / 2
353             # = [-2i + 2 + j*(j-1+2)] / 2
354             # = [-2i + 2 + j*(j+1)] / 2
355             # = 1-i + j*(j+1)/2
356             #
357             sub _pairs_order__rows_reverse__xygr_to_n {
358 65     65   1050 my ($x,$y,$g,$r) = @_;
359 65         93 $y += $x;
360 65         91 $y -= $r; # j = X+Y-r
361 65 100       115 if ($r == 0) {
362             # Case r=0 which is Y=1 becomes i=0 and that doesn't reverse to the
363             # correct place by j-i+1. Can either set $r=1,$g+=1 or leave $r==0
364             # alone and adjust $y.
365 10         16 $y -= 2;
366             }
367 65         163 return $y*($y+1)/2 - $r*$g + 1;
368             }
369              
370             # d = (i-1)+(j-1)+1
371             # = i+j-1
372             # = rg + X+Y-r - 1
373             # = X+Y + r*(g-1) - 1
374             # if r==0 Y==1 then r=1 g=X-1
375             # i = r*g = X-1
376             # j = X+Y-r = X+1-1 = X-1
377             # d = i+j-1
378             # = 2X-2
379             # N = (d*d - (d%2))/4 + X-1
380             # = ((2X-2)*(2X-2) - 0)/4 + X-1
381             # = (X-1)^2 + X-1
382             #
383             sub _pairs_order__diagonals_down__xygr_to_n {
384 80     80   1845 my ($x,$y,$g,$r) = @_;
385              
386 80         138 $y += $x + $r*($g-1) - 1; # d=X+Y + r*(g-1) - 1
387 80 100       155 if ($r == 0) {
388 7         29 $y *= 2; # d=2*g-2
389             }
390 80         188 return ($y*$y - ($y % 2))/4 + $r*$g;
391             }
392             sub _pairs_order__diagonals_up__xygr_to_n {
393 80     80   1846 my ($x,$y,$g,$r) = @_;
394              
395 80         129 $y += $x + $r*($g-1); # d=X+Y + r*(g-1)
396 80 100       150 if ($r == 0) {
397 7         14 $y = 2*$x - 1;
398             }
399 80         203 return ($y*$y - ($y % 2))/4 - $r*$g + 1;
400             }
401              
402              
403             # increase in rows, so right column
404             # in column increase within g wedge, then drop
405             #
406             # int(x2/y2) is slope of top of the wedge containing x2,y2
407             # g = int(x2/y2)+1 is the slope of the bottom of that wedge
408             # yw = floor(x2 / g) is the Y of that bottom
409             # N at x2,yw,g+1 is the top of the wedge underneath, bigger g smaller y
410             # or x2,y2,g is the top-right corner
411             #
412             # Eg.
413             # x=19 y=2 to 4
414             # g=int(19/4)+1=5
415             # yw=int(19/5)=3
416             # N(19,3,6)=
417             #
418             # at X=Y+1 g=2
419             # nhi = (y*((y-2)*g + 1) / 2 + x)*g
420             # = (y*((y-2)*2 + 1) / 2 + y+1)*2
421             # = (y*(2y-4 + 1) / 2 + y+1)*2
422             # = (y*(2y-3) / 2 + y+1)*2
423             # = y*(2y-3) + 2y+2
424             # = 2y^2 - 3y + 2y + 2
425             # = 2y^2 - y + 2
426             # = y*(2y-1) + 2
427              
428             # 11 12 13 14 47 49 51 53 108 111 114 117 194 198 202 206
429             # 7 9 30 34 69 75 124 132 195 205
430             # 4 5 17 19 39 42 70 74 110 115 159 165 217
431             # 2 8 18 32 50 72 98 128 162 200
432             # 1 3 6 10 15 21 28 36 45 55 66 78 91 105 120 136 153 171 190
433              
434             # 206=20*19/2+16 i=16,j=20 gcd=4
435             # 19,5 is slope=floor(19/5)=3 so g=4
436             #
437             # 205=20*19/2+15 i=15,j=20 gcd=5
438             # 19,4 is slope=floor(19/4)=4 so g=5
439             #
440             # 217=21*20/2 + 7, i=21,j=7 gcd=7
441             # 19,3 is slope=floor(19/3)=6 so g=7
442              
443             # not exact
444             sub rect_to_n_range {
445 44     44 1 674 my ($self, $x1,$y1, $x2,$y2) = @_;
446             ### rect_to_n_range(): "$x1,$y1 $x2,$y2"
447              
448 44         98 $x1 = round_nearest ($x1);
449 44         84 $y1 = round_nearest ($y1);
450 44         78 $x2 = round_nearest ($x2);
451 44         82 $y2 = round_nearest ($y2);
452              
453 44 50       81 ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
454 44 50       82 ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
455             ### $x2
456             ### $y2
457              
458 44 50 33     146 if ($x2 < 1 || $y2 < 1) {
459 0         0 return (1, 0); # outside quadrant
460             }
461              
462 44 100       77 if ($x1 < 1) { $x1 = 1; }
  1         2  
463 44 100       72 if ($y1 < 1) { $y1 = 1; }
  1         2  
464              
465 44 50       94 if ($self->{'pairs_order'} =~ /^diagonals/) {
466 0         0 my $d = $x2 + max($x2,$y2);
467 0         0 return (1, int($d*($d+($d%2)) / 4)); # N end of diagonal d
468             }
469              
470 44         60 my $nhi;
471             {
472 44         57 my $c = max($x2,$y2);
  44         96  
473 44         87 $nhi = _pairs_order__rows__xygr_to_n($c,$c,2,0);
474              
475             # my $rev = ($self->{'pairs_order'} eq 'rows_reverse');
476             # my $slope = int($x2/$y2);
477             # my $g = $slope + 1;
478             #
479             # # within top row
480             # {
481             # my $x;
482             # if ($rev) {
483             # if ($slope > 0) {
484             # $x = max ($x1, $y2*$slope); # left-most within this wedge
485             # } else {
486             # $x = $x1; # top-left corner
487             # }
488             # } else {
489             # # pairs_order=rows
490             # $x = $x2; # top-right corner
491             # }
492             # $nhi = $self->{'pairs_order_xygr_to_n'}->($x, $y2, $g, 0);
493             #
494             # ### $slope
495             # ### $g
496             # ### x for hi: $x
497             # ### nhi for x,y2: $nhi
498             # }
499             #
500             # # within x2 column, top of wedge below
501             # #
502             # my $yw = int(($x2+$g-1) / $g); # rounded up
503             # if ($yw >= $y1) {
504             # $nhi = max ($nhi, $self->{'pairs_order_xygr_to_n'}->($x2,$yw,$g+1,0));
505             #
506             # ### $yw
507             # ### nhi_wedge: $self->{'pairs_order_xygr_to_n'}->($x2,$yw,$g+1,0)
508             # }
509             # my $yw = int($x2 / $g) - ($g==1); # below X=Y diagonal when g==1
510             # if ($yw >= $y1) {
511             # $g = int($x2/$yw) + 1; # perhaps went across more than one wedge
512             # $nhi = max ($nhi,
513             # ($yw*(($yw-2)*($g+1) + 1) / 2 + $x2)*($g+1));
514             # ### $yw
515             # ### nhi_wedge: ($yw*(($yw-2)*($g+1) + 1) / 2 + $x2)*($g+1)
516             # }
517             }
518              
519 44         75 my $nlo;
520             {
521 44         52 $nlo = _pairs_order__rows__xygr_to_n(1,$x1, 1, $x1-1);
  44         86  
522              
523             # my $g = int($x1/$y1) + 1;
524             # $nlo = $self->{'pairs_order_xygr_to_n'}->($x1,$y1,$g,0);
525             #
526             # ### glo: $g
527             # ### $nlo
528             #
529             # if ($g > 1) {
530             # my $yw = max (int($x1 / $g),
531             # 1);
532             # ### $yw
533             # if ($yw <= $y2) {
534             # $g = int($x1/$yw); # no +1, and perhaps up across more than one wedge
535             # $nlo = min ($nlo, $self->{'pairs_order_xygr_to_n'}->($x1,$yw,$g,0));
536             # ### glo_wedge: $g
537             # ### nlo_wedge: $self->{'pairs_order_xygr_to_n'}->($x1,$yw,$g,0)
538             # }
539             # }
540             # if ($nlo < 1) {
541             # $nlo = 1;
542             # }
543             }
544              
545             ### $nhi
546             ### $nlo
547 44         96 return ($nlo, $nhi);
548             }
549              
550             sub _gcd {
551 259     259   53801 my ($x, $y) = @_;
552             #### _gcd(): "$x,$y"
553              
554             # bgcd() available in even the earliest Math::BigInt
555 259 100 66     846 if ((ref $x && $x->isa('Math::BigInt'))
      33        
      66        
556             || (ref $y && $y->isa('Math::BigInt'))) {
557 3         14 return Math::BigInt::bgcd($x,$y);
558             }
559              
560 256         349 $x = abs(int($x));
561 256         331 $y = abs(int($y));
562 256 50       463 unless ($x > 0) {
563 0         0 return $y; # gcd(0,y)=y for y>=0, giving gcd(0,0)=0
564             }
565 256 100       497 if ($y > $x) {
566 192         318 $y %= $x;
567             }
568 256         339 for (;;) {
569             ### assert: $x >= 1
570              
571 368 100       638 if ($y <= 1) {
572 256 100       585 return ($y == 0
573             ? $x # gcd(x,0)=x
574             : 1); # gcd(x,1)=1
575             }
576 112         199 ($x,$y) = ($y, $x % $y);
577             }
578             }
579              
580              
581              
582             # # old code, rows only ...
583             # sub rect_to_n_range {
584             # my ($self, $x1,$y1, $x2,$y2) = @_;
585             # ### rect_to_n_range(): "$x1,$y1 $x2,$y2"
586             #
587             # $x1 = round_nearest ($x1);
588             # $y1 = round_nearest ($y1);
589             # $x2 = round_nearest ($x2);
590             # $y2 = round_nearest ($y2);
591             #
592             # ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
593             # ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
594             # ### $x2
595             # ### $y2
596             #
597             # if ($x2 < 1 || $y2 < 1) {
598             # return (1, 0); # outside quadrant
599             # }
600             #
601             # if ($x1 < 1) { $x1 = 1; }
602             # if ($y1 < 1) { $y1 = 1; }
603             #
604             # my $g = int($x2/$y2) + 1;
605             # my $nhi = ($y2*(($y2-2)*$g + 1) / 2 + $x2)*$g;
606             # ### ghi: $g
607             # ### $nhi
608             #
609             # my $yw = int($x2 / $g) - ($g==1); # below X=Y diagonal when g==1
610             # if ($yw >= $y1) {
611             # $g = int($x2/$yw) + 1; # perhaps went across more than one wedge
612             # $nhi = max ($nhi,
613             # ($yw*(($yw-2)*($g+1) + 1) / 2 + $x2)*($g+1));
614             # ### $yw
615             # ### nhi_wedge: ($yw*(($yw-2)*($g+1) + 1) / 2 + $x2)*($g+1)
616             # }
617             #
618             # $g = int($x1/$y1) + 1;
619             # my $nlo = ($y1*(($y1-2)*$g + 1) / 2 + $x1)*$g;
620             #
621             # ### glo: $g
622             # ### $nlo
623             #
624             # if ($g > 1) {
625             # $yw = max (int($x1 / $g),
626             # 1);
627             # ### $yw
628             # if ($yw <= $y2) {
629             # $g = int($x1/$yw); # no +1, and perhaps up across more than one wedge
630             # $nlo = min ($nlo,
631             # ($yw*(($yw-2)*$g + 1) / 2 + $x1)*$g);
632             # ### glo_wedge: $g
633             # ### nlo_wedge: ($yw*(($yw-2)*$g + 1) / 2 + $x1)*$g
634             # }
635             # }
636             #
637             # return ($nlo, $nhi);
638             # }
639              
640              
641             1;
642             __END__