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.

This task has been flagged for clarification. Code on this page in its current state may be flagged incorrect once this task has been clarified. See this page's Talk page for discussion.

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 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 modified from 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 r = (sqrt(a)*b + m) / n,
       rem = r
    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)