File Coverage

blib/lib/Math/PlanePath/KochPeaks.pm
Criterion Covered Total %
statement 137 148 92.5
branch 26 34 76.4
condition 7 9 77.7
subroutine 26 27 96.3
pod 5 5 100.0
total 201 223 90.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             package Math::PlanePath::KochPeaks;
20 1     1   1113 use 5.004;
  1         2  
21 1     1   3 use strict;
  1         2  
  1         34  
22             #use List::Util 'max';
23             *max = \&Math::PlanePath::_max;
24              
25 1     1   4 use vars '$VERSION', '@ISA';
  1         2  
  1         63  
26             $VERSION = 127;
27 1     1   661 use Math::PlanePath;
  1         2  
  1         35  
28             @ISA = ('Math::PlanePath');
29              
30             use Math::PlanePath::Base::Generic
31 1         37 'is_infinite',
32 1     1   6 'round_nearest';
  1         1  
33             use Math::PlanePath::Base::Digits
34 1     1   486 'round_down_pow';
  1         2  
  1         47  
35 1     1   524 use Math::PlanePath::KochCurve;
  1         2  
  1         39  
36             *_divrem_mutate = \&Math::PlanePath::_divrem_mutate;
37              
38             # uncomment this to run the ### lines
39             #use Devel::Comments;
40              
41              
42 1     1   6 use constant class_y_negative => 0;
  1         2  
  1         47  
43 1     1   5 use constant n_frac_discontinuity => .5;
  1         1  
  1         37  
44 1     1   4 use constant x_negative_at_n => 1;
  1         2  
  1         34  
45 1     1   4 use constant sumabsxy_minimum => 1; # minimum X=1,Y=0
  1         2  
  1         32  
46 1     1   4 use constant absdiffxy_minimum => 1; # X=Y never occurs
  1         3  
  1         33  
47 1     1   5 use constant rsquared_minimum => 1; # minimum X=1,Y=0
  1         1  
  1         32  
48              
49 1     1   5 use constant dx_maximum => 2;
  1         1  
  1         33  
50 1     1   5 use constant dy_minimum => -1;
  1         1  
  1         55  
51 1     1   5 use constant dy_maximum => 1;
  1         2  
  1         35  
52 1     1   5 use constant absdx_minimum => 1; # never vertical
  1         2  
  1         41  
53 1     1   5 use constant dsumxy_maximum => 2; # diagonal NE
  1         2  
  1         34  
54 1     1   4 use constant ddiffxy_maximum => 2; # diagonal NW
  1         2  
  1         49  
55 1     1   5 use constant dir_maximum_dxdy => (1,-1); # South-East
  1         1  
  1         45  
56 1     1   6 use constant turn_any_straight => 0; # never straight
  1         1  
  1         767  
57              
58              
59             #------------------------------------------------------------------------------
60              
61             # N=1 to 3 3 of, level=0
62             # N=4 to 12 9 of, level=1
63             # N=13 to 45 33 of, level=2
64             #
65             # N=0.5 to 3.49 diff=3
66             # N=3.39 to 12.49 diff=9
67             # N=12.5 to 45.5 diff=33
68             #
69             # each length = 2*4^level + 1
70             #
71             # Nstart = 1 + 2*4^0 + 1 + 2*4^1 + 1 + ... + 2*4^(level-1) + 1
72             # = 1 + level + 2*[ 4^0 + 4^1 + ... + 4^(level-1) ]
73             # = level+1 + 2*[ (4^level - 1)/3 ]
74             # = level+1 + (2*4^level - 2)/3
75             # = level + (2*4^level - 2 + 3)/3
76             # = level + (2*4^level + 1)/3
77             #
78             # 3*n = 2*4^level + 1
79             # 3*n-1 = 2*4^level
80             # (3*n-1)/2 = 4^level
81             #
82             # Nbase = 0.5 + 2*4^0 + 1 + 2*4^1 + 1 + ... + 2*4^(level-1) + 1
83             # = level + (2*4^level + 1)/3 - 1/2
84             # = level + 2/3*4^level + 1/3 - 1/2
85             # = level + 2/3*4^level - 1/6
86             # = level + 4/6*4^level - 1/6
87             # = level + (4*4^level - 1)/6
88             # = level + (4^(level+1) - 1)/6
89             #
90             # 6*N = 4^(level+1) - 1
91             # 6*N + 1 = 4^(level+1)
92             # level+1 = log4(6*N + 1)
93             # level = log4(6*N + 1) - 1
94             #
95             ### loop 1: (2*4**1 + 1)/3
96             ### loop 2: (2*4**2 + 1)/3
97             ### loop 3: (2*4**3 + 1)/3
98              
99             # sub _n_to_level {
100             # my ($n) = @_;
101             # my ($side, $level) = round_down_pow(6*$n + 1, 4);
102             # my $base = $level + (2*$side + 1)/3 - .5;
103             # ### $level
104             # ### $base
105             # if ($base > $n) {
106             # $level--;
107             # $side /= 4;
108             # $base = $level + (2*$side + 1)/3 - .5;
109             # ### $level
110             # ### $base
111             # }
112             # return ($level, $base, $side + .5);
113             # }
114             # sub _level_to_base {
115             # my ($level) = @_;
116             # return $level + (2*$side + 1)/3 - .5;
117             # }
118              
119             sub _n_to_side_level_base {
120 381     381   479 my ($n) = @_;
121 381         723 my ($side, $level) = round_down_pow((3*$n-1)/2, 4);
122 381         585 my $base = $level + (2*$side + 1)/3;
123             ### $level
124             ### $base
125 381 100       613 if (2*$n+1 < 2*$base) {
126 11         16 $level--;
127 11         14 $side /= 4;
128 11         15 $base = $level + (2*$side + 1)/3;
129             ### $level
130             ### $base
131             }
132 381         603 return ($side, $level, $base);
133             }
134              
135             sub n_to_xy {
136 381     381 1 3325 my ($self, $n) = @_;
137             ### KochPeaks n_to_xy(): $n
138              
139             # $n<0.5 no good for Math::BigInt circa Perl 5.12, compare in integers
140 381 50       569 return if 2*$n < 1;
141              
142 381 50       609 if (is_infinite($n)) { return ($n,$n); }
  0         0  
143              
144 381         587 my ($side, $level, $base) = _n_to_side_level_base($n);
145              
146 381         500 my $rem = $n - $base;
147 381         379 my $frac;
148 381 100       620 if ($rem < 0) {
    100          
149             ### neg frac
150 2         3 $frac = $rem;
151 2         2 $rem = 0;
152             } elsif ($rem > 2*$side) {
153             ### excess frac
154 1         1 $frac = $rem - 2*$side;
155 1         2 $rem -= $frac;
156             } else {
157             ### no frac
158 378         423 $frac = 0;
159             }
160             ### $frac
161             ### $rem
162             ### $n
163              
164             ### next base would be: ($level+1) + (2*4**($level+1) + 1)/3
165             ### assert: $n-$frac >= $base
166             ### assert: $n-$frac < ($level+1) + (2*4**($level+1) + 1)/3
167              
168             ### assert: $rem>=0
169             ### assert: $rem < 2 * 4 ** $level + 1
170             ### assert: $rem <= 2*$side+1
171              
172 381         459 my $pos = 3**$level;
173 381 100       500 if ($rem < $side) {
174 361         638 my ($x, $y) = Math::PlanePath::KochCurve->n_to_xy($rem);
175             ### left side: $rem
176             ### flat: "$x,$y"
177 361         467 $x += 2*$frac;
178 361         847 return (($x-3*$y)/2 - $pos, # rotate +60
179             ($x+$y)/2);
180             } else {
181 20         44 my ($x, $y) = Math::PlanePath::KochCurve->n_to_xy($rem-$side);
182             ### right side: $rem-$side
183             ### flat: "$x,$y"
184 20         30 $x += 2*$frac;
185 20         64 return (($x+3*$y)/2, # rotate -60
186             ($y-$x)/2 + $pos);
187             }
188             }
189              
190             sub xy_to_n {
191 5631     5631 1 34232 my ($self, $x, $y) = @_;
192             ### KochPeaks xy_to_n(): "$x, $y"
193              
194 5631         8396 $x = round_nearest ($x);
195 5631         8775 $y = round_nearest ($y);
196              
197 5631 100 66     13857 if ($y < 0 || ! (($x ^ $y) & 1)) {
198             ### neg y or parity...
199 265         373 return undef;
200             }
201 5366         9121 my ($len,$level) = round_down_pow ($y+abs($x), 3);
202             ### $level
203             ### $len
204 5366 50       8574 if (is_infinite($level)) {
205 0         0 return $level;
206             }
207              
208 5366         6902 my $n;
209 5366 100       7057 if ($x < 0) {
210 259         273 $x += $len;
211 259         424 ($x,$y) = (($x+3*$y)/2, # rotate -60
212             ($y-$x)/2);
213 259         293 $n = 0;
214             ### left rotate -60 to: "x=$x,y=$y n=$n"
215             } else {
216 5107         5372 $y -= $len;
217 5107         8581 ($x,$y) = (($x-3*$y)/2, # rotate +60
218             ($x+$y)/2);
219 5107         6027 $n = 1;
220             ### right rotate +60 to: "x=$x,y=$y n=$n"
221             }
222              
223 5366         7649 foreach (1 .. $level) {
224 19311         20039 $n *= 4;
225             ### at: "level=$level len=$len x=$x,y=$y n=$n"
226 19311 100       24098 if ($x < $len) {
227 11429         11834 $len /= 3;
228 11429         12738 my $rel = 2*$len;
229 11429 100       16187 if ($x < $rel) {
230             ### digit 0
231             } else {
232             ### digit 1 sub: "$rel to x=".($x-$rel)
233 1967         2030 $x -= $rel;
234 1967         2972 ($x,$y) = (($x+3*$y)/2, # rotate -60
235             ($y-$x)/2);
236 1967         2682 $n += 1;
237             }
238             } else {
239 7882         8085 $len /= 3;
240 7882         8638 $x -= 4*$len;
241 7882 100       9743 if ($x < $y) { # before diagonal
242             ### digit 2...
243 3548         5770 ($x,$y) = (($x-3*$y)/2 + 2*$len, # rotate +60
244             ($x+$y)/2);
245 3548         4515 $n += 2;
246             } else {
247             #### digit 3...
248 4334         5411 $n += 3;
249             }
250             }
251             }
252             ### end at: "x=$x,y=$y n=$n"
253 5366 100       7551 if ($x) {
254             ### endmost point
255 4814         5141 $n += 1;
256 4814         5033 $x -= 2;
257             }
258 5366 100 100     8699 if ($x != 0 || $y != 0) {
259 4987         8323 return undef;
260             }
261 379         810 return $n + $level + (2*4**$level + 1)/3 + ($x == 2);
262             }
263              
264              
265             # level extends to x= +/- 3^level
266             # y= 0 to 3^level
267             #
268             # diagonal X=Y or Y=-X is lowest in a level, so round down abs(X)+Y to pow 3
269             #
270             # end of level is 1 before base of level+1
271             # basenext = (level+1) + (2*4^(level+1) + 1)/3
272             # basenext-1 = level + (2*4^(level+1) + 1)/3
273             # = level + (8*4^level + 1)/3
274              
275             # not exact
276             sub rect_to_n_range {
277 17     17 1 1415 my ($self, $x1,$y1, $x2,$y2) = @_;
278             ### KochPeaks rect_to_n_range(): "$x1,$y1 $x2,$y2"
279              
280 17         30 $x1 = round_nearest ($x1);
281 17         30 $y1 = round_nearest ($y1);
282 17         35 $x2 = round_nearest ($x2);
283 17         25 $y2 = round_nearest ($y2);
284             ### rounded: "$x1,$y1 $x2,$y2"
285              
286 17 50 66     36 if ($y1 < 0 && $y2 < 0) {
287 0         0 return (1,0);
288             }
289              
290             # can't make use of the len=3**$level returned by round_down_pow()
291 17         34 my ($len, $level) = round_down_pow (max(abs($x1),abs($x2))
292             + max($y1, $y2),
293             3);
294             ### $level
295 17         45 return (1, $level + (8 * 4**$level + 1)/3);
296             }
297              
298             # peak Y is at N = Nstart + (count-1)/2
299             # = level + (2*4^level + 1)/3 + (2*4^level + 1 - 1)/2
300             # = level + (2*4^level + 1)/3 + (2*4^level)/2
301             # = level + (2*4^level + 1)/3 + 4^level
302             # = level + (2*4^level + 1 + 3*4^level)/3
303             # = level + (5*4^level + 1)/3
304              
305             #------------------------------------------------------------------------------
306              
307             sub level_to_n_range {
308 10     10 1 585 my ($self, $level) = @_;
309 10         17 my $pow = 4**$level;
310 10         30 return ((2*$pow + 1)/3 + $level,
311             (8*$pow + 1)/3 + $level);
312             }
313             sub n_to_level {
314 0     0 1   my ($self, $n) = @_;
315 0 0         if ($n < 1) { return undef; }
  0            
316 0 0         if (is_infinite($n)) { return $n; }
  0            
317 0           $n = round_nearest($n);
318 0           my ($side, $level, $base) = _n_to_side_level_base($n);
319 0           return $level;
320             }
321              
322             #------------------------------------------------------------------------------
323             1;
324             __END__