On Sunday, 27 August 2017 at 19:47:59 UTC, Alex wrote:
[..]
Is there a workaround, maybe?
To expand on the earlier workaround: You can also adapt a
floating point to string algorithm in order to dynamically
determine an upper bound on the number of after decimal point
digits required. Below is an untested adaption of the reference C
implementation of errol0[1] for that purpose (MIT license as that
is what the original code is under):
---
void main()
{
assert(gcd(0.5, 32) == 0.5);
assert(gcd(0.2, 32) == 0.2);
assert(gcd(1.3e2, 3e-5) == 1e-5);
}
template gcd(T)
{
import std.traits : isFloatingPoint;
T gcd(T a, T b)
{
static if (isFloatingPoint!T)
{
return fgcd(a, b);
}
else
{
import std.numeric : igcd = gcd;
return igcd(a, b);
}
}
static if (isFloatingPoint!T)
{
import std.math : nextUp, nextDown, pow, abs, isFinite;
import std.algorithm : max;
T fgcd(T a, T b)
in
{
assert (a.isFinite);
assert (b.isFinite);
assert (a < ulong.max);
assert (b < ulong.max);
}
body
{
short a_exponent;
int a_digitCount = errol0CountOnly(abs(a),
a_exponent);
short b_exponent;
int b_digitCount = errol0CountOnly(abs(b),
b_exponent);
a_digitCount -= a_exponent;
if (a_digitCount < 0)
{
a_digitCount = 0;
}
b_digitCount -= b_exponent;
if (b_digitCount < 0)
{
b_digitCount = 0;
}
auto coeff = pow(10, max(a_digitCount, b_digitCount));
assert (a * coeff < ulong.max);
assert (b * coeff < ulong.max);
return (cast(T) euclid(cast(ulong) (a * coeff),
cast(ulong) (b * coeff))) /
coeff;
}
ulong euclid(ulong a, ulong b)
{
while (b != 0)
{
auto t = b;
b = a % b;
a = t;
}
return a;
}
struct HighPrecisionFloatingPoint
{
T base, offset;
void normalize()
{
T base = this.base;
this.base += this.offset;
this.offset += base - this.base;
}
void mul10()
{
T base = this.base;
this.base *= T(10);
this.offset *= T(10);
T offset = this.base;
offset -= base * T(8);
offset -= base * T(2);
this.offset -= offset;
normalize();
}
void div10()
{
T base = this.base;
this.base /= T(10);
this.offset /= T(10);
base -= this.base * T(8);
base -= this.base * T(2);
this.offset += base / T(10);
normalize();
}
}
alias HP = HighPrecisionFloatingPoint;
enum epsilon = T(0.0000001);
ushort errol0CountOnly(T f, out short exponent)
{
ushort digitCount;
T ten = T(1);
exponent = 1;
auto mid = HP(f, T(0));
while (((mid.base > T(10)) || ((mid.base == T(10)) &&
(mid.offset >= T(0)))) && (exponent < 308))
{
exponent += 1;
mid.div10();
ten /= T(10);
}
while (((mid.base < T(1)) || ((mid.base == T(1)) &&
(mid.offset < T(0)))) && (exponent > -307))
{
exponent -= 1;
mid.mul10();
ten *= T(10);
}
auto inhi = HP(mid.base, mid.offset + (nextUp(f) - f)
* ten / (T(2) + epsilon));
auto inlo = HP(mid.base, mid.offset + (nextDown(f) -
f) * ten / (T(2) + epsilon));
inhi.normalize();
inlo.normalize();
while (inhi.base > T(10) || (inhi.base == T(10) &&
(inhi.offset >= T(0))))
{
exponent += 1;
inhi.div10();
inlo.div10();
}
while (inhi.base < T(1) || (inhi.base == T(1) &&
(inhi.offset < T(0))))
{
exponent -= 1;
inhi.mul10();
inlo.mul10();
}
while (inhi.base != T(0) || inhi.offset != T(0))
{
auto hdig = cast(ubyte) inhi.base;
if ((inhi.base == hdig) && (inhi.offset < T(0)))
{
hdig -= 1;
}
auto ldig = cast(ubyte) inlo.base;
if ((inlo.base == ldig) && (inlo.offset < 0))
{
ldig -= 1;
}
if (ldig != hdig)
{
break;
}
digitCount += 1;
inhi.base -= hdig;
inhi.mul10();
inlo.base -= ldig;
inlo.mul10();
}
digitCount += 1;
return digitCount;
}
}
}
---
[1] https://github.com/marcandrysco/Errol
