gcd with doubles

[Moritz Maxeiner via Digitalmars-d-learn]

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