stern_brocot: Stern-Brocot Sequence

In numbers: Number-Theoretic Functions

Stern-Brocot Sequence

Description

The function generates the Stern-Brocot sequence up to length n.

Usage

  stern_brocot_seq(n)

Arguments

n	integer; length of the sequence.

Details

The Stern-Brocot sequence is a sequence S of natural numbers beginning with

1, 1, 2, 1, 3, 2, 3, 1, 4, 3, 5, 2, 5, 3, 4, 1, ...

defined with S[1] = S[2] = 1 and the following rules:

S[k] = S[k/2] if k is even
S[k] = S[(k-1)/2] + S[(k+1)/2] if k is not even

The Stern-Brocot has the remarkable properties that

(1) Consecutive values in this sequence are coprime;
(2) the list of rationals S[k+1]/S[k] (all in reduced form) covers all positive rational numbers once and once only.

Value

Returns a sequence of length n of natural numbers.

References

N. Calkin and H.S. Wilf. Recounting the rationals. The American Mathematical Monthly, Vol. 7(4), 2000.

Graham, Knuth, and Patashnik. Concrete Mathematics - A Foundation for Computer Science. Addison-Wesley, 1989.

See Also

fibonacci

Examples

( S <- stern_brocot_seq(92) )
# 1, 1, 2, 1, 3, 2, 3, 1, 4, 3, 5, 2, 5, 3, 4, 1,  5, 4, 7, 
# 3, 8, 5, 7, 2, 7, 5, 8, 3, 7, 4, 5, 1, 6, 5, 9, 4, 11, 7, 10, 
# 3, 11, 8, 13, 5, 12, 7, 9, 2, 9, 7, 12, 5, 13, 8, 11, 3, 10, 7, 11, 
# 4, 9, 5, 6, 1, 7, 6, 11, 5, 14, 9, 13, 4, 15, 11, 18, 7, 17, 10, 13, 
# 3, 14, 11, 19, 8, 21, 13, 18, 5, 17, 12, 19, 7, ...

table(S)
## S
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 17 18 19 21 
##  7  5  9  7 12  3 11  5  5  3  7  3  5  2  1  2  2  2  1 

which(S == 1)  # 1  2  4  8 16 32 64

## Not run: 
# Find the rational number p/q in S
# note that 1/2^n appears in position S[c(2^(n-1), 2^(n-1)+1)]
occurs <- function(p, q, s){
    # Find i such that (p, q) = s[i, i+1]
    inds <- seq.int(length = length(s)-1)
    inds <- inds[p == s[inds]]
    inds[q == s[inds + 1]]
}
p = 3; q = 7        # 3/7
occurs(p, q, S)     # S[28, 29]

'%//%' <- function(p, q) gmp::as.bigq(p, q)
n <- length(S)
S[1:(n-1)] %//% S[2:n]
## Big Rational ('bigq') object of length 91:
##  [1] 1     1/2  2     1/3   3/2   2/3   3     1/4   4/3   3/5   
## [11] 5/2   2/5  5/3   3/4   4     1/5   5/4   4/7   7/3   3/8   ...

as.double(S[1:(n-1)] %//% S[2:n])
## [1] 1.000000 0.500000 2.000000 0.333333 1.500000 0.666667 3.000000
## [8] 0.250000 1.333333 0.600000 2.500000 0.400000 1.666667 0.750000 ...

## End(Not run)

numbers documentation built on Dec. 29, 2022, 4:07 p.m.