Minkowski question-mark function

From Rosetta Code
Task
Minkowski question-mark function
You are encouraged to solve this task according to the task description, using any language you may know.

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



11l[edit]

Translation of: Python
-V MAXITER = 151

F minkowski(x) -> Float
   I x > 1 | x < 0
      R floor(x) + minkowski(x - floor(x))

   V p = Int(x)
   V q = 1
   V r = p + 1
   V s = 1
   V d = 1.0
   V y = Float(p)

   L
      d /= 2
      I y + d == y
         L.break

      V m = p + r
      I m < 0 | p < 0
         L.break

      V n = q + s
      I n < 0
         L.break

      I x < Float(m) / n
         r = m
         s = n
      E
         y += d
         p = m
         q = n

   R y + d

F minkowski_inv(=x) -> Float
   I x > 1 | x < 0
      R floor(x) + minkowski_inv(x - floor(x))

   I x == 1 | x == 0
      R x

   V cont_frac = [0]
   V current = 0
   V count = 1
   V i = 0

   L
      x *= 2

      I current == 0
         I x < 1
            count++
         E
            cont_frac.append(0)
            cont_frac[i] = count

            i++
            count = 1
            current = 1
            x--
      E
         I x > 1
            count++
            x--
         E
            cont_frac.append(0)
            cont_frac[i] = count

            i++
            count = 1
            current = 0

      I x == floor(x)
         cont_frac[i] = count
         L.break

      I i == :MAXITER
         L.break

   V ret = 1.0 / cont_frac[i]
   L(j) (i-1 .. 0).step(-1)
      ret = cont_frac[j] + 1.0 / ret

   R 1.0 / ret

print(‘#2.16 #2.16’.format(minkowski(0.5 * (1 + sqrt(5))), 5.0 / 3.0))
print(‘#2.16 #2.16’.format(minkowski_inv(-5.0 / 9.0), (sqrt(13) - 7) / 6))
print(‘#2.16 #2.16’.format(minkowski(minkowski_inv(0.718281828)), minkowski_inv(minkowski(0.1213141516171819))))
Output:
 1.6666666666696983  1.6666666666666667
-0.5657414540893350 -0.5657414540893351
 0.7182818280000091  0.1213141516171819

Factor[edit]

Works with: Factor version 0.99 2020-08-14
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

Haskell[edit]

In a lazy functional language Minkowski question mark function can be implemented using one of it's basic properties:

?(p+r)/(q+s) = 1/2 * ( ?(p/q) + ?(r/s) ), ?(0) = 0, ?(1) = 1.

where p/q and r/s are fractions, such that |ps - rq| = 1.

This recursive definition can be implemented as lazy corecursion, i.e. by generating two infinite binary trees: mediant-based Stern-Brocot tree, containing all rationals, and mean-based tree with corresponding values of Minkowsky ?-function. There is one-to-one correspondence between these two trees so both ?(x) and ?-1(x) may be implemented as mapping between them. For details see the paper [[1]] (in Russian).

import Data.Tree
import Data.Ratio
import Data.List

intervalTree :: (a -> a -> a) -> (a, a) -> Tree a
intervalTree node = unfoldTree $
  \(a, b) -> let m = node a b in (m, [(a,m), (m,b)])

Node a _ ==> Node b [] = const b
Node a [] ==> Node b _ = const b
Node a [l1, r1] ==> Node b [l2, r2] =
  \x -> case x `compare` a of
          LT -> (l1 ==> l2) x
          EQ -> b
          GT -> (r1 ==> r2) x

mirror :: Num a => Tree a -> Tree a
mirror t = Node 0 [reflect (negate <$> t), t]
  where
    reflect (Node a [l,r]) = Node a [reflect r, reflect l]

------------------------------------------------------------

sternBrocot :: Tree Rational
sternBrocot = toRatio <$> intervalTree mediant ((0,1), (1,0))
  where
    mediant (p, q) (r, s) = (p + r, q + s)

toRatio (p, q) = p % q

minkowski :: Tree Rational
minkowski = toRatio <$> intervalTree mean ((0,1), (1,0))

mean (p, q) (1, 0) = (p+1, q)
mean (p, q) (r, s) = (p*s + q*r, 2*q*s)


questionMark, invQuestionMark :: Rational -> Rational
questionMark    = mirror sternBrocot ==> mirror minkowski
invQuestionMark = mirror minkowski ==> mirror sternBrocot

------------------------------------------------------------
-- Floating point trees and functions

sternBrocotF :: Tree Double
sternBrocotF = mirror $ fromRational <$> sternBrocot

minkowskiF :: Tree Double
minkowskiF = mirror $ intervalTree mean (0, 1/0)
  where
    mean a b | isInfinite b = a + 1
             | otherwise = (a + b) / 2

questionMarkF, invQuestionMarkF :: Double -> Double
questionMarkF = sternBrocotF ==> minkowskiF
invQuestionMarkF = minkowskiF ==> sternBrocotF
λ> mapM_ print $ take 4 $ levels farey
[1 % 2]
[1 % 3,2 % 3]
[1 % 4,2 % 5,3 % 5,3 % 4]
[1 % 5,2 % 7,3 % 8,3 % 7,4 % 7,5 % 8,5 % 7,4 % 5]

λ> mapM_ print $ take 4 $ levels minkowski
[1 % 2]
[1 % 4,3 % 4]
[1 % 8,3 % 8,5 % 8,7 % 8]
[1 % 16,3 % 16,5 % 16,7 % 16,9 % 16,11 % 16,13 % 16,15 % 16]

λ> questionMark (1/2) 1 % 2 λ> questionMark (2/7) 3 % 16 λ> questionMark (-22/7) (-193) % 64 λ> invQuestionMark (3/16) 2 % 7 λ> invQuestionMark (13/256)

5 % 27
λ> questionMark $ (sqrt 5 + 1) / 2
1.6666666666678793
λ> 5/3
1.6666666666666667
λ> invQuestionMark (-5/9)
-0.5657414540893351
λ> (sqrt 13 - 7)/6
-0.5657414540893352

J[edit]

Implementation:

ITERCOUNT=: 52

minkowski=: {{
  f=. 1|y
  node=. *i.2 2 NB. node of Stern-Brocot tree
  B=. ''
  for. i.ITERCOUNT do.
    B=. B, b=. f>:%/t=. +/node
    node=. t (1-b)} node 
  end.
  (<.y)+B+/ .*2^-1+i.ITERCOUNT
}}
 
invmink=: {{
  f=. 1|y
  cf=. i.0
  cur=. 0 NB. 1 if generating "top" side of cf
  cnt=. 1 NB. proposed continued fraction element
  for. i.ITERCOUNT do.
    if. f=<. f do.
      cf=. cf,%cnt break.
    end.
    f=. f*2
    b=. 1 >`<@.cur f
    cf=. cf,(-.b)#cnt
    cnt=. 1+b*cnt
    cur=. cur=b
    f=. f-cur
  end.
  (+%)/(<.y),cf
}}

That said, note that this algorithm introduces significant numeric instability for √7 divided by 3:

   (minkowski@invmink - invmink@minkowski) (p:%%:)3
1.10713e_6

I see this same instability using the python implementation and appending:

    print(
        "{:19.16f} {:19.16f}".format(
            minkowski(minkowski_inv(4.04145188432738056)),
            minkowski_inv(minkowski(4.04145188432738056)),
        )
    )

Using an exact fraction for 4.04145188432738056 and bumping the iteration count from 52 up to 200 changes that difference to 1.43622e_12.

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

Mathematica / Wolfram Language[edit]

ClearAll[InverseMinkowskiQuestionMark]
InverseMinkowskiQuestionMark[val_] := Module[{x}, (x /. FindRoot[MinkowskiQuestionMark[x] == val, {x, Floor[val], Ceiling[val]}])]
MinkowskiQuestionMark[GoldenRatio]
InverseMinkowskiQuestionMark[-5/9] // RootApproximant
MinkowskiQuestionMark[InverseMinkowskiQuestionMark[0.1213141516171819]]
InverseMinkowskiQuestionMark[MinkowskiQuestionMark[0.1213141516171819]]
Output:
5/3
1/6 (-7+Sqrt[13])
0.121314
0.121314

Nim[edit]

Translation of: Go
import math, strformat

const MaxIter = 151


func minkowski(x: float): float =

  if x notin 0.0..1.0:
    return floor(x) + minkowski(x - floor(x))

  var
    p = x.uint64
    r = p + 1
    q, s = 1u64
    d = 1.0
    y = p.float

  while true:
    d /= 2
    if y + d == y: break
    let m = p + r
    if m < 0 or p < 0: break
    let n = q + s
    if n < 0: break
    if x < m.float / n.float:
      r = m
      s = n
    else:
      y += d
      p = m
      q = n

  result = y + d


func minkowskiInv(x: float): float =

  if x notin 0.0..1.0:
    return floor(x) + minkowskiInv(x - floor(x))
  if x == 1 or x == 0:
    return x

  var
    contFrac: seq[uint32]
    curr = 0u32
    count = 1u32
    i = 0
    x = x

  while true:
    x *= 2
    if curr == 0:
      if x < 1:
        inc count
      else:
        inc i
        contFrac.setLen(i + 1)
        contFrac[i - 1] = count
        count = 1
        curr = 1
        x -= 1
    else:
      if x > 1:
        inc count
        x -= 1
      else:
        inc i
        contFrac.setLen(i + 1)
        contFrac[i - 1] = count
        count = 1
        curr = 0
    if x == floor(x):
      contFrac[i] = count
      break
    if i == MaxIter:
      break

  var ret = 1 / contFrac[i].float
  for j in countdown(i - 1, 0):
    ret = contFrac[j].float + 1 / ret
  result = 1 / ret


echo &"{minkowski(0.5*(1+sqrt(5.0))):19.16f}, {5/3:19.16f}"
echo &"{minkowskiInv(-5/9):19.16f}, {(sqrt(13.0)-7)/6:19.16f}"
echo &"{minkowski(minkowskiInv(0.718281828)):19.16f}, " &
     &"{minkowskiInv(minkowski(0.1213141516171819)):19.16f}"
Output:
 1.6666666666696983,  1.6666666666666667
-0.5657414540893351, -0.5657414540893352
 0.7182818280000092,  0.1213141516171819

Perl[edit]

Translation of: Raku
use strict;
use warnings;
use feature 'say';
use POSIX qw(floor);

my $MAXITER = 50;

sub minkowski {
    my($x) = @_;

    return floor($x) + minkowski( $x - floor($x) ) if $x > 1 || $x < 0 ;

    my $y = my $p = floor($x);
    my ($q,$s,$d) = (1,1,1);
    my $r = $p + 1;

    while () {
        last if ( $y + ($d /= 2)  == $y ) or
                ( my $m = $p + $r) <  0   or
                ( my $n = $q + $s) <  0;
        $x < $m/$n ? ($r,$s) = ($m, $n) : ($y += $d and ($p,$q) = ($m, $n) );
    }
    return $y + $d
}

sub minkowskiInv {
    my($x) = @_;

    return floor($x) + minkowskiInv($x - floor($x)) if $x > 1 || $x < 0;
    return $x if $x == 1 || $x == 0 ;

    my @contFrac = 0;
    my $i = my $curr = 0 ; my $count = 1;

    while () {
        $x *= 2;
        if ($curr == 0) {
            if ($x < 1) {
                $count++
            } else {
                $i++;
                push @contFrac, 0;
                $contFrac[$i-1] = $count;
                ($count,$curr) = (1,1);
                $x--;
            }
        } else {
            if ($x > 1) {
                $count++;
                $x--;
            } else {
                $i++;
                push @contFrac, 0;
                @contFrac[$i-1] = $count;
                ($count,$curr) = (1,0);
            }
        }
        if ($x == floor($x)) { @contFrac[$i] = $count; last }
        last if $i == $MAXITER;
    }
    my $ret = 1 / $contFrac[$i];
    for (my $j = $i - 1; $j >= 0; $j--) { $ret = $contFrac[$j] + 1/$ret }
    return 1 / $ret
}

printf "%19.16f %19.16f\n", minkowski(0.5*(1 + sqrt(5))), 5/3;
printf "%19.16f %19.16f\n", minkowskiInv(-5/9), (sqrt(13)-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

Phix[edit]

Translation of: FreeBASIC
with javascript_semantics
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

Python[edit]

Translation of: Go
import math

MAXITER = 151


def minkowski(x):
    if x > 1 or x < 0:
        return math.floor(x) + minkowski(x - math.floor(x))

    p = int(x)
    q = 1
    r = p + 1
    s = 1
    d = 1.0
    y = float(p)

    while True:
        d /= 2
        if y + d == y:
            break

        m = p + r
        if m < 0 or p < 0:
            break

        n = q + s
        if n < 0:
            break

        if x < m / n:
            r = m
            s = n
        else:
            y += d
            p = m
            q = n

    return y + d


def minkowski_inv(x):
    if x > 1 or x < 0:
        return math.floor(x) + minkowski_inv(x - math.floor(x))

    if x == 1 or x == 0:
        return x

    cont_frac = [0]
    current = 0
    count = 1
    i = 0

    while True:
        x *= 2

        if current == 0:
            if x < 1:
                count += 1
            else:
                cont_frac.append(0)
                cont_frac[i] = count

                i += 1
                count = 1
                current = 1
                x -= 1
        else:
            if x > 1:
                count += 1
                x -= 1
            else:
                cont_frac.append(0)
                cont_frac[i] = count

                i += 1
                count = 1
                current = 0

        if x == math.floor(x):
            cont_frac[i] = count
            break

        if i == MAXITER:
            break

    ret = 1.0 / cont_frac[i]
    for j in range(i - 1, -1, -1):
        ret = cont_frac[j] + 1.0 / ret

    return 1.0 / ret


if __name__ == "__main__":
    print(
        "{:19.16f} {:19.16f}".format(
            minkowski(0.5 * (1 + math.sqrt(5))),
            5.0 / 3.0,
        )
    )

    print(
        "{:19.16f} {:19.16f}".format(
            minkowski_inv(-5.0 / 9.0),
            (math.sqrt(13) - 7) / 6,
        )
    )

    print(
        "{:19.16f} {:19.16f}".format(
            minkowski(minkowski_inv(0.718281828)),
            minkowski_inv(minkowski(0.1213141516171819)),
        )
    )
Output:
 1.6666666666696983  1.6666666666666667
-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 40                                /*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;    t= trunc(x);     return t   -   (x<0)  *  (x\=t)
fmt:   procedure: parse arg a;    d= digits();     return right( format(a, , d-2, 0), d+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
       cur= 0;                limit= 200;       $=               /*limit: max iterations*/
       #= 1                                                      /*#:  is the count.    */
                  do  until #==limit | words($)==limit;                        x= x * 2
                  if cur==0  then if x<1  then      #= # + 1
                                          else do;  $= $ #;  #= 1;   cur= 1;  x= x-1;  end
                             else if x>1  then do;           #= # + 1;        x= x-1;  end
                                          else do;  $= $ #;  #= 1;   cur= 0;           end
                  if x==floor(x)          then do;           $= $ #;  leave;           end
                  end   /*until*/
       z= words($)
       ret= 1 / word($, z)
                               do j=z  for z  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.66666666666666666666666666673007566392      1.66666666666666666666666666666666666667
    -0.56574145408933511781346312208825067563     -0.56574145408933511781346312208825067562
     0.71828182799999999999999999999999992890      0.12131415161718190000000000000000000833

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