Re: gcd with doubles
On Friday, 1 September 2017 at 09:33:08 UTC, Alex wrote: On Sunday, 27 August 2017 at 23:13:24 UTC, Moritz Maxeiner wrote: On Sunday, 27 August 2017 at 19:47:59 UTC, Alex wrote: [...] 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): [...] Hey, cool! Thanks for the efforts :) No problem, two corrections to myself, though: 1) It's a lower bound, not an upper bound (you need at least that many digits in order to not lose precision) 2) The code is missing `_ > ulong.min` checks along the existing `_ < ulong.max` checks
Re: gcd with doubles
On Sunday, 27 August 2017 at 23:13:24 UTC, Moritz Maxeiner wrote: On Sunday, 27 August 2017 at 19:47:59 UTC, Alex wrote: [...] 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): [...] Hey, cool! Thanks for the efforts :)
Re: gcd with doubles
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.001);
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;
Re: gcd with doubles
On Sunday, 27 August 2017 at 19:47:59 UTC, Alex wrote:
Hi, all.
Can anybody explain to me why
void main()
{
import std.numeric;
assert(gcd(0.5,32) == 0.5);
assert(gcd(0.2,32) == 0.2);
}
fails on the second assert?
I'm aware, that calculating gcd on doubles is not so obvios, as
on integers. But if the library accepts doubles, and basically
the return is correct occasionally, why it is not always the
case?
If the type isn't a builtin integral and can't be bit shifted,
the gcd algorithm falls back to using the Euclidean algorithm in
order to support custom number types, so the above gdc in the
above reduces to:
---
double gcd(double a, double b)
{
while (b != 0)
{
auto t = b;
b = a % b;
a = t;
}
return a;
}
---
The issue boils down to the fact that `32 % 0.2` yield `0.2`
instead of `0.0`, so the best answer I can give is "because
floating points calculations are approximations". I'm actually
not sure if this is a bug in fmod or expected behaviour, but I'd
tend to the latter.
Is there a workaround, maybe?
If you know how many digits of precision after the decimal dot
you can multiply beforehand, gcd in integer realm, and div
afterwards (be warned, the below is only an example
implementation for readability, it does not do the required
overflow checks for the double -> ulong conversion!):
---
import std.traits : isFloatingPoint;
T gcd(ubyte precision, T)(T a, T b) if (isFloatingPoint!T)
{
import std.numeric : _gcd = gcd;
immutable T coeff = 10 * precision;
return (cast(T) _gcd(cast(ulong) (a * coeff),
cast(ulong) (b * coeff))) / coeff;
}
---
Re: gcd with doubles
ok... googled a little bit. Seems to be the problem of the kind floating poing <--> decimal... Will try to get by with some division and lrint logic, as there is a lack of definition, how to define a gcd in general case. For example with 30 and 0.16. Never mind...
gcd with doubles
Hi, all.
Can anybody explain to me why
void main()
{
import std.numeric;
assert(gcd(0.5,32) == 0.5);
assert(gcd(0.2,32) == 0.2);
}
fails on the second assert?
I'm aware, that calculating gcd on doubles is not so obvios, as
on integers. But if the library accepts doubles, and basically
the return is correct occasionally, why it is not always the case?
Is there a workaround, maybe?
