Minkowski question-mark function

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

Factor

<lang factor>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</lang>

Output:

1.666666666668335        1.666666666666667        
-0.5657414540893351      -0.5657414540893352      
0.718281828000002        0.1213141516171819

FreeBASIC

<lang 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)) </lang>

Output:

 1.666666666669698           1.666666666666667
-0.5657414540893351         -0.5657414540893352
 0.7182818280000092          0.1213141516171819

Go

Translation of: FreeBASIC

<lang go>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)))

}</lang>

Output:

 1.6666666666696983  1.6666666666666667
-0.5657414540893351 -0.5657414540893352
 0.7182818280000092  0.1213141516171819

Julia

Translation of: FreeBASIC

<lang julia>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 true
       x *= 2
       if curr == 0
           if x < 1
               coun += 1
           else
               i += 1
               append!(contfrac, fill(0, i - length(contfrac) + 1))
               contfrac[i] = coun
               coun = 1
               curr = 1
               x -= 1
           end
       else
           if x > 1
               coun += 1
               x -= 1
           else
               i += 1
               append!(contfrac, fill(0, i - length(contfrac) + 1))
               contfrac[i] = coun
               coun = 1
               curr = 0
           end
       end
       if x == Int(floor(x))
           contfrac[i + 1] = coun
           break
       end
       i == maxiter && break
   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)))

</lang>

Output:

 1.6666666666696983   1.6666666666666667
-0.5657414540893351   -0.5657414540893352
 0.7182818280000092   0.12131415161718191

Phix

Translation of: FreeBASIC

<lang Euphoria>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))})</lang>

Output:

  1.6666666666696983   1.6666666666666668
 -0.5657414540893351  -0.5657414540893352
  0.7182818280000092   0.1213141516171819

REXX

Translation of: FreeBASIC

Translation of: Phix

<lang rexx>/*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</lang>

output when using the internal default inputs:

     1.666666666666666963      1.666666666666666667
    -0.565741454089335118     -0.565741454089335118
     0.718281828000000011      0.121314151617181900

Wren

Translation of: FreeBASIC

Library: Wren-fmt

<lang ecmascript>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)))</lang>

Output:

 1.66666666666970  1.66666666666667
-0.56574145408934 -0.56574145408934
 0.71828182800001  0.12131415161718