digitalmars.D.learn - gcd with doubles

• Alex (13/13) Aug 27 2017 Hi, all.
• Alex (6/6) Aug 27 2017 ok... googled a little bit.
• Moritz Maxeiner (37/51) Aug 27 2017 If the type isn't a builtin integral and can't be bit shifted,
• Moritz Maxeiner (178/180) Aug 27 2017 To expand on the earlier workaround: You can also adapt a
• Alex (3/12) Sep 01 2017 Hey, cool!
• Moritz Maxeiner (6/21) Sep 01 2017 No problem, two corrections to myself, though:
```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?
```
Aug 27 2017
```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...
```
Aug 27 2017
```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;
}
---
```
Aug 27 2017
```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 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;
}
}
}
---

 https://github.com/marcandrysco/Errol
```
Aug 27 2017
```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 for that purpose (MIT license as
that is what the original code is under):

[...]

Hey, cool!
Thanks for the efforts :)
```
Sep 01 2017
```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 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
```
Sep 01 2017