Math::NumSeq::SternDiatomic -- Stern's diatomic sequence
use Math::NumSeq::SternDiatomic; my $seq = Math::NumSeq::SternDiatomic->new; my ($i, $value) = $seq->next;
This is Moritz Stern's diatomic sequence
0, 1, 1, 2, 1, 3, 2, 3, ... starting i=0
It's constructed by successive levels with a recurrence
D(0) = 0 D(1) = 1 D(2*i) = D(i) D(2*i+1) = D(i) + D(i+1)
So the sequence is extended by copying the previous level to the next spead out to even indices, and at the odd indices fill in the sum of adjacent terms,
0, i=0 1, i=1 1, 2, i=2 to 3 1, 3, 2, 3, i=4 to 7 1,4,3,5,2,5,3,4, i=8 to 15
For example the i=4to7 row is a copy of the preceding row values 1,2 with sums 1+2 and 2+1 interleaved. For the new value at the end of each row the sum wraps around so as to take the copied value on the left and the first value of the next row (which is always 1).
The sequence makes a repeating pattern even,odd,odd,
0, 1, 1, 2, 1, 3, 2, 3 E O O E O O E ...
This can be seen from the copying in the recurrence above. For example the i=8 to 15 row copying to i=16 to 31,
O . E . O . O . E . O . O . E . spread O O E O O E O O sum adjacent
Adding adjacent terms odd+even and even+odd are both odd. Adding adjacent odd+odd gives even. So the pattern E,O,O in the original row when spread and added gives E,O,O again in the next row.
See "FUNCTIONS" in Math::NumSeq for behaviour common to all sequence classes.
$seq = Math::NumSeq::SternDiatomic->new ()
Create and return a new sequence object.
$value = $seq->ith($i)
$i'th value of the sequence.
$bool = $seq->pred($value)
Return true if
$value occurs in the sequence, which means simply integer
The sequence is iterated using a method by Moshe Newman.
"Recounting the Rationals, Continued", answers to problem 10906 posed by Donald E. Knuth, C. P. Rupert, Alex Smith and Richard Stong, American Mathematical Monthly, volume 110, number 7, Aug-Sep 2003, pages 642-643, http://www.jstor.org/stable/3647762
Two successive sequence values are maintained and are advanced by a simple operation.
p = seq[i] initially p=0 = seq q = seq[i+1] q=1 = seq p_next = seq[i+1] = q q_next = seq[i+2] = p+q - 2*(p mod q) where the mod operation rounds towards zero 0 <= (p mod q) <= q-1
The form by Newman uses a floor operation. This suits expressing the iteration in terms of a rational x[i]=p/q.
p_next 1 ------ = ---------------------- q_next 1 + 2*floor(p/q) - p/q
For separate p,q a little rearrangement gives it in terms of the remainder p mod q.
p = q*floor(p/q) + rem where rem = (p mod q) so p_next/q_next = q / (2*q*floor(p/q) + q - p) = q / (2*(p - rem) + q - p) = q / (p+q - 2*rem)
seek_to_i() is implemented by calculating new p,q values with
ith(i) per below.
The sequence at an arbitrary i can be calculated from the bits of i,
p = 0 q = 1 for each bit of i from low to high if bit=1 then p += q if bit=0 then q += p return p
For example i=6 is binary "110" so p,q starts 0,1 then low bit=0 q+=p leaves that unchanged as 0,1. Next bit=1 p+=q gives 1,1 and high bit=1 gives 2,1 for result 2.
Any low zero bits can be ignored since initial p=0 means their steps q+=p do nothing. The current code doesn't use this since it's not expected i would usually have many trailing zeros and the q+=p saved won't be particularly slow.
Copyright 2011, 2012, 2013 Kevin Ryde
Math-NumSeq is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3, or (at your option) any later version.
Math-NumSeq is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with Math-NumSeq. If not, see <http://www.gnu.org/licenses/>.