File Coverage

blib/lib/Math/NumSeq/NumAronson.pm
Criterion Covered Total %
statement 70 74 94.5
branch 19 22 86.3
condition 6 9 66.6
subroutine 15 15 100.0
pod 4 4 100.0
total 114 124 91.9


line stmt bran cond sub pod time code
1             # Copyright 2010, 2011, 2012, 2013, 2014 Kevin Ryde
2              
3             # This file is part of Math-NumSeq.
4             #
5             # Math-NumSeq 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-NumSeq 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-NumSeq. If not, see .
17              
18             package Math::NumSeq::NumAronson;
19 8     8   177 use 5.004;
  8         34  
  8         726  
20 8     8   52 use strict;
  8         16  
  8         333  
21              
22 8     8   49 use vars '$VERSION', '@ISA';
  8         15  
  8         559  
23             $VERSION = 71;
24 8     8   56 use Math::NumSeq;
  8         199  
  8         440  
25             @ISA = ('Math::NumSeq');
26              
27             # uncomment this to run the ### lines
28             #use Devel::Comments;
29              
30              
31             # use constant name => Math::NumSeq::__('Numerical Aronson');
32 8     8   60 use constant description => Math::NumSeq::__('Numerical version of Aronson\'s sequence');
  8         22  
  8         42  
33 8     8   43 use constant values_min => 1;
  8         23  
  8         484  
34 8     8   47 use constant characteristic_increasing => 1;
  8         13  
  8         422  
35 8     8   42 use constant characteristic_integer => 1;
  8         19  
  8         396  
36 8     8   44 use constant i_start => 1;
  8         14  
  8         419  
37              
38             # cf A080596 - a(1)=1
39             # A079253 - even
40             # A081023 - lying
41             # A014132 - lying opposite parity
42 8     8   44 use constant oeis_anum => 'A079000';
  8         13  
  8         6245  
43              
44             # a(9*2^k - 3 + j) = 12*2^k - 3 + (3/2)*j + (1/2)*abs(j)
45             # where k>=0 and -3*2^k <= j < 3*2^k
46             # step
47             # a(n+1) - 2*a(n) + a(n-1) = 1 if n=9*2^k-3, k>=0
48             # = -1 if n = 2 and 3*2^k-3, k>=1
49             # = 0 otherwise.
50             #
51             # lying
52             # g(3*2^k-1 + j) = 2*2^(k+1)-1 + (3/2)*j + (1/2)*abs(j)
53             # where -2^k <= j < 2^k and k>0
54             #
55             # then lying is d(n)=g(n+1)-1 n>=1
56              
57             sub rewind {
58 3     3 1 931 my ($self) = @_;
59 3         23 $self->{'i'} = $self->i_start;
60 3         8 $self->{'pow'} = 0;
61 3         19 $self->{'j'} = -1;
62             }
63              
64             sub next {
65 90     90 1 729 my ($self) = @_;
66 90         139 my $pow = $self->{'pow'};
67 90         108 my $j = ++ $self->{'j'};
68              
69 90 100       199 if ($pow == 0) {
    100          
70             # low special cases initial 1,4,
71 6 100       20 if ($j < 2) {
72 4         16 return ($self->{'i'}++, $j*3 + 1);
73             }
74 2         7 $pow = $self->{'pow'} = 1; # 2**k for k=0
75 2         4 $j = $self->{'j'} = -3; # -3*(2**k) for k=0
76             } elsif ($j >= 3 * $pow) {
77 6         8 $pow = ($self->{'pow'} *= 2);
78 6         9 $j = $self->{'j'} = -3 * $pow;
79             }
80              
81             ### assert: -3 * $pow <= $j
82             ### assert: $j < 3 * $pow
83              
84 86         241 return ($self->{'i'}++, 12*$pow - 3 + (3*$j + abs($j))/2);
85             }
86              
87             # i = 9*2^k - 3 + j
88             # base at j=-3*2^k
89             # is i >= 9*2^k - 3 - 3*2^k
90             # i >= 6*2^k - 3
91             # i+3 >= 6*2^k
92             # (i+3)/6 >= 2^k
93             # 2^k <= (i+3)/6
94             # then i = 9*2^k - 3 + j
95             # j = i - 9*2^k + 3
96             #
97             sub ith {
98 138     138 1 281 my ($self, $i) = @_;
99             ### NumAronson ith(): $i
100              
101             # special cases ith(1)=1, ith(2)=4
102 138 100       236 if ($i <= 2) {
103 9 100       17 if ($i < 0) {
104 3         7 return undef;
105             } else {
106 6         18 return $i*$i;
107             }
108             }
109              
110 129         366 my ($pow, $k) = _round_down_pow (int(($i+3)/6), 2);
111 129         174 my $j = $i - 9*$pow + 3;
112              
113             ### round down for k: ($i+3)/6
114             ### $k
115             ### $pow
116             ### $j
117             ### assert: $k >= 0
118             ### assert: -3 * 2 ** $k <= $j
119             ### assert: $j < 3 * 2 ** $k
120              
121 129         352 return 12*$pow - 3 + (3*$j + abs($j))/2;
122             }
123              
124             # value = 12*2^k - 3 + (3/2)*j + (1/2)*abs(j)
125             # minimum j=-3*2^k
126             # value = 12*2^k - 3 + (3*j + abs(j))/2
127             # = 12*2^k - 3 + (3*-3*2^k + 3*2^k)/2
128             # = 12*2^k - 3 + (-9 + 3)*2^k/2
129             # = 12*2^k - 3 + -6*2^k/2
130             # = 12*2^k - 3 + -3*2^k
131             # = 9*2^k - 3
132             # value >= 9*2^k - 3
133             # value+3 >= 9*2^k
134             # 9*2^k <= value+3
135             # 2^k <= (value+3)/9
136             #
137             # from which
138             # value = 12*2^k - 3 + (3*j + abs(j))/2
139             # (3*j + abs(j))/2 = value - 12*2^k + 3
140             # 3*j + abs(j) = 2*(value - 12*2^k + 3)
141             #
142             # j>=0 4*j = a, j=a/4
143             # j<0 3*$j-$j = 2*$j = a, j=a/2
144             #
145             sub pred {
146 182     182 1 1084 my ($self, $value) = @_;
147             ### NumAronson pred(): $value
148              
149             # special cases pred(1) true, pred(4) true
150 182 100       362 if ($value < 6) {
151 12   100     64 return ($value == 1 || $value == 4);
152             }
153              
154 170         452 my $k = _round_down_pow (int(($value+3)/9), 2);
155 170         237 my $pow = 2**$k;
156 170         285 my $aj = 2*($value - 12*$pow + 3);
157 170 100       579 return ($aj % ($aj < 0 ? 2 : 4)) == 0;
158             }
159              
160             #------------------------------------------------------------------------------
161             # generic
162              
163             # return ($pow, $exp) with $pow = $base**$exp <= $n,
164             # the next power of $base at or below $n
165             #
166             sub _round_down_pow {
167 2859     2859   4606 my ($n, $base) = @_;
168             ### _round_down_pow(): "$n base $base"
169              
170 2859 100       7219 if ($n < $base) {
171 644         1673 return (1, 0);
172             }
173              
174             # Math::BigInt and Math::BigRat overloaded log() return NaN, use integer
175             # based blog()
176 2215 50 33     6541 if (ref $n && ($n->isa('Math::BigInt') || $n->isa('Math::BigRat'))) {
      66        
177 13         45 my $exp = $n->copy->blog($base);
178 13         9873 return (Math::BigInt->new(1)->blsft($exp,$base),
179             $exp);
180             }
181              
182 2202         5409 my $exp = int(log($n)/log($base));
183 2202         3380 my $pow = $base**$exp;
184              
185             ### n: ref($n)." $n"
186             ### exp: ref($exp)." $exp"
187             ### pow: ref($pow)." $pow"
188              
189             # check how $pow actually falls against $n, not sure should trust float
190             # rounding in log()/log($base)
191             # Crib: $n as first arg in case $n==BigFloat and $pow==BigInt
192 2202 50       13455 if ($n < $pow) {
    50          
193             ### hmm, int(log) too big, decrease...
194 0         0 $exp -= 1;
195 0         0 $pow = $base**$exp;
196             } elsif ($n >= $base*$pow) {
197             ### hmm, int(log) too small, increase...
198 0         0 $exp += 1;
199 0         0 $pow *= $base;
200             }
201 2202         7030 return ($pow, $exp);
202             }
203              
204             1;
205             __END__