Sturmian word: Difference between revisions

From Rosetta Code
Content added Content deleted
(added Raku programming solution)
(Added various BASIC dialects (BASIC256,FreeBASIC and Yabasic))
Line 25: Line 25:
In summary, <math>floor(k\sqrt a) = floor(mk/n)</math> where <math>m/n</math> is the first continued fraction approximant to <math>\sqrt a</math> with a denominator <math>n \geq k</math>
In summary, <math>floor(k\sqrt a) = floor(mk/n)</math> where <math>m/n</math> is the first continued fraction approximant to <math>\sqrt a</math> with a denominator <math>n \geq k</math>



=={{header|BASIC}}==
==={{header|BASIC256}}===
{{trans|FreeBASIC}}
<syntaxhighlight lang="vbnet">fib = fibWord(10)
sturmian = sturmian_word(13, 21)
if left(fib, length(sturmian)) = sturmian then print sturmian
end

function invert(cadena)
for i = 1 to length(cadena)
b = mid(cadena, i, 1)
if b = "0" then
inverted += "1"
else
if b = "1" then
inverted += "0"
end if
end if
next
return inverted
end function

function sturmian_word(m, n)
sturmian = ""
if m > n then 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 mod n = 0 then exit while
if previous_floor = current_floor then
sturmian += "0"
else
sturmian += "10"
end if
k += 1
end while
return sturmian
end function

function fibWord(n)
Sn_1 = "0"
Sn = "01"

for i = 2 to n
tmp = Sn
Sn += Sn_1
Sn_1 = tmp
next
return Sn
end function</syntaxhighlight>
{{out}}
<pre>01001010010010100101001001010010</pre>

==={{header|FreeBASIC}}===
{{trans|Python}}
<syntaxhighlight lang="vbnet">Function invert(cadena As String) As String
Dim As String inverted, b
For i As Integer = 1 To Len(cadena)
b = Mid(cadena, i, 1)
inverted += Iif(b = "0", "1", "0")
Next
Return inverted
End Function

Function sturmian_word(m As Integer, n As Integer) As String
Dim sturmian As String
If m > n Then Return invert(sturmian_word(n, m))
Dim As Integer k = 1
Do
Dim As Integer current_floor = Int(k * m / n)
Dim As Integer previous_floor = Int((k - 1) * m / n)
If k * m Mod n = 0 Then Exit Do
sturmian += Iif(previous_floor = current_floor, "0", "10")
k += 1
Loop
Return sturmian
End Function

Function fibWord(n As Integer) As String
Dim As String Sn_1 = "0"
Dim As String Sn = "01"
For i As Integer = 2 To n
Dim As String tmp = Sn
Sn += Sn_1
Sn_1 = tmp
Next
Return Sn
End Function

Dim As String fib = fibWord(10)
Dim As String sturmian = sturmian_word(13, 21)
If Left(fib, Len(sturmian)) = sturmian Then Print sturmian

Sleep</syntaxhighlight>
{{out}}
<pre>01001010010010100101001001010010</pre>

==={{header|Yabasic}}===
{{trans|FreeBASIC}}
<syntaxhighlight lang="vbnet">fib$ = fibword$(10)
sturmian$ = sturmian_word$(13,21)
if left$(fib$,len(sturmian$)) = sturmian$ print sturmian$
end
sub invert$(cadena$)
for i = 1 to len(cadena$)
b$ = mid$(cadena$,i,1)
if b$ = "0" then
inverted$ = inverted$ +"1"
elsif b$ = "1" then
inverted$ = inverted$ +"0"
fi
next i
return inverted$
end sub
sub sturmian_word$(m,n)
sturmian$ = ""
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 mod(k*m, n) = 0 break
if previous_floor = current_floor then
sturmian$ = sturmian$ +"0"
else
sturmian$ = sturmian$ +"10"
fi
k = k +1
end while
return sturmian$
end sub
sub fibword$(n)
sn_1$ = "0"
sn$ = "01"
for i = 2 to n
tmp$ = sn$
sn$ = sn$ + sn_1$
sn_1$ = tmp$
next i
return sn$
end sub</syntaxhighlight>
{{out}}
<pre>01001010010010100101001001010010</pre>


=={{header|Julia}}==
=={{header|Julia}}==

Revision as of 00:26, 22 February 2024

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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x>1} , then it is the inverse of the Sturmian word for Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 1/x} . So we have reduced to the case of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 0 < x \leq 1} .
  • Iterate over Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle floor(1x), floor(2x), floor(3x), \dots }
  • If Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle kx } is an integer, then the program terminates. Else, if Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle floor((k-1)x) = floor(kx)} , then the program outputs 0, else, it outputs 10.

The problem:

  • Given a positive rational number Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac mn} , specified by two positive integers Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m, n} , output its entire Sturmian word.
  • Given a quadratic real number Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{b\sqrt{a} + m}{n} > 0} , specified by integers Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a, b, m, n } , where Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a} is not a perfect square, output the first Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} letters of Sturmian words when given a positive number Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\sqrt{5}-1}{2}} , that is, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle floor(k\sqrt a) } for large Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} . Floating point arithmetic would lose precision. One can either do this simply by directly searching for some integer Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a'} such that Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a'^2 \leq k^2a < (a'+1)^2} , or by more trickly methods, such as the continued fraction approach.

First calculate the continued fraction convergents to Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sqrt a} . Let Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac mn } be a convergent to Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sqrt a} , such that Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{0}{k}, \frac{1}{k}, \dots} , the sequence would stride right across the gap Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle (m/n, 2x - m/n)} . Thus, we can take the largest Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle l} such that Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle l/k \leq m/n} , and we would know for sure that Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle l = floor(k\sqrt a)} .

In summary, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle floor(k\sqrt a) = floor(mk/n)} where Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m/n} is the first continued fraction approximant to Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sqrt a} with a denominator Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle n \geq k}


BASIC

BASIC256

Translation of: FreeBASIC
fib = fibWord(10)
sturmian = sturmian_word(13, 21)
if left(fib, length(sturmian)) = sturmian then print sturmian
end

function invert(cadena)
	for i = 1 to length(cadena)
		b = mid(cadena, i, 1)
		if b = "0" then
			inverted += "1"
		else
			if b = "1" then
				inverted += "0"
			end if
		end if
	next
	return inverted
end function

function sturmian_word(m, n)
	sturmian = ""
	if m > n then 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 mod n = 0 then exit while
		if previous_floor = current_floor then
			sturmian += "0"
		else
			sturmian += "10"
		end if
		k += 1
	end while
	return sturmian
end function

function fibWord(n)
	Sn_1 = "0"
	Sn = "01"

	for i = 2 to n
		tmp = Sn
		Sn += Sn_1
		Sn_1 = tmp
	next
	return Sn
end function
Output:
01001010010010100101001001010010

FreeBASIC

Translation of: Python
Function invert(cadena As String) As String
    Dim As String inverted, b
    For i As Integer = 1 To Len(cadena)
        b = Mid(cadena, i, 1)
        inverted += Iif(b = "0", "1", "0")
    Next
    Return inverted
End Function

Function sturmian_word(m As Integer, n As Integer) As String
    Dim sturmian As String
    If m > n Then Return invert(sturmian_word(n, m))
    
    Dim As Integer k = 1
    Do
        Dim As Integer current_floor = Int(k * m / n)
        Dim As Integer previous_floor = Int((k - 1) * m / n)
        If k * m Mod n = 0 Then Exit Do
        sturmian += Iif(previous_floor = current_floor, "0", "10")
        k += 1
    Loop
    Return sturmian
End Function

Function fibWord(n As Integer) As String
    Dim As String Sn_1 = "0"
    Dim As String Sn = "01"
    
    For i As Integer = 2 To n
        Dim As String tmp = Sn
        Sn += Sn_1
        Sn_1 = tmp
    Next
    Return Sn
End Function

Dim As String fib = fibWord(10)
Dim As String sturmian = sturmian_word(13, 21)
If Left(fib, Len(sturmian)) = sturmian Then Print sturmian

Sleep
Output:
01001010010010100101001001010010

Yabasic

Translation of: FreeBASIC
fib$ = fibword$(10)
sturmian$ = sturmian_word$(13,21)
if left$(fib$,len(sturmian$)) = sturmian$  print sturmian$
end
sub invert$(cadena$)
    for i = 1 to len(cadena$)
        b$ = mid$(cadena$,i,1)
        if b$ = "0" then
            inverted$ = inverted$ +"1"
        elsif b$ = "1" then
            inverted$ = inverted$ +"0"
        fi
    next i
    return inverted$
end sub
sub sturmian_word$(m,n)
    sturmian$ = ""
    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 mod(k*m, n) = 0  break
        if previous_floor = current_floor then
            sturmian$ = sturmian$ +"0"
        else
            sturmian$ = sturmian$ +"10"
        fi
        k = k +1
    end while
    return sturmian$
end sub
sub fibword$(n)
    sn_1$ = "0"
    sn$ = "01"
    for i = 2 to n
        tmp$ = sn$
        sn$ = sn$ + sn_1$
        sn_1$ = tmp$
    next i
    return sn$
end sub
Output:
01001010010010100101001001010010

Julia

Translation of: Phix
function sturmian_word(m, n)
    m > n && return replace(sturmian_word(n, m), '0' => '1', '1' => '0')
    res = ""
    k, prev = 1, 0
    while rem(k * m, n) > 0
        curr = (k * m) ÷ n
        res *= prev == curr ? "0" : "10"
        prev = curr
        k += 1
    end
    return res
end

function fibWord(n)
    Sn_1, Sn, tmp = "0", "01", ""
    for _ in 2:n 
        tmp = Sn
        Sn *= Sn_1
        Sn_1 = tmp
    end
    return Sn
end

const fib = fibWord(7)
const sturmian = sturmian_word(13, 21)
@assert fib[begin:length(sturmian)] == sturmian
println(" $sturmian <== 13/21")

""" return the kth convergent """
function cfck(a, b, m, n, k)
    p = [0, 1]
    q = [1, 0]
    r = (sqrt(a) * b + m) / n
    for _ in 1:k
        whole = Int(trunc(r))
        pn = whole * p[end] + p[end-1]
        qn = whole * q[end] + q[end-1]
        push!(p, pn)
        push!(q, qn)
        r = 1/(r - whole)
    end
    return [p[end], q[end]]
end

println(" $(sturmian_word(cfck(5, 1, -1, 2, 8)...)) <== 1/phi (8th convergent golden ratio)")
Output:
 01001010010010100101001001010010 <== 13/21
 01001010010010100101001001010010 <== 1/phi (8th convergent golden ratio)


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

Raku

Translation of: Julia
# 20240215 Raku programming solution

sub chenfoxlyndonfactorization(Str $s) {
 sub sturmian-word($m, $n) {
   return sturmian-word($n, $m).trans('0' => '1', '1' => '0') if $m > $n; 
   my ($res, $k, $prev) = '', 1, 0;
   while ($k * $m) % $n > 0 {
      my $curr = ($k * $m) div $n;
      $res ~= $prev == $curr ?? '0' !! '10';
      $prev = $curr;
      $k++;
   }
   return $res;
}

sub fib-word($n) {
   my ($Sn_1, $Sn) = '0', '01';
   for 2..$n { ($Sn, $Sn_1) = ($Sn~$Sn_1, $Sn) }
   return $Sn;
}

my $fib = fib-word(7);
my $sturmian = sturmian-word(13, 21);
say "{$sturmian} <== 13/21" if $fib.substr(0, $sturmian.chars) eq $sturmian;

sub cfck($a, $b, $m, $n, $k) {
   my (@p, @q) := [0, 1], [1, 0]; 
   my $r = (sqrt($a) * $b + $m) / $n;
   for ^$k {
      my $whole = $r.Int;
      my ($pn, $qn) = $whole * @p[*-1] + @p[*-2], $whole * @q[*-1] + @q[*-2];
      @p.push($pn);
      @q.push($qn);
      $r = 1/($r - $whole);
   }
   return [@p[*-1], @q[*-1]];
}

my $cfck-result = cfck(5, 1, -1, 2, 8);
say "{sturmian-word($cfck-result[0], $cfck-result[1])} <== 1/phi (8th convergent golden ratio)";
Output:

Same as Julia example.

You may Attempt This Online!

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)