Minkowski question-mark function: Difference between revisions
Add C# implementation |
|||
(11 intermediate revisions by 4 users not shown) | |||
Line 128: | Line 128: | ||
</pre> |
</pre> |
||
=={{header|C#}}== |
|||
{{trans|Go}} |
|||
<syntaxhighlight lang="C#"> |
|||
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))); |
|||
} |
|||
} |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
1.6666666666697000 1.6666666666666700 |
|||
-0.5657414540893350 -0.5657414540893350 |
|||
0.7182818280000090 0.1213141516171820 |
|||
</pre> |
|||
=={{header|C++}}== |
|||
<syntaxhighlight lang="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; |
|||
} |
|||
</syntaxhighlight> |
|||
{{ out }} |
|||
<pre> |
|||
1.6666666666696983 1.6666666666666667 |
|||
-0.5657414540893351 -0.5657414540893352 |
|||
0.7182818281828269 0.7182818281828180 |
|||
0.1213141516171819 0.1213141516171819 |
|||
</pre> |
|||
=={{header|F_Sharp|F#}}== |
|||
<syntaxhighlight lang="fsharp"> |
|||
// 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}" |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
?(φ) = 5/3 is true |
|||
?⁻¹(-5/9) = (√13-7)/6 is true |
|||
?⁻¹(?(n)) = n is true |
|||
?(?⁻¹(n)) = n is true |
|||
</pre> |
|||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
{{works with|Factor|0.99 2020-08-14}} |
{{works with|Factor|0.99 2020-08-14}} |
||
Line 516: | Line 805: | ||
Using an exact fraction for 4.04145188432738056 and bumping the iteration count from 52 up to 200 changes that difference to 1.43622e_12. |
Using an exact fraction for 4.04145188432738056 and bumping the iteration count from 52 up to 200 changes that difference to 1.43622e_12. |
||
=={{header|Java}}== |
|||
<syntaxhighlight lang="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; |
|||
} |
|||
</syntaxhighlight> |
|||
{{ out }} |
|||
<pre> |
|||
1.6666666666696983 1.6666666666666667 |
|||
-0.5657414540893351 -0.5657414540893352 |
|||
0.7182818280000092 0.7182818280000000 |
|||
0.1213141516271825 0.1213141516171819 |
|||
</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
Line 532: | Line 949: | ||
y, bits = fldmod(x, 1) |
y, bits = fldmod(x, 1) |
||
lo, hi = [0, 1], [1, 1] |
lo, hi = [0, 1], [1, 1] |
||
while (y + /(lo...)) |
while (y + /(lo...)) < (y + /(hi...)) |
||
bit, bits = fldmod( |
bit, bits = fldmod(2bits, 1) |
||
bit > 0 ? (lo .+= hi) : (hi .+= lo) |
bit > 0 ? (lo .+= hi) : (hi .+= lo) |
||
end |
end |
||
y + / |
y + /(lo...) |
||
end |
end |
||
x, y = 0. |
x, y = 0.7182818281828, 0.1213141516171819 |
||
for (a, b) ∈ [ |
for (a, b) ∈ [ |
||
(5/3, questionmark((1 + √5)/2)), |
(5/3, questionmark((1 + √5)/2)), |
||
Line 549: | Line 966: | ||
</syntaxhighlight>{{out}} |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
1.6666666666666667 ≈ 1.666666666667894 |
|||
-0.5657414540893352 ≈ -0.5657414540893351 |
-0.5657414540893352 ≈ -0.5657414540893351 |
||
0. |
0.7182818281828 ≈ 0.7182818281828183 |
||
0.1213141516171819 ≈ 0.12131415161718095 |
0.1213141516171819 ≈ 0.12131415161718095 |
||
</pre> |
</pre> |
||
Line 1,069: | Line 1,486: | ||
{{trans|FreeBASIC}} |
{{trans|FreeBASIC}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
<syntaxhighlight lang=" |
<syntaxhighlight lang="wren">import "./fmt" for Fmt |
||
var MAXITER = 151 |
var MAXITER = 151 |
Latest revision as of 11:44, 12 February 2024
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
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
- Wikipedia entry: Minkowski's question-mark function
11l
-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#
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
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
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
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
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
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
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
# 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
/*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
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