File Coverage

blib/lib/Math/PlanePath/QuadricIslands.pm
Criterion Covered Total %
statement 75 121 61.9
branch 6 26 23.0
condition 0 9 0.0
subroutine 22 26 84.6
pod 6 6 100.0
total 109 188 57.9


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             # math-image --path=QuadricIslands --lines --scale=10
20             # math-image --path=QuadricIslands --all --output=numbers_dash --size=132x50
21              
22              
23             package Math::PlanePath::QuadricIslands;
24 1     1   9022 use 5.004;
  1         9  
25 1     1   6 use strict;
  1         2  
  1         48  
26             #use List::Util 'max';
27             *max = \&Math::PlanePath::_max;
28              
29 1     1   7 use vars '$VERSION', '@ISA';
  1         2  
  1         64  
30             $VERSION = 128;
31 1     1   641 use Math::PlanePath;
  1         2  
  1         44  
32             @ISA = ('Math::PlanePath');
33              
34             use Math::PlanePath::Base::Generic
35 1         47 'is_infinite',
36             'round_nearest',
37 1     1   7 'floor';
  1         2  
38             use Math::PlanePath::Base::Digits
39 1     1   458 'round_down_pow';
  1         2  
  1         57  
40              
41 1     1   477 use Math::PlanePath::QuadricCurve;
  1         3  
  1         31  
42              
43             # uncomment this to run the ### lines
44             #use Smart::Comments;
45              
46              
47 1     1   7 use constant n_frac_discontinuity => 0;
  1         2  
  1         49  
48 1     1   6 use constant x_negative_at_n => 1;
  1         2  
  1         38  
49 1     1   5 use constant y_negative_at_n => 1;
  1         2  
  1         37  
50 1     1   5 use constant sumabsxy_minimum => 1; # minimum X=1/2,Y=1/2
  1         2  
  1         48  
51 1     1   7 use constant rsquared_minimum => 0.5; # minimum X=1/2,Y=1/2
  1         2  
  1         41  
52              
53 1     1   5 use constant dx_maximum => 1;
  1         2  
  1         35  
54 1     1   5 use constant dy_maximum => 1;
  1         2  
  1         36  
55 1     1   5 use constant dsumxy_maximum => 1;
  1         1  
  1         65  
56 1     1   7 use constant ddiffxy_minimum => -1; # dDiffXY=+1 or -1
  1         1  
  1         52  
57 1     1   7 use constant ddiffxy_maximum => 1;
  1         1  
  1         43  
58 1     1   6 use constant dir_maximum_dxdy => (0,-1); # South
  1         9  
  1         58  
59              
60             # N=1,2,3,4 gcd(1/2,1/2) = 1/2
61 1     1   6 use constant gcdxy_minimum => 1/2;
  1         2  
  1         845  
62              
63             #------------------------------------------------------------------------------
64              
65             # N=1 to 4 4 of, level=0
66             # N=5 to 36 12 of, level=1
67             # N=37 to .. 48 of, level=3
68             #
69             # each loop = 4*8^level
70             #
71             # n_base = 1 + 4*8^0 + 4*8^1 + ... + 4*8^(level-1)
72             # = 1 + 4*[ 8^0 + 8^1 + ... + 8^(level-1) ]
73             # = 1 + 4*[ (8^level - 1)/7 ]
74             # = 1 + 4*(8^level - 1)/7
75             # = (4*8^level - 4 + 7)/7
76             # = (4*8^level + 3)/7
77             #
78             # n >= (4*8^level + 3)/7
79             # 7*n = 4*8^level + 3
80             # (7*n - 3)/4 = 8^level
81             #
82             # nbase(k+1)-nbase(k)
83             # = (4*8^(k+1)+3 - (4*8^k+3)) / 7
84             # = (4*8*8^k - 4*8^k) / 7
85             # = (4*8-4) * 8^k / 7
86             # = 28 * 8^k / 7
87             # = 4 * 8^k
88             #
89             # nbase(0) = (4*8^0 + 3)/7 = (4+3)/7 = 1
90             # nbase(1) = (4*8^1 + 3)/7 = (4*8+3)/7 = (32+3)/7 = 35/7 = 5
91             # nbase(2) = (4*8^2 + 3)/7 = (4*64+3)/7 = (256+3)/7 = 259/7 = 37
92             #
93             ### loop 1: 4* 8**1
94             ### loop 2: 4* 8**2
95             ### loop 3: 4* 8**3
96              
97             # sub _level_to_base {
98             # my ($level) = @_;
99             # return (4*8**$level + 3) / 7;
100             # }
101             # ### level_to_base(1): _level_to_base(1)
102             # ### level_to_base(2): _level_to_base(2)
103             # ### level_to_base(3): _level_to_base(3)
104              
105             # base = (4 * 8**$level + 3)/7
106             # = (4 * 8**($level+1) / 8 + 3)/7
107             # = (8**($level+1) / 2 + 3)/7
108             sub _n_to_base_and_level {
109 14     14   27 my ($n) = @_;
110 14         36 my ($base,$level) = round_down_pow ((7*$n - 3)*2, 8);
111 14         38 return (($base/2 + 3)/7,
112             $level - 1);
113             }
114              
115             sub n_to_xy {
116 14     14 1 65 my ($self, $n) = @_;
117             ### QuadricIslands n_to_xy(): "$n"
118 14 50       32 if ($n < 1) { return; }
  0         0  
119 14 50       34 if (is_infinite($n)) { return ($n,$n); }
  0         0  
120              
121 14         28 my ($base, $level) = _n_to_base_and_level($n);
122             ### $level
123             ### $base
124             ### level: "$level"
125             ### next base would be: (4 * 8**($level+1) + 3)/7
126              
127 14         27 my $rem = $n - $base;
128             ### $rem
129              
130             ### assert: $n >= $base
131             ### assert: $n < 8**($level+1)
132             ### assert: $rem>=0
133             ### assert: $rem < 4 * 8 ** $level
134              
135 14         19 my $sidelen = 8**$level;
136 14         25 my $side = int($rem / $sidelen);
137             ### $sidelen
138             ### $side
139             ### $rem
140 14         19 $rem -= $side*$sidelen;
141             ### assert: $side >= 0 && $side < 4
142 14         32 my ($x, $y) = Math::PlanePath::QuadricCurve::n_to_xy ($self, $rem);
143              
144 14         27 my $pos = 4**$level / 2;
145             ### side calc: "$x,$y for pos $pos"
146             ### $x
147             ### $y
148              
149 14 100       35 if ($side < 1) {
    50          
    50          
150             ### horizontal rightwards
151 7         22 return ($x - $pos,
152             $y - $pos);
153             } elsif ($side < 2) {
154             ### right vertical upwards
155 0         0 return (-$y + $pos, # rotate +90, offset
156             $x - $pos);
157             } elsif ($side < 3) {
158             ### horizontal leftwards
159 0         0 return (-$x + $pos, # rotate 180, offset
160             -$y + $pos);
161             } else {
162             ### left vertical downwards
163 7         24 return ($y - $pos, # rotate -90, offset
164             -$x + $pos);
165             }
166             }
167              
168             # +-------+-------+-------+
169             # |31 | 24 0,1| 23|
170             # | | | |
171             # | +-------+-------+ |
172             # | |4 | |3 | | |
173             # | | | | | | |
174             # +---|--- ---|--- ---|---+ Y=0.5
175             # |32 | | | | | 16|
176             # | | | | | | |
177             # | +=======+=======+ | Y=0
178             # | |1 | |2 | | |
179             # | | | | | | |
180             # +---|--- ---|--- ---|---+ Y=-0.5
181             # | | | | | | |
182             # | | | | | | |
183             # | +-------+-------+ | Y=-1
184             # | | | |
185             # |7 |8 | 15|
186             # +-------+-------+-------+
187             #
188             # -2 <= 2*x < 2, round to -2,-1,0,1
189             # then 4*yround -8,-4,0,4
190             # total -10 to 5 inclusive
191              
192             my @inner_n_list = ([1,7], [1,8], [2,8], [2,15], # Y=-1
193             [1,32], [1], [2], [2,16], # Y=-0.5
194             [4,32], [4], [3], [3,16], # Y=0
195             [4,31],[4,24],[3,24],[3,23]); # Y=0.5
196              
197             sub xy_to_n {
198 0     0 1 0 return scalar((shift->xy_to_n_list(@_))[0]);
199             }
200             sub xy_to_n_list {
201 0     0 1 0 my ($self, $x, $y) = @_;
202             ### QuadricIslands xy_to_n(): "$x, $y"
203              
204 0 0 0     0 if ($x >= -1 && $x < 1
      0        
      0        
205             && $y >= -1 && $y < 1) {
206              
207             ### round 2x: floor(2*$x)
208             ### round 2y: floor(2*$y)
209             ### index: floor(2*$x) + 4*floor(2*$y) + 10
210             ### assert: floor(2*$x) + 4*floor(2*$y) + 10 >= 0
211             ### assert: floor(2*$x) + 4*floor(2*$y) + 10 <= 15
212              
213 0         0 return @{$inner_n_list[ floor(2*$x)
  0         0  
214             + 4*floor(2*$y)
215             + 10 ]};
216             }
217              
218 0         0 $x = round_nearest($x);
219 0         0 $y = round_nearest($y);
220              
221 0         0 my $high;
222 0 0       0 if ($x >= $y + ($y>0)) {
223             # +($y>0) to exclude the downward bump of the top side
224             ### below leading diagonal ...
225 0 0       0 if ($x < -$y) {
226             ### bottom quarter ...
227 0         0 $high = 0;
228             } else {
229             ### right quarter ...
230 0         0 $high = 1;
231 0         0 ($x,$y) = ($y, -$x); # rotate -90
232             }
233             } else {
234             ### above leading diagonal
235 0 0       0 if ($y > -$x) {
236             ### top quarter ...
237 0         0 $high = 2;
238 0         0 $x = -$x; # rotate 180
239 0         0 $y = -$y;
240             } else {
241             ### right quarter ...
242 0         0 $high = 3;
243 0         0 ($x,$y) = (-$y, $x); # rotate +90
244             }
245             }
246             ### rotate to: "$x,$y high=$high"
247              
248             # ymax = (10*4^(l-1)-1)/3
249             # ymax < (10*4^(l-1)-1)/3+1
250             # (10*4^(l-1)-1)/3+1 > ymax
251             # (10*4^(l-1)-1)/3 > ymax-1
252             # 10*4^(l-1)-1 > 3*(ymax-1)
253             # 10*4^(l-1) > 3*(ymax-1)+1
254             # 10*4^(l-1) > 3*(ymax-1)+1
255             # 10*4^(l-1) > 3*ymax-3+1
256             # 10*4^(l-1) > 3*ymax-2
257             # 4^(l-1) > (3*ymax-2)/10
258             #
259             # (2*4^(l-1) + 1)/3 = ymin
260             # 2*4^(l-1) + 1 = 3*y
261             # 2*4^(l-1) = 3*y-1
262             # 4^(l-1) = (3*y-1)/2
263             #
264             # ypos = 4^l/2 = 2*4^(l-1)
265              
266              
267             # z = -2*y+x
268             # (2*4**($level-1) + 1)/3 = z
269             # 2*4**($level-1) + 1 = 3*z
270             # 2*4**($level-1) = 3*z - 1
271             # 4**($level-1) = (3*z - 1)/2
272             # = (3*(-2y+x)-1)/2
273             # = (-6y+3x - 1)/2
274             # = -3*y + (3x-1)/2
275              
276             # 2*4**($level-1) = -2*y-x
277             # 4**($level-1) = -y-x/2
278             # 4**$level = -4y-2x
279             #
280             # line slope y/x = 1/2 as an index
281 0         0 my $z = -$y-$x/2;
282 0         0 my ($len,$level) = round_down_pow ($z, 4);
283             ### $z
284             ### amin: 2*4**($level-1)
285             ### $level
286             ### $len
287 0 0       0 if (is_infinite($level)) {
288 0         0 return $level;
289             }
290              
291 0         0 $len *= 2;
292 0         0 $x += $len;
293 0         0 $y += $len;
294             ### shift to: "$x,$y"
295 0         0 my $n = Math::PlanePath::QuadricCurve::xy_to_n($self, $x, $y);
296              
297             # Nmin = (4*8^l+3)/7
298             # Nmin+high = (4*8^l+3)/7 + h*8^l
299             # = (4*8^l + 3 + 7h*8^l)/7 +
300             # = ((4+7h)*8^l + 3)/7
301             #
302             ### plain curve on: ($x+2*$len).",".($y+2*$len)." give n=".(defined $n && $n)
303             ### $high
304             ### high: (8**$level)*$high
305             ### base: (4 * 8**($level+1) + 3)/7
306             ### base with high: ((4+7*$high) * 8**($level+1) + 3)/7
307              
308 0 0       0 if (defined $n) {
309 0         0 return ((4+7*$high) * 8**($level+1) + 3)/7 + $n;
310             } else {
311 0         0 return;
312             }
313             }
314              
315             # level width extends
316             # side = 4^level
317             # ypos = 4^l / 2
318             # width = 1 + 4 + ... + 4^(l-1)
319             # = (4^l - 1)/3
320             # ymin = ypos(l) - 4^(l-1) - width(l-1)
321             # = 4^l / 2 - 4^(l-1) - (4^(l-1) - 1)/3
322             # = 4^(l-1) * (2 - 1 - 1/3) + 1/3
323             # = (2*4^(l-1) + 1) / 3
324             #
325             # (2*4^(l-1) + 1) / 3 = y
326             # 2*4^(l-1) + 1 = 3*y
327             # 2*4^(l-1) = 3*y-1
328             # 4^(l-1) = (3*y-1)/2
329             #
330             # ENHANCE-ME: slope Y=X/2+1 or thereabouts for sides
331             #
332             # not exact
333             sub rect_to_n_range {
334 0     0 1 0 my ($self, $x1,$y1, $x2,$y2) = @_;
335             ### QuadricIslands rect_to_n_range(): "$x1,$y1 $x2,$y2"
336              
337             # $x1 = round_nearest ($x1);
338             # $y1 = round_nearest ($y1);
339             # $x2 = round_nearest ($x2);
340             # $y2 = round_nearest ($y2);
341              
342 0         0 my $m = max(abs($x1), abs($x2),
343             abs($y1), abs($y2));
344              
345 0         0 my ($len,$level) = round_down_pow (6*$m-2, 4);
346             ### $len
347             ### $level
348 0         0 return (1,
349             (32*8**$level - 4)/7);
350             }
351              
352             # ymax = ypos(l) + 4^(l-1) + width(l-1)
353             # = 4^l / 2 + 4^(l-1) + (4^(l-1) - 1)/3
354             # = 4^(l-1) * (4/2 + 1 + 1/3) - 1/3
355             # = 4^(l-1) * (2 + 1 + 1/3) - 1/3
356             # = 4^(l-1) * 10/3 - 1/3
357             # = (10*4^(l-1) - 1) / 3
358             #
359             # (10*4^(l-1) - 1) / 3 = y
360             # 10*4^(l-1) - 1 = 3*y
361             # 10*4^(l-1) = 3*y+1
362             # 4^(l-1) = (3*y+1)/10
363             #
364             # based on max ??? ...
365             #
366             # my ($power,$level) = round_down_pow ((3*$m+1-3)/10, 4);
367             # ### $power
368             # ### $level
369             # return (1,
370             # (4*8**($level+3) + 3)/7 - 1);
371              
372             #------------------------------------------------------------------------------
373             # Nstart(k) = (4*8^k + 3)/7
374             # Nend(k) = Nstart(k+1) - 1
375             # = (4*8*8^k + 3)/7 - 1
376             # = (4*8*8^k + 8*3 - 8*3 + 3)/7 - 1
377             # = (4*8*8^k + 8*3)/7 + (-8*3 + 3)/7 - 1
378             # = 8*Nstart(k) + (-8*3 + 3)/7 - 1
379             # = 8*Nstart(k) - 4
380              
381             sub level_to_n_range {
382 10     10 1 886 my ($self, $level) = @_;
383 10         29 my $n_lo = (4 * 8**$level + 3)/7;
384 10         30 return ($n_lo, 8*$n_lo - 4);
385             }
386             sub n_to_level {
387 0     0 1   my ($self, $n) = @_;
388 0 0         if ($n < 1) { return undef; }
  0            
389 0 0         if (is_infinite($n)) { return $n; }
  0            
390 0           $n = round_nearest($n);
391 0           my ($base,$level) = _n_to_base_and_level($n);
392 0           return $level;
393             }
394              
395             #------------------------------------------------------------------------------
396             1;
397             __END__