Minkowski question-mark function: Difference between revisions
Line 266: | Line 266: | ||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
{{trans| |
{{trans|FreeBASIC}} |
||
<lang julia>function minkowski(x) |
<lang julia>function minkowski(x) |
||
p = Int(floor(x)) |
p = Int(floor(x)) |
Revision as of 21:42, 14 November 2020
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
- 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
<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
<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
<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
<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
<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