## 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[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

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:Hey, cool! Thanks for the efforts :)[...]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): [...]

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: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` checksOn Sunday, 27 August 2017 at 19:47:59 UTC, Alex wrote:Hey, cool! Thanks for the efforts :)[...]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): [...]

Sep 01 2017