Thread overview
Exact arithmetic with quadratic irrationals
Apr 19, 2017
H. S. Teoh
Apr 19, 2017
Stanislav Blinov
Apr 19, 2017
Timon Gehr
April 19, 2017
I alluded to this in D.learn some time ago, and finally decided to take
the dip and actually write the code. So here it is: exact arithmetic
with numbers of the form (a+b√r)/c where a, b, c are integers, c!=0, and
r is a (fixed) square-free integer.

Code:	https://github.com/quickfur/qrat

I wrote this code in just a little over a day, a testament to D's productivity-increasing features, among which include:

1) Built-in unittests: I couldn't have had the confidence that my code was correct if I hadn't been able to verify it using unittests, that also ensured that there would be no regressions, *and* also provided nice ddoc'd examples for the user. This is a killer combo, IMO.

2) Sane(r) operator overloading: even though there's still a certain amount of boilerplate, operator overloading in D is far saner than in C++, thanks to:

	a) Templated opUnary / opBinary / opOpAssign with the operator
	as a string template argument that can be used in mixins. This
	saved me quite a bit of copy-pasta that would have otherwise
	been necessary, e.g., in a C++ implementation.

	b) Unification of <, <=, >, >= under opCmp, and !=, == under
	opEquals. Again, tons of copy-pasta were avoided where they
	would have been necessary in C++.

	c) opBinaryRight, an underrated genius design decision, that
	allowed easy support of expressions of the form `1 + q` without
	C++ hacks like friend functions and what-not.

3) Sane template syntax, that, combined with opBinaryRight and the other
overloading features, made expressions like this possible, and
readable(!):

	// So clear, so readable!
	auto goldenRatio = (1 + surd!5)/2;

	// This would have been a mess in C++ syntax due to template<x>
	// clashing visually with operator <.
	assert((10 + surd!5)/20 < (surd!5 - 1)/2);

4) Code coverage built into the compiler with -cov, that caught at least one bug in a piece of code that my unittests missed. Now with 100% code coverage, I feel far more confident that there are no nasty bugs left!

5) Pay-as-you-go template methods: I deliberately turned all QRat
methods into template functions, because of (a) template attribute
inference (see below), and (b) reducing template bloat from
instantiating methods that are never used in user code.

6) Template attribute inference: this allowed me, *after the fact*, to slap 'pure nothrow @safe @nogc' onto my module ddoc'd unittest, and instantly get compiler feedback on where exactly I'd inadvertently broke purity / nothrowness / safety / etc., where I didn't intend.  And it was very gratifying, once I'd isolated those bits of code, to have the confidence that the compiler has verified that the core operations (unary / binary operators, comparisons, etc.) were all pure nothrow @safe @nogc, and thus widely usable.  If I'd had to wrangle with explicitly writing attributes, I would've been greatly hampered in productivity (not to mention frustrated and distracted from focusing on the actual algorithms rather than the fussy nitty-gritties of the language).  Attribute inference is definitely the way of the future, and I support Walter in pushing it as far as it can possibly go. And statically-verified guarantees too.  That's another win for D.

Having said that, though, there were a few roadblocks:

1) std.numeric.gcd doesn't support BigInt. This in itself wouldn't have been overly horrible, if it weren't for the fact that:

	a) The algorithm in std.numeric.gcd actually *does* work for
	BigInt, but BigInt support wasn't implemented because there are
	ostensibly better-performing BigInt GCD algorithms out there --
	but they are too complex to implement so nobody has done it yet.
	An unfortunate example of letting the perfect being the enemy of
	the good, that's been plaguing D for a while now.

		https://issues.dlang.org/show_bug.cgi?id=7102

	b) There are no sig constraints in std.numeric.gcd even though
	it doesn't support BigInts (or, for that matter, custom
	numerical types). This means that once you import std.numeric,
	you can't even provide your own overloads of gcd to handle
	BigInt or whatever other types you wish to handle, without
	running into overload conflicts. Not without unnecessarily
	convoluted schemes of static imports, symbol aliasing, and all
	the rest of that churn.

	c) The only thing that's really standing in the way of BigInt in
	std.numeric.gcd is an ill-advised (IMO) way to test if a numeric
	type has sign, by assuming that that's equivalent to having a
	.min property.  That really only works for built-in types, which
	again screams "missing sig constraints!", and basically fails
	for everything else.

2) std.numeric.gcd isn't variadic, so I had to roll my own. Not a biggie, but it was still annoying and wasted my time writing code that calls gcd(a,b) multiple times when I first started coding, only to later realize I'd be doing this a *lot* and deciding to write variadic gcd myself.

Haha, it seems that the only roadblocks were related to the implementation quality of std.numeric.gcd... nothing that a few relatively-simple PRs couldn't fix.  So overall, D is still awesome.


T

-- 
Those who don't understand Unix are condemned to reinvent it, poorly.
April 19, 2017
Awesome! Congrats and thanks for sharing.

On Wednesday, 19 April 2017 at 19:32:14 UTC, H. S. Teoh wrote:

> Haha, it seems that the only roadblocks were related to the implementation quality of std.numeric.gcd... nothing that a few relatively-simple PRs couldn't fix.  So overall, D is still awesome.

There's another one, which is more about dmd: dmd does not inline gcd, which, when arguments are const, turns gcd into a double function call :D

April 19, 2017
On 19.04.2017 21:32, H. S. Teoh via Digitalmars-d wrote:
> I alluded to this in D.learn some time ago, and finally decided to take
> the dip and actually write the code. So here it is: exact arithmetic
> with numbers of the form (a+b√r)/c where a, b, c are integers, c!=0, and
> r is a (fixed) square-free integer.
>
> Code:	https://github.com/quickfur/qrat
>
> ...

Nice. :)

Some suggestions:

- You might want to support ^^ (it is useful for examples like the one below).

- constructor parameter _b should have a default value of 0.

- formatting should special case b==-1 like it special cases b==1.
  (also: as you are using Unicode anyway, you could also use · instead of *. Another cute thing to do is to add a vinculum.)


Example application: Computing large Fibonacci numbers efficiently:

import qrat;
import std.bigint;
alias ℕ=BigInt;
enum φ=(1+surd!(5,ℕ))/2,ψ=(1-surd!(5,ℕ))/2;

auto pow(T,S)(T a,S n){
    T r=T(ℕ(1),ℕ(0));
    for(auto x=a;n;n>>=1,a*=a)
        if(n&1) r*=a;
    return r;
}

auto fib(long n){
    return (pow(φ,n)-pow(ψ,n))/surd!(5,ℕ);
}
void main(){
    import std.stdio;
    foreach(i;0..40) writeln(fib(i));
    writeln(fib(100000));
}