++ed by:
Kevin Ryde
and 1 contributors

# NAME

Math::NumSeq::SternDiatomic -- Stern's diatomic sequence

# SYNOPSIS

`````` use Math::NumSeq::SternDiatomic;
my \$seq = Math::NumSeq::SternDiatomic->new;
my (\$i, \$value) = \$seq->next;``````

# DESCRIPTION

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 last copied value and the first value of the next row, which is always 1. This means the last value in each row increments 1,2,3,4,5,etc.

## Odd and Even

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``````

# FUNCTIONS

See "FUNCTIONS" in Math::NumSeq for behaviour common to all sequence classes.

`\$seq = Math::NumSeq::SternDiatomic->new ()`

Create and return a new sequence object.

## Random Access

`\$value = \$seq->ith(\$i)`

Return the `\$i`'th value of the sequence.

`(\$v0, \$v1) = \$seq->ith_pair(\$i)`

Return two values `ith(\$i)` and `ith(\$i+1)` from the sequence. As described below ("Ith Pair") two values can be calculated with the same work as one.

`\$bool = \$seq->pred(\$value)`

Return true if `\$value` occurs in the sequence, which means simply integer `\$value>=0`.

# FORMULAS

## Next

The sequence is iterated using a method by Moshe Newman in

"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]    q = seq[i+1]
initially p=seq[0]=0 and q=seq[1]=1

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.

``````    division p = q*floor(p/q) + rem      where rem = (p mod q)
then
p_next/q_next = 1 / (1 + 2*floor(p/q) - p/q)    per Newman
= q / (2*q*floor(p/q) + q - p)
= q / (2*(p - rem) + q - p)
= q / (p+q - 2*rem)               using p,q``````

`seek_to_i()` is implemented by calculating new p,q values with an Ith Pair per below.

## Ith Pair

The two sequence values at an arbitrary i,i+1 can be calculated from the bits of i,

``````    p = 0
q = 1
for each bit of i from high to low
if bit=1 then p += q
if bit=0 then q += p
return p,q      # are ith(i) and ith(i+1)``````

For example i=6 is binary "110" so

``````                         p,q
---
initial               0,1
high bit=1 so p+=q    1,1
next bit=1 so p+=q    2,1
low  bit=0 so q+=p    2,3   is ith(6),ith(7)``````

This is the same as the Calkin-Wilf tree descent, per "Calkin-Wilf Tree" in Math::PlanePath::RationalsTree. Its X/Y fractions are successive Stern diatomic sequence values.

## Ith Alone

If only a single ith() value is desired then the bits of i can be taken from low to high with the same loop as above. In that case p=ith(i) but q is not ith(i+1). Any low zero bits can be ignored for this method since initial p=0 means their steps q+=p do nothing.

Math::NumSeq

Math::PlanePath::RationalsTree