File Coverage

blib/lib/Math/PlanePath/PeanoCurve.pm
Criterion Covered Total %
statement 174 198 87.8
branch 38 64 59.3
condition 22 28 78.5
subroutine 19 26 73.0
pod 7 7 100.0
total 260 323 80.5


line stmt bran cond sub pod time code
1             # Copyright 2010, 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             # cf
20             #
21             # http://www.cut-the-knot.org/Curriculum/Geometry/PeanoComplete.shtml
22             # applet, directions in 9 sub-parts
23             #
24             # math-image --path=PeanoCurve,radix=5 --all --output=numbers
25             # math-image --path=PeanoCurve,radix=5 --lines
26             #
27             # -----------
28             # Peano:
29             # T = 0.a1 a2 a3 a4 ...
30             # x y x y
31             #
32             # X = 0.b1 b2 ...
33             # a1 a3.k(a2)
34             #
35             # Y = 0.c1 c2 ...
36             # a2.k(a1) a4.k(a1,a3)
37             #
38             # b1=a1
39             # c1 = a2 comp(a1)
40             # b2 = a3 comp(a2)
41             # c2 = a4 comp(a1+a3)
42             #
43             # bn = a[2n-1] comp a2+a4+...+a[2n-2]
44             # cn = a[2n] comp a1+a3+...+a[2n-1]
45             #
46             # Brouwer(?) no continuous one-to-one between R and RxR, so line and plane
47             # are distinguished.
48             #
49              
50              
51             package Math::PlanePath::PeanoCurve;
52 5     5   4891 use 5.004;
  5         25  
53 5     5   28 use strict;
  5         10  
  5         230  
54             #use List::Util 'max';
55             *max = \&Math::PlanePath::_max;
56              
57 5     5   31 use vars '$VERSION', '@ISA';
  5         19  
  5         370  
58             $VERSION = 128;
59 5     5   1208 use Math::PlanePath;
  5         11  
  5         205  
60             @ISA = ('Math::PlanePath');
61              
62             use Math::PlanePath::Base::Generic
63 5         272 'is_infinite',
64 5     5   30 'round_nearest';
  5         43  
65             use Math::PlanePath::Base::Digits
66 5         312 'round_down_pow',
67             'digit_split_lowtohigh',
68 5     5   1314 'digit_join_lowtohigh';
  5         12  
69 5     5   1783 use Math::PlanePath::Base::NSEW;
  5         12  
  5         140  
70              
71             # uncomment this to run the ### lines
72             # use Smart::Comments;
73              
74              
75 5     5   28 use constant n_start => 0;
  5         11  
  5         239  
76 5     5   27 use constant class_x_negative => 0;
  5         10  
  5         224  
77 5     5   26 use constant class_y_negative => 0;
  5         10  
  5         354  
78             *xy_is_visited = \&Math::PlanePath::Base::Generic::xy_is_visited_quad1;
79              
80 5         8440 use constant parameter_info_array =>
81             [ { name => 'radix',
82             display => 'Radix',
83             share_key => 'radix_3',
84             type => 'integer',
85             minimum => 2,
86             default => 3,
87             width => 3,
88 5     5   29 } ];
  5         10  
89              
90             # shared by WunderlichSerpentine
91             sub dx_minimum {
92 0     0 1 0 my ($self) = @_;
93 0 0       0 return ($self->{'radix'} % 2
94             ? -1 # odd
95             : undef); # even, unlimited
96             }
97             sub dx_maximum {
98 0     0 1 0 my ($self) = @_;
99 0 0       0 return ($self->{'radix'} % 2
100             ? 1 # odd
101             : undef); # even, unlimited
102             }
103              
104             # shared by WunderlichSerpentine
105             sub _UNDOCUMENTED__dxdy_list {
106 0     0   0 my ($self) = @_;
107 0 0       0 return ($self->{'radix'} % 2
108             ? Math::PlanePath::Base::NSEW->_UNDOCUMENTED__dxdy_list
109             : ()); # even, unlimited
110             }
111             # *--- b^2-1 -- b^2 ---- b^2+b-1 = (b+1)b-1
112             # | |
113             # *-------
114             # |
115             # 0 ----- b
116             #
117             sub _UNDOCUMENTED__dxdy_list_at_n {
118 0     0   0 my ($self) = @_;
119 0         0 return ($self->{'radix'} + 1) * $self->{'radix'} - 1;
120             }
121              
122             # shared by WunderlichSerpentine
123             *dy_minimum = \&dx_minimum;
124             *dy_maximum = \&dx_maximum;
125              
126             *dsumxy_minimum = \&dx_minimum;
127             *dsumxy_maximum = \&dx_maximum;
128              
129             *ddiffxy_minimum = \&dx_minimum;
130             *ddiffxy_maximum = \&dx_maximum;
131              
132             sub dir_maximum_dxdy {
133 0     0 1 0 my ($self) = @_;
134 0 0       0 return ($self->{'radix'} % 2
135             ? (0,-1) # odd, South
136             : (0,0)); # even, supremum
137             }
138              
139             sub _UNDOCUMENTED__turn_any_left_at_n {
140 0     0   0 my ($self) = @_;
141 0         0 return $self->{'radix'} - 1;
142             }
143             sub _UNDOCUMENTED__turn_any_right_at_n {
144 0     0   0 my ($self) = @_;
145             return ($self->{'radix'} == 2 ? 5
146 0 0       0 : 2*$self->{'radix'} - 1);
147             }
148              
149              
150             #------------------------------------------------------------------------------
151              
152             sub new {
153 11     11 1 1817 my $self = shift->SUPER::new(@_);
154              
155 11 100 66     80 if (! $self->{'radix'} || $self->{'radix'} < 2) {
156 7         17 $self->{'radix'} = 3;
157             }
158 11         26 return $self;
159             }
160              
161             sub _n_to_xykk {
162 60585     60585   96088 my ($self, $n) = @_;
163 60585         93536 my $radix = $self->{'radix'};
164 60585         84895 my $radix_minus_1 = $radix - 1;
165              
166 60585         117943 my @ndigits = digit_split_lowtohigh($n,$radix);
167 60585 100       121394 if (scalar(@ndigits) & 1) {
168 52494         74805 push @ndigits, 0; # so even number of entries
169             }
170             ### @ndigits
171              
172 60585         88743 my $xk = 0;
173 60585         77566 my $yk = 0;
174 60585         88268 my @ydigits;
175             my @xdigits;
176              
177 60585         123616 for (my $i = $#ndigits >> 1; @ndigits; $i--) { # high to low
178             ### $i
179             {
180 172473         227655 my $ndigit = pop @ndigits; # high to low
181 172473         241710 $xk ^= $ndigit;
182 172473 100       493929 $ydigits[$i] = ($yk & 1 ? $radix_minus_1-$ndigit : $ndigit);
183             }
184             {
185 172473         440252 my $ndigit = pop @ndigits;
  172473         567520  
  172473         222167  
186 172473         226061 $yk ^= $ndigit;
187 172473 100       638821 $xdigits[$i] = ($xk & 1 ? $radix_minus_1-$ndigit : $ndigit);
188             }
189             }
190 60585         89540 my $zero = $n*0;
191 60585         97157 return ((map {digit_join_lowtohigh($_, $radix, $zero)} \@xdigits, \@ydigits),
  121170         227851  
192             $xk,$yk);
193             }
194              
195              
196             sub n_to_xy {
197 20183     20183 1 116602 my ($self, $n) = @_;
198             ### PeanoCurve n_to_xy(): $n
199 20183 50       38257 if ($n < 0) { # negative
200 0         0 return;
201             }
202 20183 50       40594 if (is_infinite($n)) {
203 0         0 return ($n,$n);
204             }
205              
206             {
207             # ENHANCE-ME: for odd radix the ends join and the direction can be had
208             # without a full N+1 calculation
209 20183         38208 my $int = int($n);
  20183         28111  
210             ### $int
211             ### $n
212 20183 100       35884 if ($n != $int) {
213 1         360 my ($x1,$y1) = $self->n_to_xy($int);
214 1         7 my ($x2,$y2) = $self->n_to_xy($int+1);
215 1         12 my $frac = $n - $int; # inherit possible BigFloat
216 1         600 my $dx = $x2-$x1;
217 1         319 my $dy = $y2-$y1;
218 1         253 return ($frac*$dx + $x1, $frac*$dy + $y1);
219             }
220 20182         30976 $n = $int; # BigFloat int() gives BigInt, use that
221             }
222              
223 20182         37210 my ($x,$y) = _n_to_xykk($self,$n);
224 20182         53702 return ($x,$y);
225             }
226              
227             sub _xykk_to_n {
228 853     853   1547 my ($self, $x,$y, $offset_xk,$offset_yk) = @_;
229             ### PeanoCurve _xykk_to_n(): "$x, $y offset $offset_xk,$offset_yk"
230              
231 853 100 100     2698 if (($offset_xk && ($x-=$offset_xk) < 0)
      100        
      100        
232             || ($offset_yk && ($y-=$offset_yk) < 0)) {
233 11         36 return; # offset goes negative
234             }
235              
236 842         1417 my $radix = $self->{'radix'};
237 842         1730 my @x = digit_split_lowtohigh ($x, $radix);
238 842         1607 my @y = digit_split_lowtohigh ($y, $radix);
239              
240 842         1280 my $radix_minus_1 = $radix - 1;
241 842         1158 my $xk = 0;
242 842         1136 my $yk = 0;
243              
244 842         1063 my @n; # stored low to high, generated from high to low
245 842         2627 my $i_high = max($#x,$#y);
246 842         1319 my $npos = 2*$i_high+1;
247              
248 842         1809 foreach my $i (reverse 0 .. $i_high) { # high to low
249             {
250 4159   100     7984 my $digit = $y[$i] || 0;
251 4159 100       7041 if ($yk & 1) {
252 1670         2096 $digit = $radix_minus_1 - $digit; # reverse digit
253             }
254 4159         6003 $n[$npos--] = $digit;
255 4159         5607 $xk ^= $digit;
256             }
257             {
258 4159   100     5014 my $digit = $x[$i] || 0;
  4159         5133  
  4159         7840  
259 4159 100       7124 if ($xk & 1) {
260 2084         2727 $digit = $radix_minus_1 - $digit; # reverse digit
261             }
262 4159         6061 $n[$npos--] = $digit;
263 4159         5843 $yk ^= $digit;
264             }
265             }
266             ### final n: @n
267             ### final xkyk: ($xk&1).' '.($yk&1)
268 842 100 100     4584 return ((! defined $offset_xk || ($xk&1) == $offset_xk)
269             && (! defined $offset_yk || ($yk&1) == $offset_yk)
270             ? (digit_join_lowtohigh (\@n, $radix,
271             $x*0*$y)) # inherit bignum 0
272             : ());
273             }
274              
275             sub xy_to_n {
276 731     731 1 1482 my ($self, $x, $y) = @_;
277             ### PeanoCurve xy_to_n(): "$x, $y"
278              
279 731         1795 $x = round_nearest ($x);
280 731         1445 $y = round_nearest ($y);
281              
282 731 50 33     2713 if ($x < 0 || $y < 0) { return undef; }
  0         0  
283 731 50       1605 if (is_infinite($x)) { return $x; }
  0         0  
284 731 50       1602 if (is_infinite($y)) { return $y; }
  0         0  
285              
286 731         1658 return _xykk_to_n($self, $x,$y);
287             }
288              
289             # exact
290             sub rect_to_n_range {
291 1     1 1 8 my ($self, $x1,$y1, $x2,$y2) = @_;
292              
293 1         3 $x1 = round_nearest ($x1);
294 1         3 $y1 = round_nearest ($y1);
295 1         80 $x2 = round_nearest ($x2);
296 1         4 $y2 = round_nearest ($y2);
297 1 50       4 ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
298 1 50       3 ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
299             ### rect_to_n_range(): "$x1,$y1 to $x2,$y2"
300              
301 1 50 33     6 if ($x2 < 0 || $y2 < 0) {
302 0         0 return (1, 0);
303             }
304              
305 1         3 my $radix = $self->{'radix'};
306              
307 1         4 my ($power, $level) = round_down_pow (max($x2,$y2), $radix);
308 1 50       3 if (is_infinite($level)) {
309 0         0 return (0, $level);
310             }
311              
312 1         3 my $n_power = $power * $power * $radix;
313 1         2 my $max_x = 0;
314 1         2 my $max_y = 0;
315 1         2 my $max_n = 0;
316 1         1 my $max_xk = 0;
317 1         2 my $max_yk = 0;
318              
319 1         2 my $min_x = 0;
320 1         2 my $min_y = 0;
321 1         1 my $min_n = 0;
322 1         2 my $min_xk = 0;
323 1         3 my $min_yk = 0;
324              
325             # l<=c
326             # l>c2 or h-1
327             # l>c2 or h<=c1
328             # so does overlap if
329             # l<=c2 and h>c1
330             #
331 1         2 my $radix_minus_1 = $radix - 1;
332             my $overlap = sub {
333 5     5   17 my ($c,$ck,$digit, $c1,$c2) = @_;
334 5 100       12 if ($ck & 1) {
335 1         2 $digit = $radix_minus_1 - $digit;
336             }
337             ### overlap consider: "inv".($ck&1)."digit=$digit ".($c+$digit*$power)."<=c<".($c+($digit+1)*$power)." cf $c1 to $c2 incl"
338 5   66     27 return ($c + $digit*$power <= $c2
339             && $c + ($digit+1)*$power > $c1);
340 1         6 };
341              
342 1         4 while ($level-- >= 0) {
343             ### $power
344             ### $n_power
345             ### $max_n
346             ### $min_n
347             {
348 1         2 my $digit;
349 1         4 for ($digit = $radix_minus_1; $digit > 0; $digit--) {
350 2 100       5 last if &$overlap ($max_y,$max_yk,$digit, $y1,$y2);
351             }
352 1         11 $max_n += $n_power * $digit;
353 1         3 $max_xk ^= $digit;
354 1 50       3 if ($max_yk&1) { $digit = $radix_minus_1 - $digit; }
  0         0  
355 1         2 $max_y += $power * $digit;
356             ### max y digit (complemented): $digit
357             ### $max_y
358             ### $max_n
359             }
360             {
361 1         1 my $digit;
  1         1  
  1         8  
362 1         4 for ($digit = 0; $digit < $radix_minus_1; $digit++) {
363 1 50       3 last if &$overlap ($min_y,$min_yk,$digit, $y1,$y2);
364             }
365 1         3 $min_n += $n_power * $digit;
366 1         2 $min_xk ^= $digit;
367 1 50       4 if ($min_yk&1) { $digit = $radix_minus_1 - $digit; }
  0         0  
368 1         2 $min_y += $power * $digit;
369             ### min y digit (complemented): $digit
370             ### $min_y
371             ### $min_n
372             }
373              
374 1         2 $n_power = int($n_power/$radix);
375             {
376 1         2 my $digit;
377 1         3 for ($digit = $radix_minus_1; $digit > 0; $digit--) {
378 1 50       2 last if &$overlap ($max_x,$max_xk,$digit, $x1,$x2);
379             }
380 1         3 $max_n += $n_power * $digit;
381 1         2 $max_yk ^= $digit;
382 1 50       3 if ($max_xk&1) { $digit = $radix_minus_1 - $digit; }
  1         2  
383 1         2 $max_x += $power * $digit;
384             ### max x digit (complemented): $digit
385             ### $max_x
386             ### $max_n
387             }
388             {
389 1         2 my $digit;
  1         2  
  1         1  
390 1         3 for ($digit = 0; $digit < $radix_minus_1; $digit++) {
391 1 50       3 last if &$overlap ($min_x,$min_xk,$digit, $x1,$x2);
392             }
393 1         3 $min_n += $n_power * $digit;
394 1         2 $min_yk ^= $digit;
395 1 50       10 if ($min_xk&1) { $digit = $radix_minus_1 - $digit; }
  0         0  
396 1         3 $min_x += $power * $digit;
397             ### min x digit (complemented): $digit
398             ### $min_x
399             ### $min_n
400             }
401              
402 1         3 $power = int($power/$radix);
403 1         3 $n_power = int($n_power/$radix);
404             }
405             ### is: "$min_n at $min_x,$min_y to $max_n at $max_x,$max_y"
406 1         7 return ($min_n, $max_n);
407             }
408              
409             #------------------------------------------------------------------------------
410             # levels
411              
412 5     5   2083 use Math::PlanePath::ZOrderCurve;
  5         14  
  5         286  
413             *level_to_n_range = \&Math::PlanePath::ZOrderCurve::level_to_n_range;
414             *n_to_level = \&Math::PlanePath::ZOrderCurve::n_to_level;
415              
416             #------------------------------------------------------------------------------
417             1;
418             __END__