Minkowski question-mark function

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 [a₀; a₁, a₂, a₃, ...] of a number into a binary decimal representation in which the integer part a₀ is unchanged and the a₁, a₂, ... 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.0111₂.

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:

•   [a₀;a₁,... a_n−1,a_n]     and
•   [a₀;a₁,... a_n−1,a_n−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

• Wikipedia entry: Minkowski's question-mark function

11l

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

C#

Translation of: Go

using System;

class Program
{
    const int MAXITER = 151;

    static double Minkowski(double x)
    {
        if (x > 1 || x < 0)
        {
            return Math.Floor(x) + Minkowski(x - Math.Floor(x));
        }
        ulong p = (ulong)x;
        ulong q = 1;
        ulong r = p + 1;
        ulong s = 1;
        double d = 1.0;
        double y = (double)p;
        while (true)
        {
            d = d / 2;
            if (y + d == y)
            {
                break;
            }
            ulong m = p + r;
            if (m < 0 || p < 0)
            {
                break;
            }
            ulong n = q + s;
            if (n < 0)
            {
                break;
            }
            if (x < (double)m / (double)n)
            {
                r = m;
                s = n;
            }
            else
            {
                y = y + d;
                p = m;
                q = n;
            }
        }
        return y + d;
    }

    static double MinkowskiInv(double x)
    {
        if (x > 1 || x < 0)
        {
            return Math.Floor(x) + MinkowskiInv(x - Math.Floor(x));
        }
        if (x == 1 || x == 0)
        {
            return x;
        }
        uint[] contFrac = new uint[] { 0 };
        uint curr = 0;
        uint count = 1;
        int i = 0;
        while (true)
        {
            x *= 2;
            if (curr == 0)
            {
                if (x < 1)
                {
                    count++;
                }
                else
                {
                    i++;
                    Array.Resize(ref contFrac, i + 1);
                    contFrac[i - 1] = count;
                    count = 1;
                    curr = 1;
                    x--;
                }
            }
            else
            {
                if (x > 1)
                {
                    count++;
                    x--;
                }
                else
                {
                    i++;
                    Array.Resize(ref contFrac, i + 1);
                    contFrac[i - 1] = count;
                    count = 1;
                    curr = 0;
                }
            }
            if (x == Math.Floor(x))
            {
                contFrac[i] = count;
                break;
            }
            if (i == MAXITER)
            {
                break;
            }
        }
        double ret = 1.0 / contFrac[i];
        for (int j = i - 1; j >= 0; j--)
        {
            ret = contFrac[j] + 1.0 / ret;
        }
        return 1.0 / ret;
    }

    static void Main(string[] args)
    {
        Console.WriteLine("{0,19:0.0000000000000000} {1,19:0.0000000000000000}", Minkowski(0.5 * (1 + Math.Sqrt(5))), 5.0 / 3.0);
        Console.WriteLine("{0,19:0.0000000000000000} {1,19:0.0000000000000000}", MinkowskiInv(-5.0 / 9.0), (Math.Sqrt(13) - 7) / 6);
        Console.WriteLine("{0,19:0.0000000000000000} {1,19:0.0000000000000000}", Minkowski(MinkowskiInv(0.718281828)),
            MinkowskiInv(Minkowski(0.1213141516171819)));
    }
}

Output:

 1.6666666666697000  1.6666666666666700
-0.5657414540893350 -0.5657414540893350
 0.7182818280000090  0.1213141516171820

C++

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>

constexpr int32_t MAX_ITERATIONS = 151;

double minkowski(const double& x) {
	if ( x < 0 || x > 1 ) {
		return floor(x) + minkowski(x - floor(x));
	}

	int64_t p = (int64_t) x;
	int64_t q = 1;
	int64_t r = p + 1;
	int64_t s = 1;
	double d = 1.0;
	double y = (double) p;

	while ( true ) {
		d /= 2;
		if ( d == 0.0 ) {
			break;
		}

		int64_t m = p + r;
		if ( m < 0 || p < 0 ) {
			break;
		}

		int64_t n = q + s;
		if ( n < 0 ) {
			break;
		}

		if ( x < (double) m / n ) {
			r = m;
			s = n;
		} else {
			y += d;
			p = m;
			q = n;
		}
	}
	return y + d;
}

double minkowski_inverse(double x) {
	if ( x < 0 || x > 1 ) {
		return floor(x) + minkowski_inverse(x - floor(x));
	}

	if ( x == 0 || x == 1 ) {
		return x;
	}

	std::vector<int32_t> continued_fraction(1, 0);
	int32_t current = 0;
	int32_t count = 1;
	int32_t i = 0;

	while ( true ) {
		x *= 2;
		if ( current == 0 ) {
			if ( x < 1 ) {
				count += 1;
			} else {
				 continued_fraction.emplace_back(0);
				 continued_fraction[i] = count;

				 i += 1;
				 count = 1;
				 current = 1;
				 x -= 1;
			}
		} else {
			 if ( x > 1 ) {
				 count += 1;
				 x -= 1;
			 } else {
				 continued_fraction.emplace_back(0);
				 continued_fraction[i] = count;

				 i += 1;
				 count = 1;
				 current = 0;
			 }
		}

		if ( x == floor(x) ) {
			 continued_fraction[i] = count;
			 break;
		}

		if ( i == MAX_ITERATIONS ) {
			 break;
		}
	}

	double reciprocal = 1.0 / continued_fraction[i];
	for ( int32_t j = i - 1; j >= 0; --j ) {
		 reciprocal = continued_fraction[j] + 1.0 / reciprocal;
	}

	return 1.0 / reciprocal;
}

int main() {
	std::cout << std::setw(20) << std::fixed << std::setprecision(16) << minkowski(0.5 * ( 1 + sqrt(5) ))
			  << std::setw(20) << 5.0 / 3.0 << std::endl;
	std::cout << std::setw(20) << minkowski_inverse(-5.0 / 9.0)
			  << std::setw(20) << ( sqrt(13) - 7 ) / 6  << std::endl;
	std::cout << std::setw(20) << minkowski(minkowski_inverse(0.718281828182818))
			  << std::setw(20) << 0.718281828182818 << std::endl;
	std::cout << std::setw(20) << minkowski_inverse(minkowski(0.1213141516271819))
			  << std::setw(20) << 0.1213141516171819 << std::endl;
}

Output:

  1.6666666666696983  1.6666666666666667
 -0.5657414540893351 -0.5657414540893352
  0.7182818281828269  0.7182818281828180
  0.1213141516171819  0.1213141516171819

F#

// Minkowski question-mark function. Nigel Galloway: July 14th., 2023
let fN g=let n=(int>>float)g in ((if g<0.0 then -1.0 else 1.0),abs n,abs (g-n))
let fI(n,_,(nl,nh))(g,_,(gl,gh))=let l,h=nl+gl,nh+gh in ((n+g)/2.0,(float l)/(float h),(l,h))
let fG n g=(max n g)-(min n g)
let fE(s,z,l)=Seq.unfold(fun(i,e)->let (n,g,_) as r=fI i e in Some((s*(z+n),s*(g+z)),if l<n then (i,r) else if l=n then (r,r) else (r,e)))((0.0,0.0,(0,1)),(1.0,1.0,(1,1)))
let fL(s,z,l)=Seq.unfold(fun(i,e)->let (n,g,_) as r=fI i e in Some((s*(z+n),s*(g+z)),if l<g then (i,r) else if l=g then (r,r) else (r,e)))((0.0,0.0,(0,1)),(1.0,1.0,(1,1)))
let f2M g=let _,(n,_)=fL(fN g)|>Seq.pairwise|>Seq.find(fun((n,_),(g,_))->(fG n g)<2.328306437e-11) in n
let m2F g=let _,(_,n)=fE(fN g)|>Seq.pairwise|>Seq.find(fun((_,n),(_,g))->(fG n g)<2.328306437e-11) in n

printfn $"?(φ) = 5/3 is %A{fG(f2M 1.61803398874989490253)(5.0/3.0)<2.328306437e-10}"
printfn $"?⁻¹(-5/9) = (√13-7)/6 is %A{fG(m2F(-5.0/9.0))((sqrt(13.0)-7.0)/6.0)<2.328306437e-10}"
let n=42.0/23.0 in printfn $"?⁻¹(?(n)) = n is %A{(fG(m2F(f2M n)) n)<2.328306437e-10}"
let n= -3.0/13.0 in printfn $"?(?⁻¹(n)) = n is %A{(fG(f2M(m2F n)) n)<2.328306437e-10}"

Output:

?(φ) = 5/3 is true
?⁻¹(-5/9) = (√13-7)/6 is true
?⁻¹(?(n)) = n is true
?(?⁻¹(n)) = n is true

Factor

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

#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

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

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

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.

Java

import java.util.ArrayList;
import java.util.List;

public final class MinkowskiQuestionMarkFunction {

	public static void main(String[] aArgs) {
		System.out.println(String.format("%20.16f%20.16f", minkowski(0.5 * ( 1 + Math.sqrt(5) )), 5.0 / 3.0));
		System.out.println(String.format("%20.16f%20.16f", minkowskiInverse(-5.0 / 9.0), ( Math.sqrt(13) - 7 ) / 6 ));
		System.out.println(String.format("%20.16f%20.16f", minkowski(minkowskiInverse(0.718281828)), 0.718281828));
		System.out.println(String.format("%20.16f%20.16f",
			minkowskiInverse(minkowski(0.1213141516271819)), 0.1213141516171819));
	}
	
	private static double minkowski(double aX) {
	    if ( aX < 0 || aX > 1 ) {
	        return Math.floor(aX) + minkowski(aX - Math.floor(aX));
	    }
	    
	    long p = (long) aX;
	    long q = 1;
	    long r = p + 1;
	    long s = 1;
	    double d = 1.0;
	    double y = (double) p;
	    
	    while ( true ) {
	        d /= 2;
	        if ( d == 0.0 ) {
	            break;
	        }
	        
	        long m = p + r;
	        if ( m < 0 || p < 0 ) {
	            break;
	        }
	        
	        long n = q + s;
	        if ( n < 0 ) {
	            break;
	        }
	        
	        if ( aX < (double) m / n ) {
	            r = m;
	            s = n;
	        } else {
	            y += d;
	            p = m;
	            q = n;
	        }
	    }	    
	    return y + d;
	}
	
	private static double minkowskiInverse(double aX) {
	    if ( aX < 0 || aX > 1 ) {
	        return Math.floor(aX) + minkowskiInverse(aX - Math.floor(aX));
	    }
	    
	    if ( aX == 0 || aX == 1 ) {
	        return aX;
	    }
	    
	    List<Integer> continuedFraction = new ArrayList<Integer>();
	    continuedFraction.add(0);
	    int current = 0;
	    int count = 1;
	    int i = 0;
	    
	    while ( true ) {
	        aX *= 2;
	        if ( current == 0 ) {
	            if ( aX < 1 ) {
	                count += 1;
	            } else {
	            	 continuedFraction.add(0);
	                 continuedFraction.set(i, count);

	                 i += 1;
	                 count = 1;
	                 current = 1;
	                 aX -= 1;
	            }
	        } else {
	        	 if ( aX > 1 ) {
	                 count += 1;
	                 aX -= 1;
	        	 } else {
	                 continuedFraction.add(0);
	                 continuedFraction.set(i, count);

	                 i += 1;
	                 count = 1;
	                 current = 0;
	        	 }
	        }

	        if ( aX == Math.floor(aX) ) {
	             continuedFraction.set(i, count);
	             break;
	        }

	        if ( i == MAX_ITERATIONS ) {
	             break;
	        }
	    }

	    double reciprocal = 1.0 / continuedFraction.get(i);
	    for ( int j = i - 1; j >= 0; j-- ) {
	         reciprocal = continuedFraction.get(j) + 1.0 / reciprocal;
	    }

	    return 1.0 / reciprocal;	           
	}
	
	private static final int MAX_ITERATIONS = 150;

}

Output:

  1.6666666666696983  1.6666666666666667
 -0.5657414540893351 -0.5657414540893352
  0.7182818280000092  0.7182818280000000
  0.1213141516271825  0.1213141516171819

Julia

function questionmark(x)
    y, p = fldmod(x, 1)
    q, d = 1 - p, .5
    while y + d > y
        p < q ? (q -= p) : (p -= q; y += d)
        d /= 2
    end
    y
end

function questionmark_inv(x)
    y, bits = fldmod(x, 1)
    lo, hi = [0, 1], [1, 1]
    while (y + /(lo...)) < (y + /(hi...))
        bit, bits = fldmod(2bits, 1) 
        bit > 0 ? (lo .+= hi) : (hi .+= lo)
    end
    y + /(lo...)
end

x, y = 0.7182818281828, 0.1213141516171819
for (a, b) ∈ [
        (5/3, questionmark((1 + √5)/2)),
        ((√13-7)/6, questionmark_inv(-5/9)),
        (x, questionmark_inv(questionmark(x))),
        (y, questionmark(questionmark_inv(y)))]
    println(a, a ≈ b ? " ≈ " : " != ", b)
end

Output:

1.6666666666666667 ≈ 1.666666666667894
-0.5657414540893352 ≈ -0.5657414540893351
0.7182818281828 ≈ 0.7182818281828183
0.1213141516171819 ≈ 0.12131415161718095

Mathematica / Wolfram Language

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

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

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

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

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

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; $j ≥ 0; $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

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

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