File Coverage

blib/lib/Math/PlanePath/DivisibleColumns.pm
Criterion Covered Total %
statement 112 172 65.1
branch 26 66 39.3
condition 7 39 17.9
subroutine 20 24 83.3
pod 7 7 100.0
total 172 308 55.8


line stmt bran cond sub pod time code
1             # Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 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             # A006218 - cumulative count of divisors
20             #
21             # Dirichlet:
22             # n * (log(n) + 2*gamma - 1) + O(sqrt(n)) gamma=0.57721... Euler-Mascheroni
23             #
24             # n * (log(n) + 2*gamma - 1) + O(log(n)*n^(1/3))
25             #
26             # Chandrasekharan: bounds
27             # n log(n) + (2 gamma - 1) n - 4 sqrt(n) - 1
28             # <= a(n) <=
29             # n log(n) + (2 gamma - 1) n + 4 sqrt(n)
30             #
31             # a(n)=2 * sum[ i=1 to floor(sqrt(n)) of floor(n/i) ] - floor(sqrt(n))^2
32             #
33             # cf A003988,A010766 - triangle with values floor(i/j)
34             #
35             # http://mathworld.wolfram.com/DirichletDivisorProblem.html
36             #
37             # compile-command: "math-image --path=DivisibleColumns --all"
38             #
39             # math-image --path=DivisibleColumns --output=numbers --all
40             #
41              
42             package Math::PlanePath::DivisibleColumns;
43 1     1   9087 use 5.004;
  1         9  
44 1     1   14 use strict;
  1         4  
  1         30  
45              
46 1     1   6 use vars '$VERSION', '@ISA';
  1         2  
  1         68  
47             $VERSION = 128;
48 1     1   697 use Math::PlanePath;
  1         3  
  1         50  
49             *_sqrtint = \&Math::PlanePath::_sqrtint;
50             @ISA = ('Math::PlanePath');
51              
52             use Math::PlanePath::Base::Generic
53 1         68 'is_infinite',
54 1     1   7 'round_nearest';
  1         2  
55              
56             # uncomment this to run the ### lines
57             # use Smart::Comments;
58              
59 1         51 use constant parameter_info_array =>
60             [ { name => 'divisor_type',
61             share_key => 'divisor_type_allproper',
62             display => 'Divisor Type',
63             type => 'enum',
64             choices => ['all','proper'],
65             default => 'all',
66             description => 'Divisor type, with "proper" meaning divisors d
67             },
68             # { name => 'direction',
69             # share_key => 'direction_updown',
70             # display => 'Direction',
71             # type => 'enum',
72             # default => 'up',
73             # choices => ['up','down'],
74             # choices_display => ['Down','Up'],
75             # description => 'Number points upwards or downwards in the columns.',
76             # },
77             Math::PlanePath::Base::Generic::parameter_info_nstart0(),
78 1     1   6 ];
  1         2  
79              
80 1     1   5 use constant default_n_start => 0;
  1         2  
  1         39  
81 1     1   6 use constant class_x_negative => 0;
  1         2  
  1         38  
82 1     1   5 use constant class_y_negative => 0;
  1         2  
  1         54  
83 1     1   6 use constant n_frac_discontinuity => .5;
  1         2  
  1         128  
84              
85             # X=2,Y=1 when proper
86             # X=1,Y=1 when not
87             sub x_minimum {
88 0     0 1 0 my ($self) = @_;
89 0 0       0 return ($self->{'proper'} ? 2 : 1);
90             }
91 1     1   8 use constant y_minimum => 1;
  1         2  
  1         92  
92              
93             sub diffxy_minimum {
94 0     0 1 0 my ($self) = @_;
95             # octant Y<=X so X-Y>=0
96 0 0       0 return ($self->{'proper'} ? 1 : 0);
97             }
98              
99 1     1   7 use constant dx_minimum => 0;
  1         2  
  1         67  
100 1     1   6 use constant dx_maximum => 1;
  1         2  
  1         46  
101 1     1   6 use constant dir_maximum_dxdy => (1,-1); # South-East
  1         1  
  1         1269  
102              
103              
104             #------------------------------------------------------------------------------
105              
106             sub new {
107 4     4 1 1799 my $self = shift->SUPER::new (@_);
108              
109 4   50     30 my $divisor_type = ($self->{'divisor_type'} ||= 'all');
110 4         13 $self->{'proper'} = ($divisor_type eq 'proper'); # bool
111              
112 4   50     19 $self->{'direction'} ||= 'up';
113              
114 4 100       15 if (! defined $self->{'n_start'}) {
115 3         13 $self->{'n_start'} = $self->default_n_start;
116             }
117 4         11 return $self;
118             }
119              
120             my @x_to_n = (0,0,1);
121             sub _extend {
122             ### _extend(): $#x_to_n
123 198     198   297 my $x = $#x_to_n;
124 198         394 push @x_to_n, $x_to_n[$x] + _count_divisors($x);
125              
126             # if ($x > 2) {
127             # if (($x & 3) == 2) {
128             # $x >>= 1;
129             # $next_n += $x_to_n[$x] - $x_to_n[$x-1];
130             # } else {
131             # $next_n +=
132             # }
133             # }
134             ### last x: $#x_to_n
135             ### second last: $x_to_n[$#x_to_n-2]
136             ### last: $x_to_n[$#x_to_n-1]
137             ### diff: $x_to_n[$#x_to_n-1] - $x_to_n[$#x_to_n-2]
138             ### divisors of: $#x_to_n - 2
139             ### divisors: _count_divisors($#x_to_n-2)
140             ### assert: $x_to_n[$#x_to_n-1] - $x_to_n[$#x_to_n-2] == _count_divisors($#x_to_n-2)
141             }
142              
143             sub n_to_xy {
144 0     0 1 0 my ($self, $n) = @_;
145             ### DivisibleColumns n_to_xy(): "$n"
146              
147 0         0 $n = $n - $self->{'n_start'}; # to N=0 basis, and warn on undef
148              
149             # $n<-0.5 works with Math::BigInt circa Perl 5.12, it seems
150 0 0       0 if ($n < -0.5) {
151 0         0 return;
152             }
153 0 0       0 if (is_infinite($n)) {
154 0         0 return ($n,$n);
155             }
156              
157 0         0 my $frac;
158             {
159 0         0 my $int = int($n);
  0         0  
160 0 0       0 if ($n == $int) {
161 0         0 $frac = 0;
162             } else {
163 0         0 $frac = $n - $int; # -.5 <= $frac < 1
164 0         0 $n = $int; # BigFloat int() gives BigInt, use that
165 0 0       0 if ($frac > .5) {
166 0         0 $frac--;
167 0         0 $n += 1;
168             # now -.5 <= $frac < .5
169             }
170             }
171             ### $n
172             ### n: "$n"
173             ### $frac
174             ### assert: $frac >= -.5
175             ### assert: $frac < .5
176             }
177 0   0     0 my $proper = $self->{'proper'} || 0; # cannot add false '' to BigInt
178             ### $proper
179              
180 0         0 my $x;
181 0 0       0 if ($proper) {
182 0         0 $x = 2;
183             ### proper adjusted n: $n
184             } else {
185 0         0 $x = 1;
186             }
187              
188 0         0 for (;;) {
189 0         0 while ($x > $#x_to_n) {
190 0         0 _extend();
191             }
192 0         0 $n += $proper;
193             ### consider: "n=$n x=$x x_to_n=".$x_to_n[$x]
194 0 0       0 if ($x_to_n[$x] > $n) {
195 0         0 $x--;
196 0         0 last;
197             }
198 0         0 $x++;
199             }
200 0         0 $n -= $x_to_n[$x];
201 0         0 $n -= $proper;
202             ### $x
203             ### x_to_n: $x_to_n[$x]
204             ### x_to_n next: $x_to_n[$x+1]
205             ### remainder: $n
206              
207 0         0 my $y = 1;
208 0         0 for (;;) {
209 0 0       0 unless ($x % $y) {
210 0 0       0 if (--$n < 0) {
211 0         0 return ($x, $frac + $y);
212             }
213             }
214 0 0       0 if (++$y > $x) {
215             ### oops, not enough in this column
216 0         0 return;
217             }
218             }
219             }
220              
221             # Feturn a count of the number of integers dividing $x, including 1 and $x
222             # itself. Cf Math::Factor::XS maybe.
223             sub _count_divisors {
224 705     705   153477 my ($x) = @_;
225 705         1067 my $ret = 1;
226 705 100       1626 unless ($x % 2) {
227 352         485 my $count = 1;
228 352         469 do {
229 692         976 $x /= 2;
230 692         1344 $count++;
231             } until ($x % 2);
232 352         527 $ret *= $count;
233             }
234 705         1515 my $limit = _sqrtint($x);
235 705         1653 for (my $d = 3; $d <= $limit; $d+=2) {
236 2007 100       4246 unless ($x % $d) {
237 414         547 my $count = 1;
238 414         530 do {
239 583         816 $x /= $d;
240 583         1088 $count++;
241             } until ($x % $d);
242 414         589 $limit = sqrt($x);
243 414         938 $ret *= $count;
244             }
245             }
246 705 100       1317 if ($x > 1) {
247 616         890 $ret *= 2;
248             }
249 705         1510 return $ret;
250             }
251              
252             sub xy_is_visited {
253 0     0 1 0 my ($self, $x, $y) = @_;
254 0         0 $x = round_nearest ($x);
255 0         0 $y = round_nearest ($y);
256             return ($y >= 1
257 0   0     0 && ($self->{'proper'}
258             ? $x >= 2 && $y <= int($x/2)
259             : $x >= 1 && $y <= $x)
260             && ($x%$y) == 0);
261             }
262              
263             sub xy_to_n {
264 400     400 1 1373 my ($self, $x, $y) = @_;
265             ### DivisibleColumns xy_to_n(): "$x,$y"
266              
267 400         776 $x = round_nearest ($x);
268 400         745 $y = round_nearest ($y);
269 400 50       755 if (is_infinite($x)) { return $x; }
  0         0  
270 400 50       847 if (is_infinite($y)) { return $y; }
  0         0  
271              
272 400         770 my $proper = $self->{'proper'};
273 400 50       692 if ($proper) {
274 0 0 0     0 if ($x < 2
      0        
      0        
275             || $y < 1
276             || $y > int($x/2)
277             || ($x%$y)) {
278 0         0 return undef;
279             }
280             } else {
281 400 50 33     1927 if ($x < 1
      33        
      33        
282             || $y < 1
283             || $y > $x
284             || ($x%$y)) {
285 0         0 return undef;
286             }
287             }
288              
289 400         878 while ($#x_to_n < $x) {
290 198         392 _extend();
291             }
292             ### x_to_n: $x_to_n[$x]
293              
294 400 50       761 my $n = $x_to_n[$x] - ($proper ? $x-1 : 1);
295             ### base n: $n
296              
297 400         778 for (my $i = 1+$proper; $i <= $y; $i++) {
298 20300 100       40083 unless ($x % $i) {
299 1298         2389 $n += 1;
300             }
301             }
302 400         825 return $n + $self->{'n_start'};
303             }
304              
305             # not exact
306             sub rect_to_n_range {
307 200     200 1 26007 my ($self, $x1,$y1, $x2,$y2) = @_;
308             ### DivisibleColumns rect_to_n_range(): "$x1,$y1 $x2,$y2"
309              
310 200         565 $x1 = round_nearest($x1);
311 200         399 $y1 = round_nearest($y1);
312 200         385 $x2 = round_nearest($x2);
313 200         365 $y2 = round_nearest($y2);
314              
315 200 50       451 ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
316 200 50       351 ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
317              
318             ### rounded ...
319             ### $x2
320             ### $y2
321              
322 200 50       431 if ($self->{'proper'}) {
323 0 0 0     0 if ($x2 < 2 # rect all negative
      0        
324             || $y2 < 1 # rect all negative
325             || 2*$y1 > $x2) { # rect all above X=2Y octant
326             ### outside proper divisors ...
327 0         0 return (1, 0);
328             }
329 0 0       0 if ($x1 < 2) { $x1 = 2; }
  0         0  
330             } else {
331 200 50 33     887 if ($x2 < 1 # rect all negative
      33        
332             || $y2 < 1 # rect all negative
333             || $y1 > $x2) { # rect all above X=Y diagonal
334             ### outside all divisors ...
335 0         0 return (1, 0);
336             }
337 200 50       353 if ($x1 < 1) { $x1 = 1; }
  0         0  
338             }
339 200 50       454 if (is_infinite($x2)) {
340 0         0 return ($self->{'n_start'}, $x2);
341             }
342              
343 200         352 my ($n_lo, $n_hi);
344 200 100       398 if ($x1 <= $#x_to_n) {
345 2         5 $n_lo = $x_to_n[$x1];
346             } else {
347 198         399 $n_lo = _count_divisors_cumulative($x1-1);
348             }
349 200 100       393 if ($x2 < $#x_to_n) {
350 1         3 $n_hi = $x_to_n[$x2+1];
351             } else {
352 199         301 $n_hi = _count_divisors_cumulative($x2);
353             }
354 200         436 $n_hi -= 1;
355              
356             ### rect at: "x=".($x2+1)." x_to_n=".($x_to_n[$x2+1]||'none')
357              
358 200 50       422 if ($self->{'proper'}) {
359 0         0 $n_lo -= $x1-1;
360 0         0 $n_hi -= $x2;
361             }
362             return ($n_lo + $self->{'n_start'},
363 200         482 $n_hi + $self->{'n_start'});
364             }
365              
366             # Return a total count of all the divisors of all the integers 1 to $x
367             # inclusive.
368             sub _count_divisors_cumulative {
369 904     904   155040 my ($x) = @_;
370 904         1293 my $total = 0;
371 904         1941 my $limit = _sqrtint($x);
372 904         1759 foreach my $i (1 .. $limit) {
373 10820         16319 $total += int($x/$i);
374             }
375 904         1833 return 2*$total - $limit*$limit;
376             }
377              
378             1;
379             __END__