File Coverage

blib/lib/Math/PlanePath/CoprimeColumns.pm
Criterion Covered Total %
statement 129 144 89.5
branch 38 48 79.1
condition 28 41 68.2
subroutine 24 25 96.0
pod 5 5 100.0
total 224 263 85.1


line stmt bran cond sub pod time code
1             # Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 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             # math-image --path=CoprimeColumns --all --scale=10
20             # math-image --path=CoprimeColumns --output=numbers --all
21              
22             package Math::PlanePath::CoprimeColumns;
23 10     10   27448 use 5.004;
  10         54  
24 10     10   58 use strict;
  10         19  
  10         312  
25              
26 10     10   75 use vars '$VERSION', '@ISA', '@_x_to_n';
  10         18  
  10         776  
27             $VERSION = 127;
28 10     10   2006 use Math::PlanePath;
  10         17  
  10         384  
29             @ISA = ('Math::PlanePath');
30              
31             use Math::PlanePath::Base::Generic
32 10         496 'is_infinite',
33 10     10   57 'round_nearest';
  10         34  
34              
35             # uncomment this to run the ### lines
36             # use Smart::Comments;
37              
38              
39 10     10   59 use constant default_n_start => 0;
  10         20  
  10         651  
40 10     10   59 use constant class_x_negative => 0;
  10         27  
  10         480  
41 10     10   58 use constant class_y_negative => 0;
  10         20  
  10         500  
42 10     10   54 use constant n_frac_discontinuity => .5;
  10         20  
  10         475  
43              
44 10     10   59 use constant x_minimum => 1;
  10         22  
  10         412  
45 10     10   55 use constant y_minimum => 1;
  10         27  
  10         571  
46 10     10   65 use constant diffxy_minimum => 0; # octant Y<=X so X-Y>=0
  10         19  
  10         477  
47 10     10   60 use constant gcdxy_maximum => 1; # no common factor
  10         21  
  10         493  
48              
49 10     10   58 use constant dx_minimum => 0;
  10         20  
  10         468  
50 10     10   60 use constant dx_maximum => 1;
  10         18  
  10         567  
51 10     10   62 use constant dir_maximum_dxdy => (1,-1); # South-East
  10         26  
  10         889  
52              
53 10         9972 use constant parameter_info_array =>
54             [
55             { name => 'direction',
56             share_key => 'direction_updown',
57             display => 'Direction',
58             type => 'enum',
59             default => 'up',
60             choices => ['up','down'],
61             choices_display => ['Down','Up'],
62             description => 'Number points upwards or downwards in the columns.',
63             },
64             Math::PlanePath::Base::Generic::parameter_info_nstart0(),
65 10     10   68 ];
  10         32  
66              
67             #------------------------------------------------------------------------------
68              
69             sub new {
70 8     8 1 2499 my $self = shift->SUPER::new (@_);
71 8   100     60 $self->{'direction'} ||= 'up';
72 8 100       30 if (! defined $self->{'n_start'}) {
73 7         41 $self->{'n_start'} = $self->default_n_start;
74             }
75 8         21 return $self;
76             }
77              
78              
79             # shared with DiagonalRationals
80             @_x_to_n = (0,0,1);
81             sub _extend {
82             ### _extend(): $#_x_to_n
83 298     298   475 my $x = $#_x_to_n;
84 298         531 push @_x_to_n, $_x_to_n[$x] + _totient($x);
85              
86             # if ($x > 2) {
87             # if (($x & 3) == 2) {
88             # $x >>= 1;
89             # $next_n += $_x_to_n[$x] - $_x_to_n[$x-1];
90             # } else {
91             # $next_n +=
92             # }
93             # }
94             ### last x: $#_x_to_n
95             ### second last: $_x_to_n[$#_x_to_n-2]
96             ### last: $_x_to_n[$#_x_to_n-1]
97             ### diff: $_x_to_n[$#_x_to_n-1] - $_x_to_n[$#_x_to_n-2]
98             ### totient of: $#_x_to_n - 2
99             ### totient: _totient($#_x_to_n-2)
100             ### assert: $_x_to_n[$#_x_to_n-1] - $_x_to_n[$#_x_to_n-2] == _totient($#_x_to_n-2)
101             }
102              
103             sub n_to_xy {
104 46     46 1 2636 my ($self, $n) = @_;
105             ### CoprimeColumns n_to_xy(): $n
106              
107 46         80 $n = $n - $self->{'n_start'}; # to N=0 basis, and warn on undef
108              
109             # $n<-0.5 is ok for Math::BigInt circa Perl 5.12, it seems
110 46 100       679 if (2*$n < -1) {
111 2         525 return;
112             }
113 44 50       350 if (is_infinite($n)) {
114 0         0 return ($n,$n);
115             }
116              
117 44         266 my $frac;
118             {
119 44         57 my $int = int($n);
  44         78  
120 44         105 $frac = $n - $int; # -.5 <= $frac < 1
121 44         120 $n = $int; # BigFloat int() gives BigInt, use that
122              
123 44 50       94 if (2*$frac >= 1) {
124 0         0 $frac--;
125 0         0 $n += 1;
126             # now -.5 <= $frac < .5
127             }
128             ### $n
129             ### $frac
130             ### assert: 2*$frac >= -1
131             ### assert: 2*$frac < 1
132             }
133              
134 44         314 my $x = 1;
135 44         69 for (;;) {
136 252         485 while ($x > $#_x_to_n) {
137 13         26 _extend();
138             }
139 252 100       440 if ($_x_to_n[$x] > $n) {
140 44         129 $x--;
141 44         70 last;
142             }
143 208         408 $x++;
144             }
145 44         70 $n -= $_x_to_n[$x];
146             ### $x
147             ### n base: $_x_to_n[$x]
148             ### n next: $_x_to_n[$x+1]
149             ### remainder: $n
150              
151 44         224 my $y = 1;
152 44         64 for (;;) {
153 104 100       171 if (_coprime($x,$y)) {
154 90 100       170 if (--$n < 0) {
155 44         258 last;
156             }
157             }
158 60 50       131 if (++$y >= $x) {
159             ### oops, not enough in this column ...
160 0         0 return;
161             }
162             }
163              
164 44         80 $y += $frac;
165 44 100 100     355 if ($x >= 2 && $self->{'direction'} eq 'down') {
166 27         40 $y = $x - $y;
167             }
168 44         128 return ($x, $y);
169             }
170              
171             sub xy_is_visited {
172 0     0 1 0 my ($self, $x, $y) = @_;
173 0         0 $x = round_nearest ($x);
174 0         0 $y = round_nearest ($y);
175 0 0 0     0 if ($x < 1
      0        
      0        
176             || $y < 1
177             || $y >= $x+($x==1) # Y
178             || ! _coprime($x,$y)) {
179 0         0 return 0;
180             }
181 0         0 return 1;
182             }
183              
184             sub xy_to_n {
185 11283     11283 1 55184 my ($self, $x, $y) = @_;
186             ### CoprimeColumns xy_to_n(): "$x,$y"
187 11283         20535 $x = round_nearest ($x);
188 11283         20383 $y = round_nearest ($y);
189 11283 50       20438 if (is_infinite($x)) { return $x; }
  0         0  
190 11283 50       22588 if (is_infinite($y)) { return $y; }
  0         0  
191 11283 100 100     47991 if ($x < 1
      100        
      100        
192             || $y < 1
193             || $y >= $x+($x==1) # Y
194             || ! _coprime($x,$y)) {
195 8231         16809 return undef;
196             }
197              
198 3052 100 100     10415 if ($x >= 2 && $self->{'direction'} eq 'down') {
199 2034         3082 $y = $x - $y;
200             }
201              
202 3052         6508 while ($#_x_to_n < $x) {
203 0         0 _extend();
204             }
205 3052         5090 my $n = $_x_to_n[$x];
206             ### base n: $n
207 3052 100       5392 if ($y != 1) {
208 2935         5761 foreach my $i (1 .. $y-1) {
209 130132 100       190609 if (_coprime($x,$i)) {
210 68320         109912 $n += 1;
211             }
212             }
213             }
214 3052         7482 return $n + $self->{'n_start'};
215             }
216              
217             # Asymptotically
218             # phisum(x) ~ 1/(2*zeta(2)) * x^2 + O(x ln x)
219             # = 3/pi^2 * x^2 + O(x ln x)
220             # or by Walfisz
221             # phisum(x) ~ 3/pi^2 * x^2 + O(x * (ln x)^(2/3) * (ln ln x)^4/3)
222             #
223             # but want an upper bound, so that for a given X at least enough N is
224             # covered ...
225             #
226             # Note: DiagonalRationals depends on this working only to column resolution.
227             # not exact
228             sub rect_to_n_range {
229 24     24 1 2039 my ($self, $x1,$y1, $x2,$y2) = @_;
230             ### CoprimeColumns rect_to_n_range(): "$x1,$y1 $x2,$y2"
231              
232 24 100       69 ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
233 24 50       42 ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
234 24         62 $x2 = round_nearest($x2);
235 24         47 $y2 = round_nearest($y2);
236             ### rounded ...
237             ### $x2
238             ### $y2
239              
240 24 50 33     108 if ($x2 < 1 || $y2 < 1
      33        
241             # bottom right corner above X=Y diagonal, except X=1,Y=1 included
242             || ($y1 >= $x2 + ($x2 == 1))) {
243             ### outside ...
244 0         0 return (1, 0);
245             }
246 24 50       48 if (is_infinite($x2)) {
247 0         0 return ($self->{'n_start'}, $x2);
248             }
249              
250 24         67 while ($#_x_to_n <= $x2) {
251 93         170 _extend();
252             }
253              
254             ### rect use xy_to_n at: "x=".($x2+1)." y=1"
255 24 100       46 if ($x1 < 0) { $x1 = 0; }
  2         6  
256             return ($_x_to_n[$x1] + $self->{'n_start'},
257 24         84 $_x_to_n[$x2+1] - 1 + $self->{'n_start'});
258              
259             # asympototically ?
260             # return ($self->{'n_start'}, $self->{'n_start'} + .304*$x2*$x2 + 20);
261             }
262              
263             # A000010
264             sub _totient {
265 404     404   1096 my ($x) = @_;
266 404   100     2015 my $count = (1 # y=1 always
      100        
      100        
267             + ($x > 2 && ($x&1)) # y=2 if $x odd
268             + ($x > 3 && ($x % 3) != 0) # y=3
269             + ($x > 4 && ($x&1)) # y=4 if $x odd
270             );
271 404         833 for (my $y = 5; $y < $x; $y++) {
272 28231         43303 $count += _coprime($x,$y);
273             }
274 404         1232 return $count;
275             }
276              
277             # code here only uses X>=Y but allow for any X,Y>=0 for elsewhere
278             sub _coprime {
279 171859     171859   278092 my ($x, $y) = @_;
280             #### _coprime(): "$x,$y"
281              
282 171859 100       291413 if ($y > $x) {
283 424 100       1079 if ($x <= 1) {
284             ### result yes ...
285 56         687 return 1;
286             }
287 368         1179 $y %= $x;
288             }
289              
290 171803         217614 for (;;) {
291 697175 100       1121233 if ($y <= 1) {
292             ### result: ($y == 1)
293 171803         386642 return ($y == 1);
294             }
295 525372         801427 ($x,$y) = ($y, $x % $y);
296             }
297             }
298              
299             1;
300             __END__