Sturmian word

Sturmian word is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

A Sturmian word is a binary sequence, finite or infinite, that makes up the cutting sequence for a positive real number x, as shown in the picture.

Example Sturmian word when x = 0.618..., the golden ratio.

The Sturmian word can be computed thus as an algorithm:

• If $x>1$, then it is the inverse of the Sturmian word for $1/x$. So we have reduced to the case of $0<x\leq 1$.
• Iterate over $floor(1x),floor(2x),floor(3x),\dots$
• If $kx$ is an integer, then the program terminates. Else, if $floor((k-1)x)=floor(kx)$, then the program outputs 0, else, it outputs 10.

The problem:

• Given a positive rational number ${\frac {m}{n}}$, specified by two positive integers $m,n$, output its entire Sturmian word.
• Given a quadratic real number ${\frac {b{\sqrt {a}}+m}{n}}>0$, specified by integers $a,b,m,n$, where $a$ is not a perfect square, output the first $k$ letters of Sturmian words when given a positive number $k$.

(If the programming language can represent infinite data structures, then that works too.)

A simple check is to do this for the inverse golden ratio ${\frac {{\sqrt {5}}-1}{2}}$, that is, $a=5,b=1,m=-1,n=2$, which would just output the Fibonacci word.

Stretch goal: calculate the Sturmian word for other kinds of definable real numbers, such as cubic roots.

The key difficulty is accurately calculating $floor(k{\sqrt {a}})$ for large $k$. Floating point arithmetic would lose precision. One can either do this simply by directly searching for some integer $a'$ such that $a'^{2}\leq k^{2}a<(a'+1)^{2}$, or by more trickly methods, such as the continued fraction approach.

First calculate the continued fraction convergents to ${\sqrt {a}}$. Let ${\frac {m}{n}}$ be a convergent to ${\sqrt {a}}$, such that $n\geq k$, then since the convergent sequence is the best rational approximant for denominators up to that point, we know for sure that, if we write out ${\frac {0}{k}},{\frac {1}{k}},\dots$, the sequence would stride right across the gap $(m/n,2x-m/n)$. Thus, we can take the largest $l$ such that $l/k\leq m/n$, and we would know for sure that $l=floor(k{\sqrt {a}})$.

In summary, $floor(k{\sqrt {a}})=floor(mk/n)$ where $m/n$ is the first continued fraction approximant to ${\sqrt {a}}$ with a denominator $n\geq k$

Phix

Rational part translated from Python, quadratic part a modified copy of the Continued fraction convergents task.

with javascript_semantics
function sturmian_word(integer m, n)
    if m > n then
        return sq_sub('0'+'1',sturmian_word(n,m))
    end if
    string res = ""
    integer k = 1, prev = 0
    while remainder(k*m,n) do
        integer curr = floor(k*m/n)
        res &= iff(prev=curr?"0":"10")
        prev = curr
        k += 1
    end while
    return res
end function

function fibWord(integer n)
    string Sn_1 = "0",
             Sn = "01",
            tmp = ""
    for i=2 to n do
        tmp = Sn
        Sn &= Sn_1
        Sn_1 = tmp
    end for
    return Sn
end function

string fib = fibWord(7),
  sturmian = sturmian_word(13, 21)
assert(fib[1..length(sturmian)] == sturmian)
printf(1," %s <== 13/21\n",sturmian)

function cfck(integer a, b, m, n, k)
    -- return the kth convergent
    sequence p = {0, 1},
             q = {1, 0}
    atom rem = (sqrt(a)*b + m) / n
    for i=1 to k do
        integer whole = trunc(rem),
                pn = whole * p[-1] + p[-2],
                qn = whole * q[-1] + q[-2]
        p &= pn
        q &= qn
        rem = 1/(rem-whole)
    end for
    return {p[$],q[$]}
end function

integer {m,n} = cfck(5, 1, -1, 2, 8) -- (inverse golden ratio)
printf(1," %s <== 1/phi (8th convergent)\n",sturmian_word(m,n))

Output:

 01001010010010100101001001010010 <== 13/21
 01001010010010100101001001010010 <== 1/phi (8th convergent)

Python

For rational numbers:

def sturmian_word(m, n):
    sturmian = ""
    def invert(string):
      return ''.join(list(map(lambda b: {"0":"1", "1":"0"}[b], string)))
    if m > n:
      return invert(sturmian_word(n, m))

    k = 1
    while True:
        current_floor = int(k * m / n)
        previous_floor = int((k - 1) * m / n)
        if k * m % n == 0:
            break
        if previous_floor == current_floor:
            sturmian += "0"
        else:
            sturmian += "10"
        k += 1
    return sturmian

Checking that it works on the finite Fibonacci word:

def fibWord(n):
    Sn_1 = "0"
    Sn = "01"
    tmp = ""
    for i in range(2, n + 1):
        tmp = Sn
        Sn += Sn_1
        Sn_1 = tmp
    return Sn

fib = fibWord(10)
sturmian = sturmian_word(13, 21)
assert fib[:len(sturmian)] == sturmian
print(sturmian)

# Output:
# 01001010010010100101001001010010

Wren

Library: Wren-assert

The 'rational number' function is a translation of the Python entry.

The 'quadratic number' function is a modified version of the one used in the Continued fraction convergents task.

import "./assert" for Assert

var SturmianWordRat = Fn.new { |m, n|
    if (m > n) return SturmianWordRat.call(n, m).reduce("") { |acc, c| 
        return acc + (c == "0" ? "1" : "0") 
    }
    var sturmian = ""
    var k = 1
    while (k * m % n != 0) {
        var currFloor = (k * m / n).floor
        var prevFloor = ((k - 1) * m / n).floor
        sturmian = sturmian + (prevFloor == currFloor ? "0" : "10")
        k  = k + 1
    }
    return sturmian
}

var SturmianWordQuad = Fn.new { |a, b, m, n, k|
    var p = [0, 1]
    var q = [1, 0]
    var rem = (a.sqrt * b + m) / n
    for (i in 1..k) {
        var whole = rem.truncate
        var frac  = rem.fraction
        var pn = whole * p[-1] + p[-2]
        var qn = whole * q[-1] + q[-2]
        p.add(pn)
        q.add(qn)
        rem = 1 / frac
    }
    return SturmianWordRat.call(p[-1], q[-1])
}

var fibWord = Fn.new { |n|
    var sn1 = "0"
    var sn  = "01"
    for (i in 2..n) {
        var tmp = sn
        sn = sn + sn1
        sn1 = tmp
    }
    return sn
}

var fib = fibWord.call(10)
var sturmian = SturmianWordRat.call(13, 21)
Assert.equal(fib[0...sturmian.count], sturmian)
System.print("%(sturmian) from rational number 13/21")
var sturmian2 = SturmianWordQuad.call(5, 1, -1, 2, 8)
Assert.equal(sturmian, sturmian2)
System.print("%(sturmian2) from quadratic number (√5 - 1)/2 (k = 8)")

Output:

01001010010010100101001001010010 from rational number 13/21
01001010010010100101001001010010 from quadratic number (√5 - 1)/2 (k = 8)