www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.learn - gcd with doubles

reply Alex <sascha.orlov gmail.com> writes:
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
next sibling parent Alex <sascha.orlov gmail.com> writes:
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
prev sibling next sibling parent Moritz Maxeiner <moritz ucworks.org> writes:
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
prev sibling parent reply Moritz Maxeiner <moritz ucworks.org> writes:
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
parent reply Alex <sascha.orlov gmail.com> writes:
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 :)
Sep 01 2017
parent Moritz Maxeiner <moritz ucworks.org> writes:
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
Sep 01 2017