I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

Minkowski question-mark function

From Rosetta Code
Minkowski question-mark function 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.

The Minkowski question-mark function converts the continued fraction representation [a0; a1, a2, a3, ...] of a number into a binary decimal representation in which the integer part a0 is unchanged and the a1, a2, ... become alternating runs of binary zeroes and ones of those lengths. The decimal point takes the place of the first zero.

Thus, ?(31/7) = 71/16 because 31/7 has the continued fraction representation [4;2,3] giving the binary expansion 4 + 0.01112.

Among its interesting properties is that it maps roots of quadratic equations, which have repeating continued fractions, to rational numbers, which have repeating binary digits.

The question-mark function is continuous and monotonically increasing, so it has an inverse.

  • Produce a function for ?(x).   Be careful: rational numbers have two possible continued fraction representations:
  •   [a0;a1,... an−1,an]     and
  •   [a0;a1,... an−1,an−1,1]
  • Choose one of the above that will give a binary expansion ending with a   1.
  • Produce the inverse function ?-1(x)
  • Verify that ?(φ) = 5/3, where φ is the Greek golden ratio.
  • Verify that ?-1(-5/9) = (√13 - 7)/6
  • Verify that the two functions are inverses of each other by showing that ?-1(?(x))=x and ?(?-1(y))=y for x, y of your choice


Don't worry about precision error in the last few digits.

See also



Factor[edit]

USING: formatting kernel make math math.constants
math.continued-fractions math.functions math.parser
math.statistics sequences sequences.extras splitting.monotonic
vectors ;
 
CONSTANT: max-iter 151
 
: >continued-fraction ( x -- seq )
0 swap 1vector
[ dup last integer? pick max-iter > or ]
[ dup next-approx [ 1 + ] dip ] until nip
dup last integer? [ but-last-slice ] unless ;
 
: ? ( x -- y )
>continued-fraction unclip swap cum-sum
[ max-iter < ] take-while
[ even? 1 -1 kernel:? swap 2^ / ] map-index
sum 2 * + >float ;
 
: (float>bin) ( x -- y )
[ dup 0 > ]
[ 2 * dup >integer # dup 1 >= [ 1 - ] when ] while ;
 
: float>bin ( x -- n str )
>float dup >integer [ - ] keep swap abs
[ 0 # (float>bin) ] "" make nip ;
 
: ?⁻¹ ( x -- y )
dup float>bin [ = ] monotonic-split
[ length ] map swap prefix >ratio swap copysign ;
 
: compare ( x y -- ) "%-25u%-25u\n" printf ;
 
phi ? 5 3 /f compare
-5/9 ?⁻¹ 13 sqrt 7 - 6 /f compare
0.718281828 ?⁻¹ ? 0.1213141516171819 ? ?⁻¹ compare
Output:
1.666666666668335        1.666666666666667        
-0.5657414540893351      -0.5657414540893352      
0.718281828000002        0.1213141516171819       

FreeBASIC[edit]

#define MAXITER 151
 
function minkowski( x as double ) as double
if x>1 or x<0 then return int(x)+minkowski(x-int(x))
dim as ulongint p = int(x)
dim as ulongint q = 1, r = p + 1, s = 1, m, n
dim as double d = 1, y = p
while true
d = d / 2.0
if y + d = y then exit while
m = p + r
if m < 0 or p < 0 then exit while
n = q + s
if n < 0 then exit while
if x < cast(double,m) / n then
r = m
s = n
else
y = y + d
p = m
q = n
end if
wend
return y + d
end function
 
function minkowski_inv( byval x as double ) as double
if x>1 or x<0 then return int(x)+minkowski_inv(x-int(x))
if x=1 or x=0 then return x
redim as uinteger contfrac(0 to 0)
dim as uinteger curr=0, count=1, i = 0
do
x *= 2
if curr = 0 then
if x<1 then
count += 1
else
i += 1
redim preserve contfrac(0 to i)
contfrac(i-1)=count
count = 1
curr = 1
x=x-1
endif
else
if x>1 then
count += 1
x=x-1
else
i += 1
redim preserve contfrac(0 to i)
contfrac(i-1)=count
count = 1
curr = 0
endif
end if
if x = int(x) then
contfrac(i)=count
exit do
end if
loop until i = MAXITER
dim as double ret = 1.0/contfrac(i)
for j as integer = i-1 to 0 step -1
ret = contfrac(j) + 1.0/ret
next j
return 1./ret
end function
 
print minkowski( 0.5*(1+sqr(5)) ), 5./3
print minkowski_inv( -5./9 ), (sqr(13)-7)/6
print minkowski(minkowski_inv(0.718281828)), minkowski_inv(minkowski(0.1213141516171819))
 
Output:
 1.666666666669698           1.666666666666667
-0.5657414540893351         -0.5657414540893352
 0.7182818280000092          0.1213141516171819

Go[edit]

Translation of: FreeBASIC
package main
 
import (
"fmt"
"math"
)
 
const MAXITER = 151
 
func minkowski(x float64) float64 {
if x > 1 || x < 0 {
return math.Floor(x) + minkowski(x-math.Floor(x))
}
p := uint64(x)
q := uint64(1)
r := p + 1
s := uint64(1)
d := 1.0
y := float64(p)
for {
d = d / 2
if y+d == y {
break
}
m := p + r
if m < 0 || p < 0 {
break
}
n := q + s
if n < 0 {
break
}
if x < float64(m)/float64(n) {
r = m
s = n
} else {
y = y + d
p = m
q = n
}
}
return y + d
}
 
func minkowskiInv(x float64) float64 {
if x > 1 || x < 0 {
return math.Floor(x) + minkowskiInv(x-math.Floor(x))
}
if x == 1 || x == 0 {
return x
}
contFrac := []uint32{0}
curr := uint32(0)
count := uint32(1)
i := 0
for {
x *= 2
if curr == 0 {
if x < 1 {
count++
} else {
i++
t := contFrac
contFrac = make([]uint32, i+1)
copy(contFrac, t)
contFrac[i-1] = count
count = 1
curr = 1
x--
}
} else {
if x > 1 {
count++
x--
} else {
i++
t := contFrac
contFrac = make([]uint32, i+1)
copy(contFrac, t)
contFrac[i-1] = count
count = 1
curr = 0
}
}
if x == math.Floor(x) {
contFrac[i] = count
break
}
if i == MAXITER {
break
}
}
ret := 1.0 / float64(contFrac[i])
for j := i - 1; j >= 0; j-- {
ret = float64(contFrac[j]) + 1.0/ret
}
return 1.0 / ret
}
 
func main() {
fmt.Printf("%19.16f %19.16f\n", minkowski(0.5*(1+math.Sqrt(5))), 5.0/3.0)
fmt.Printf("%19.16f %19.16f\n", minkowskiInv(-5.0/9.0), (math.Sqrt(13)-7)/6)
fmt.Printf("%19.16f %19.16f\n", minkowski(minkowskiInv(0.718281828)),
minkowskiInv(minkowski(0.1213141516171819)))
}
Output:
 1.6666666666696983  1.6666666666666667
-0.5657414540893351 -0.5657414540893352
 0.7182818280000092  0.1213141516171819

Julia[edit]

Translation of: FreeBASIC
function minkowski(x)
p = Int(floor(x))
(x > 1 || x < 0) && return p + minkowski(x)
q, r, s, m, n = 1, p + 1, 1, 0, 0
d, y = 1.0, Float64(p)
while true
d /= 2.0
y + d == y && break
m = p + r
(m < 0 || p < 0) && break
n = q + s
n < 0 && break
if x < (m / n)
r, s = m, n
else
y, p, q = y + d, m, n
end
end
return y + d
end
 
function minkowski_inv(x, maxiter=151)
p = Int(floor(x))
(x > 1 || x < 0) && return p + minkowski_inv(x - p, maxiter)
(x == 1 || x == 0) && return x
contfrac = [0]
curr, coun, i = 0, 1, 0
while i < maxiter
x *= 2
if curr == 0
if x < 1
coun += 1
else
i += 1
push!(contfrac, 0)
contfrac[i] = coun
coun = 1
curr = 1
x -= 1
end
else
if x > 1
coun += 1
x -= 1
else
i += 1
push!(contfrac, 0)
contfrac[i] = coun
coun = 1
curr = 0
end
end
if x == Int(floor(x))
contfrac[i + 1] = coun
break
end
end
ret = 1.0 / contfrac[i + 1]
for j in i:-1:1
ret = contfrac[j] + 1.0 / ret
end
return 1.0 / ret
end
 
println(" ", minkowski((1 + sqrt(5)) / 2), " ", 5 / 3)
println(minkowski_inv(-5/9), " ", (sqrt(13) - 7) / 6)
println(" ", minkowski(minkowski_inv(0.718281828)), " ",
minkowski_inv(minkowski(0.1213141516171819)))
 
Output:
 1.6666666666696983   1.6666666666666667
-0.5657414540893351   -0.5657414540893352
 0.7182818280000092   0.12131415161718191

Phix[edit]

Translation of: FreeBASIC
constant MAXITER = 151
 
function minkowski(atom x)
atom p = floor(x)
if x>1 or x<0 then return p+minkowski(x-p) end if
atom q = 1, r = p + 1, s = 1, m, n, d = 1, y = p
while true do
d = d/2
if y + d = y then exit end if
m = p + r
if m < 0 or p < 0 then exit end if
n = q + s
if n < 0 then exit end if
if x < m/n then
r = m
s = n
else
y = y + d
p = m
q = n
end if
end while
return y + d
end function
 
function minkowski_inv(atom x)
if x>1 or x<0 then return floor(x)+minkowski_inv(x-floor(x)) end if
if x=1 or x=0 then return x end if
sequence contfrac = {}
integer curr = 0, count = 1
while true do
x *= 2
if curr = 0 then
if x<1 then
count += 1
else
contfrac &= count
count = 1
curr = 1
x -= 1
end if
else
if x>1 then
count += 1
x -= 1
else
contfrac &= count
count = 1
curr = 0
end if
end if
if x = floor(x) then
contfrac &= count
exit
end if
if length(contfrac)=MAXITER then exit end if
end while
atom ret = 1/contfrac[$]
for i = length(contfrac)-1 to 1 by -1 do
ret = contfrac[i] + 1.0/ret
end for
return 1/ret
end function
 
printf(1,"%20.16f %20.16f\n",{minkowski(0.5*(1+sqrt(5))), 5/3})
printf(1,"%20.16f %20.16f\n",{minkowski_inv(-5/9), (sqrt(13)-7)/6})
printf(1,"%20.16f %20.16f\n",{minkowski(minkowski_inv(0.718281828)),
minkowski_inv(minkowski(0.1213141516171819))})
Output:
  1.6666666666696983   1.6666666666666668
 -0.5657414540893351  -0.5657414540893352
  0.7182818280000092   0.1213141516171819

Raku[edit]

Translation of: Go
# 20201120 Raku programming solution
 
my \MAXITER = 151;
 
sub minkowski(\x) {
 
return x.floor + minkowski( x - x.floor ) if x > 1 || x < 0 ;
 
my $y = my $p = x.floor;
my ($q,$s,$d) = 1 xx 3;
my $r = $p + 1;
 
loop {
last if ( $y + ($d /= 2) == $y ) ||
( my $m = $p + $r) < 0 | $p < 0 ||
( my $n = $q + $s) < 0 ;
x < $m/$n ?? ( ($r,$s) = ($m, $n) ) !! ( $y += $d; ($p,$q) = ($m, $n) );
}
return $y + $d
}
 
sub minkowskiInv($x is copy) {
 
return $x.floor + minkowskiInv($x - $x.floor) if $x > 1 || $x < 0 ;
 
return $x if $x == 1 || $x == 0 ;
 
my @contFrac = 0;
my $i = my $curr = 0 ; my $count = 1;
 
loop {
$x *= 2;
if $curr == 0 {
if $x < 1 {
$count++
} else {
$i++;
@contFrac.append: 0;
@contFrac[$i-1] = $count;
($count,$curr) = 1,1;
$x--;
}
} else {
if $x > 1 {
$count++;
$x--;
} else {
$i++;
@contFrac.append: 0;
@contFrac[$i-1] = $count;
($count,$curr) = 1,0;
}
}
if $x == $x.floor { @contFrac[$i] = $count ; last }
last if $i == MAXITER;
}
my $ret = 1 / @contFrac[$i];
loop (my $j = $i - 1; $j0; $j--) { $ret = @contFrac[$j] + 1/$ret }
return 1 / $ret
}
 
printf "%19.16f %19.16f\n", minkowski(0.5*(1 + 5.sqrt)), 5/3;
printf "%19.16f %19.16f\n", minkowskiInv(-5/9), (13.sqrt-7)/6;
printf "%19.16f %19.16f\n", minkowski(minkowskiInv(0.718281828)),
minkowskiInv(minkowski(0.1213141516171819))
Output:
 1.6666666666696983  1.6666666666666667
-0.5657414540893351 -0.5657414540893352
 0.7182818280000092  0.1213141516171819

REXX[edit]

Translation of: FreeBASIC
Translation of: Phix
/*REXX program uses the Minkowski question─mark function to convert a continued fraction*/
numeric digits 20 /*use enough dec. digits for precision.*/
say fmt(mink(0.5*(1+sqrt(5)))) fmt(5/3 )
say fmt(minkI(-5/9)) fmt((sqrt(13)-7)/6)
say fmt(mink(minkI(0.718281828))) fmt(mink(minkI(.1213141516171819)))
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
floor: procedure; parse arg x; _= trunc(x); return _ - (x<0) * (x\=_)
fmt: procedure: parse arg z; return right( format(z, , digits() - 2, 0), digits() +5)
/*──────────────────────────────────────────────────────────────────────────────────────*/
mink: procedure: parse arg x; p= x % 1; if x>1 | x<0 then return p + mink(x-p)
q= 1; s= 1; m= 0; n= 0; d= 1; y= p
r= p + 1
do forever; d= d * 0.5; if y+d=y | d=0 then leave /*d= d÷2*/
m= p + r; if m<0 | p<0 then leave
n= q + s; if n<0 then leave
if x<m/n then do; r= m; s= n; end
else do; y= y + d; p= m; q= n; end
end /*forever*/
return y + d
/*──────────────────────────────────────────────────────────────────────────────────────*/
minkI: procedure; parse arg x; p= floor(x); if x>1 | x<0 then return p + minkI(x-p)
if x=1 | x=0 then return x
curr= 0; count= 1; maxIter= 200; $=
 
do until count==maxIter | words($)==maxIter; x= x + x /*a fast double*/
if curr==0 then if x<1 then count= count + 1
else do; $= $ count; count= 1; curr= 1; x= x-1; end
else if x>1 then do; count= count + 1; x= x-1; end
else do; $= $ count; count= 1; curr= 0; end
if x==floor(x) then do; $= $ count; leave; end
end /*until*/
 
#= words($)
ret= 1 / word($, #)
do j=# for # by -1; ret= word($, j) + 1 / ret
end /*j*/
return 1 / ret
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6
numeric form; m.=9; parse value format(x,2,1,,0) 'E0' with g "E" _ .; g=g *.5'e'_ %2
do j=0 while h>9; m.j= h; h= h % 2 + 1; end /*j*/
do k=j+5 to 0 by -1; numeric digits m.k; g= (g + x/g) * .5; end /*k*/; return g
output   when using the internal default inputs:
     1.666666666666666963      1.666666666666666667
    -0.565741454089335118     -0.565741454089335118
     0.718281828000000011      0.121314151617181900

Wren[edit]

Translation of: FreeBASIC
Library: Wren-fmt
import "/fmt" for Fmt
 
var MAXITER = 151
 
var minkowski // predeclare as recursive
minkowski = Fn.new { |x|
if (x > 1 || x < 0) return x.floor + minkowski.call(x - x.floor)
var p = x.floor
var q = 1
var r = p + 1
var s = 1
var d = 1
var y = p
while (true) {
d = d / 2
if (y + d == y) break
var m = p + r
if (m < 0 || p < 0 ) break
var n = q + s
if (n < 0) break
if (x < m/n) {
r = m
s = n
} else {
y = y + d
p = m
q = n
}
}
return y + d
}
 
var minkowskiInv
minkowskiInv = Fn.new { |x|
if (x > 1 || x < 0) return x.floor + minkowskiInv.call(x - x.floor)
if (x == 1 || x == 0) return x
var contFrac = [0]
var curr = 0
var count = 1
var i = 0
while (true) {
x = x * 2
if (curr == 0) {
if (x < 1) {
count = count + 1
} else {
i = i + 1
var t = contFrac
contFrac = List.filled(i + 1, 0)
for (j in 0...t.count) contFrac[j] = t[j]
contFrac[i-1] = count
count = 1
curr = 1
x = x - 1
}
} else {
if (x > 1) {
count = count + 1
x = x - 1
} else {
i = i + 1
var t = contFrac
contFrac = List.filled(i + 1, 0)
for (j in 0...t.count) contFrac[j] = t[j]
contFrac[i-1] = count
count = 1
curr = 0
}
}
if (x == x.floor) {
contFrac[i] = count
break
}
if (i == MAXITER) break
}
var ret = 1/contFrac[i]
for (j in i-1..0) ret = contFrac[j] + 1/ret
return 1/ret
}
 
Fmt.print("$17.16f $17.14f", minkowski.call(0.5 * (1 + 5.sqrt)), 5/3)
Fmt.print("$17.14f $17.14f", minkowskiInv.call(-5/9), (13.sqrt - 7)/6)
Fmt.print("$17.14f $17.14f", minkowski.call(minkowskiInv.call(0.718281828)),
minkowskiInv.call(minkowski.call(0.1213141516171819)))
Output:
 1.66666666666970  1.66666666666667
-0.56574145408934 -0.56574145408934
 0.71828182800001  0.12131415161718