www.digitalmars.com         C & C++   DMDScript  

digitalmars.D - Exact arithmetic with quadratic irrationals

reply "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
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.
Apr 19 2017
next sibling parent reply Stanislav Blinov <stanislav.blinov gmail.com> writes:
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
Apr 19 2017
parent "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
On Wed, Apr 19, 2017 at 07:54:02PM +0000, Stanislav Blinov via Digitalmars-d
wrote:

 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


If I weren't such a sucker for the bleeDing edge with dmd (I actually
compile even my serious projects with git HEAD dmd, except when
performance matters), I'd be standardizing on ldc or gdc, which have far
superior optimizers.

I consistently find that my CPU-intensive projects perform at least
20-30% worse on dmd than gdc (and I presume ldc), sometimes even as bad
as 40-50%, due to dmd's inliner giving up far too easily. I don't know
if this has been fixed yet, but the last time I checked, if you wrote
this:

	int func(int x, int y) {
		if (x<0)
			return y;
		return x;
	}

instead of:

	int func(int x, int y) {
		if (x<0)
			return y;
		else
			return x;
	}

then the dmd inliner would not inline the function.

Because of sensitivities like this, the inliner gives up far too early
in the process, thus the optimizer misses out on further optimization
opportunities that would have opened up, had the function been inlined.

The last time I checked, I also found that dmd was rather weak at loop
optimizations (and loops are very important in performance as we all
know) compared to gdc. Again, the same domino effect (or rather, the
missing thereof) applies: by failing to, for example, hoist a constant
expression out of the loop, further optimization opportunities are
missed, whereas gdc, after hoisting the expression out, would discover
that the loop can be reduced further, perhaps via a strength reduction,
and then unrolled, and then inlined inside the caller, then vectorized,
etc.. This chain of optimizations were missed because of one missed
optimization early in the process. Hence the suboptimal resulting code.


T

-- 
If blunt statements had a point, they wouldn't be blunt...
Apr 19 2017
prev sibling next sibling parent reply Timon Gehr <timon.gehr gmx.ch> writes:
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));
}
Apr 19 2017
parent reply "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
On Wed, Apr 19, 2017 at 10:47:04PM +0200, Timon Gehr via Digitalmars-d wrote:

 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).


I would, except that I doubt it would perform any better than an actual
recursive or iterative algorithm for computing Fibonacci sequences,
because I don't know of any simple way to exponentiate a quadratic
rational using only integer arithmetic other than repeated
multiplication.  (For all I know, it may perform even worse, because
multiplying n quadratic rationals involves a lot more than just summing
n+1 integers as in an iterative implementation of Fibonacci.)

Hmm, come to think of it, I *could* expand the numerator using the
binomial theorem, treating (a+b√r) as a binomial (a+bx) where x=√r, and
folding even powers into the integral part (since x^2 = r, so x^(2k) =
r^k). The denominator could be exponentiated using plain ole integer
exponentiation.  Then it's just a matter of summing coefficients.

But it still seems to be about the same amount of work as (or more than)
summing n+1 integers in an iterative implementation of Fibonacci.  Am I
missing something?



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


Good point, done.



 - formatting should special case b==-1 like it special cases b==1.


Done, good catch!



   (also: as you are using Unicode anyway, you could also use · instead
   of *.  Another cute thing to do is to add a vinculum.)


Well, I would, but it gets a bit too fancy for my tastes and may not
render well on all displays. But I'll put it on my list of things to
consider.

Another module I'm thinking about is an extension of QRat that allows
you to mix different radicals in the same expression. For example,
(√3+√5)/√7 and so forth.  I have discovered algorithms that, given n
distinct radicals, allow a closed-form expression with 2^n coefficients
(+1 denominator), closed under field operations.  The only trouble in
this case is that reciprocating such things will be very slow, as will
comparisons, and both have a high chance of overflow (so BigInt is
probably a necessity).  And 2^n+1 coefficients for n radicals quickly
gets expensive space-wise as n increases.

Yesterday I also discovered an algorithm for expressing the reciprocal
of numbers of the form:

	(a + b∛r + c∛r^2)/d

in the same form. I.e., for rewriting:

	d/(a + b∛r + c∛r^2)

in the first form.  Which means it's possible to implement a QRat-like
representation for cubic rationals.  (The actual computation is rather
expensive, as it involves quite a lot of multiplications, squaring, and
cubing. But it's *possible*.)  I'm still trying to verify the
correctness of the formula I obtained, since while checking a concrete
example last night I discovered a possible error.

If this works out, I might consider 4th roots as well -- though I'm
expecting that might be near the limit of the usefulness of these
representations, since the complexity becomes so great that symbolic
manipulation like in Mathematica may turn out to be more feasible after
all.  But it may be of some theoretical interest whether such
representations are possible, even if they are ultimately impractical. A
particularly interesting question is whether such representations exist
for *all* algebraic numbers (of bounded degree).

Currently I have a conjecture that given a rational extension of n
radicals of degree k, field closure can be achieved with a
representation of k^n+1 coefficients. But it's still too early to say
whether algorithms exist for inverting radicals of degree k for large k
-- I have a creeping suspicion that perhaps somewhere around k=5 or k=6
the unsolvability of the general quintic may raise its ugly head and
prevent further progress.


T

-- 
INTEL = Only half of "intelligence".
Apr 19 2017
parent reply Timon Gehr <timon.gehr gmx.ch> writes:
On 19.04.2017 23:39, H. S. Teoh via Digitalmars-d wrote:

 On Wed, Apr 19, 2017 at 10:47:04PM +0200, Timon Gehr via Digitalmars-d wrote:

 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).


 I would, except that I doubt it would perform any better than an actual
 recursive or iterative algorithm for computing Fibonacci sequences,
 because I don't know of any simple way to exponentiate a quadratic
 rational using only integer arithmetic other than repeated
 multiplication.  (For all I know, it may perform even worse, because
 multiplying n quadratic rationals involves a lot more than just summing
 n+1 integers as in an iterative implementation of Fibonacci.)

 Hmm, come to think of it, I *could* expand the numerator using the
 binomial theorem, treating (a+b√r) as a binomial (a+bx) where x=√r, and
 folding even powers into the integral part (since x^2 = r, so x^(2k) =
 r^k). The denominator could be exponentiated using plain ole integer
 exponentiation.  Then it's just a matter of summing coefficients.

 But it still seems to be about the same amount of work as (or more than)
 summing n+1 integers in an iterative implementation of Fibonacci.  Am I
 missing something?
 ...


Yes, there is in fact a beautifully simple way to do better. :)

Assume we want to compute some power of x. With a single multiplication, 
we obtain x². Multiplying x² by itself, we obtain x⁴. Repeating this a 
few times, we get:

x, x², x⁴, x⁸, x¹⁶, x³², etc.

In general, we only need n operations to compute x^(2ⁿ).

To compute xʸ, it basically suffices to express y as a sum of powers of 
two (i.e., we write it in binary).

For example, 22 = 16 + 4 + 2, and x²² = x¹⁶·x⁴·x².

My last post includes an implementation of this algorithm. ;)

In particular, this leads to multiple ways to compute the n-th Fibonacci 
number using O(log n) basic operations. (One way is to use your QRat 
type, but we can also use the matrix (1 1; 1 0).)


 ...

 Another module I'm thinking about is an extension of QRat that allows
 you to mix different radicals in the same expression. For example,
 (√3+√5)/√7 and so forth. ...


That would certainly be nice. Note that QRat is basically already there 
for different quadratic radicals, the only reason it does not work 
already is that we cannot use a QRat as the base field instead of ℚ 
(because ℚ is hardcoded).


This is the relevant concept from algebra:
https://en.wikipedia.org/wiki/Splitting_field


All your conjectures are true, except the last one. (Galois theory is 
not an obstacle, because here, we only need to consider splitting fields 
of particularly simple polynomials that are solvable in radicals.) You 
can even mix radicals of different degrees.


To get the formula for multiplicative inverses, one possible algorithm is:
https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Polynomial_extended_Euclidean_algorithm
Apr 19 2017
next sibling parent Timon Gehr <timon.gehr gmx.ch> writes:
On 20.04.2017 02:01, Timon Gehr wrote:

 To get the formula for multiplicative inverses, one possible algorithm is:
 https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Polynomial_extended_Euclidean_algorithm


Better reference: 
https://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor#Arithmetic_of_algebraic_extensions
Apr 19 2017
prev sibling next sibling parent Timon Gehr <timon.gehr gmx.ch> writes:
On 20.04.2017 02:01, Timon Gehr wrote:

 My last post includes an implementation of this algorithm. ;)


But in that implementation I used the parameter 'a' instead of the 
variable 'x' as a result of being tired, which makes it slightly more 
confusing than necessary even though it is correct. More readable version:

auto pow(T,S)(T a,S n){
     T r=T(ℕ(1),ℕ(0));
     for(auto x=a;n;n>>=1,x*=x)
         if(n&1) r*=x;
     return r;
}
Apr 19 2017
prev sibling parent reply "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
On Thu, Apr 20, 2017 at 02:01:20AM +0200, Timon Gehr via Digitalmars-d wrote:
[...]

 Yes, there is in fact a beautifully simple way to do better. :)
 
 Assume we want to compute some power of x. With a single
 multiplication, we obtain x². Multiplying x² by itself, we obtain x⁴.
 Repeating this a few times, we get:
 
 x, x², x⁴, x⁸, x¹⁶, x³², etc.
 
 In general, we only need n operations to compute x^(2ⁿ).
 
 To compute xʸ, it basically suffices to express y as a sum of powers
 of two (i.e., we write it in binary).
 
 For example, 22 = 16 + 4 + 2, and x²² = x¹⁶·x⁴·x².
 
 My last post includes an implementation of this algorithm. ;)


Ahh, so *that's* what it's all about. I figured that's what I was
missing. :-D  Thanks, I'll include this in QRat soon.



 In particular, this leads to multiple ways to compute the n-th
 Fibonacci number using O(log n) basic operations. (One way is to use
 your QRat type, but we can also use the matrix (1 1; 1 0).)


True, though I'm still jealous that with transcendental functions like
with floating-point, one could ostensibly compute that in O(1).



 Another module I'm thinking about is an extension of QRat that
 allows you to mix different radicals in the same expression. For
 example, (√3+√5)/√7 and so forth. ...
 


 
 That would certainly be nice. Note that QRat is basically already
 there for different quadratic radicals, the only reason it does not
 work already is that we cannot use a QRat as the base field instead of
 ℚ (because ℚ is hardcoded).


Oh?  I didn't try it myself, but if QRat itself passes isArithmeticType,
I'd venture to say QRat!(n, QRat!m) ought to work... There are some
hidden assumptions about properties of the rationals, though, but I
surmise none that couldn't be replaced by prerequisites about the
relative linear dependence of the mixed radicals over Q.

What I had in mind, though, was a more direct approach that perhaps may
reduce the total number of operations, since if the code is aware that
multiple radicals are involved, it could potentially factor out some
commonalities to minimize recomputations.

Also, the current implementation of QRat fixes the radical at
compile-time; I wanted to see if I could dynamically handle arbitrary
radicals at runtime. It would have to be restricted by only allowing
operations between two QRats of the same extension, of course, but if
the code could handle arbitrary extensions dynamically, then that
restriction could be lifted and we could (potentially, anyway) support
arbitrary combinations of expressions involving radicals. That would be
far more useful than QRat, for some of the things I'd like to use it
for.



 This is the relevant concept from algebra:
 https://en.wikipedia.org/wiki/Splitting_field
 
 
 All your conjectures are true, except the last one. (Galois theory is
 not an obstacle, because here, we only need to consider splitting
 fields of particularly simple polynomials that are solvable in
 radicals.)


That's nice to know. But I suppose Galois theory *would* become an
obstacle if I wanted to implement, for example, Q(x) for an arbitrary
algebraic x?



 You can even mix radicals of different degrees.


Yes, I've thought about that before. So it should be possible to
represent Q(r1,r2,...rn) using deg(r1)*deg(r2)*...*deg(rn)+1
coefficients?



 To get the formula for multiplicative inverses, one possible algorithm is:
 https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Polynomial_extended_Euclidean_algorithm


[...]

Thanks, will look into this at some point. :-)


T

-- 
Some ideas are so stupid that only intellectuals could believe them. -- George
Orwell
Apr 19 2017
parent reply Timon Gehr <timon.gehr gmx.ch> writes:
On 20.04.2017 03:00, H. S. Teoh via Digitalmars-d wrote:

 On Thu, Apr 20, 2017 at 02:01:20AM +0200, Timon Gehr via Digitalmars-d wrote:
 [...]

 Yes, there is in fact a beautifully simple way to do better. :)
 ...


 Ahh, so *that's* what it's all about. I figured that's what I was
 missing. :-D  Thanks, I'll include this in QRat soon.



 In particular, this leads to multiple ways to compute the n-th
 Fibonacci number using O(log n) basic operations. (One way is to use
 your QRat type, but we can also use the matrix (1 1; 1 0).)


 True, though I'm still jealous that with transcendental functions like
 with floating-point, one could ostensibly compute that in O(1).
 ...


BTW, you are right that with std.bigint, computation using a linear 
number of additions is actually faster for my example (100000th 
Fibonacci number). The asymptotic running time of the version with pow 
on QRats is lower though, so there ought to be a crossover point. (It is 
Θ(n^2) vs. O(n^log₂(3)·log(n)). std.bigint does not implement anything 
that is asymptotically faster than Karatsuba.)

For computations over field extensions of (small) finite fields, pow is 
a lot faster though.


 Another module I'm thinking about is an extension of QRat that
 allows you to mix different radicals in the same expression. For
 example, (√3+√5)/√7 and so forth. ...


 That would certainly be nice. Note that QRat is basically already
 there for different quadratic radicals, the only reason it does not
 work already is that we cannot use a QRat as the base field instead of
 ℚ (because ℚ is hardcoded).


 Oh?  I didn't try it myself, but if QRat itself passes isArithmeticType,
 I'd venture to say QRat!(n, QRat!m) ought to work... There are some
 hidden assumptions about properties of the rationals, though, but I
 surmise none that couldn't be replaced by prerequisites about the
 relative linear dependence of the mixed radicals over Q.
 ...


The issue is that gcd does not work on QRats. If QRat had two 
coefficients from an arbitrary (possibly ordered) field instead of 
encoding rationals explicitly, I think it would work.



 What I had in mind, though, was a more direct approach that perhaps may
 reduce the total number of operations, since if the code is aware that
 multiple radicals are involved, it could potentially factor out some
 commonalities to minimize recomputations.
 ...


This is probably the case.


 Also, the current implementation of QRat fixes the radical at
 compile-time; I wanted to see if I could dynamically handle arbitrary
 radicals at runtime. It would have to be restricted by only allowing
 operations between two QRats of the same extension, of course, but if
 the code could handle arbitrary extensions dynamically, then that
 restriction could be lifted and we could (potentially, anyway) support
 arbitrary combinations of expressions involving radicals. That would be
 far more useful than QRat, for some of the things I'd like to use it
 for.
 ...


What applications do you have in mind? Computational geometry?



 This is the relevant concept from algebra:
 https://en.wikipedia.org/wiki/Splitting_field


 All your conjectures are true, except the last one. (Galois theory is
 not an obstacle, because here, we only need to consider splitting
 fields of particularly simple polynomials that are solvable in
 radicals.)


 That's nice to know. But I suppose Galois theory *would* become an
 obstacle if I wanted to implement, for example, Q(x) for an arbitrary
 algebraic x?
 ...


All that the result about the quintic really says is that you will not, 
in general, be able to express x using field operations on radicals. It 
is still possible to compute the roots to arbitrary precision.
Computing the field operations in ℚ(x) will actually still be quite 
straightforward but you'd have to think about what to do with toString 
and opCmp. (Or more generally, you'd have to think about how to pick one 
of the roots of a given polynomial.)


 You can even mix radicals of different degrees.


 Yes, I've thought about that before. So it should be possible to
 represent Q(r1,r2,...rn) using deg(r1)*deg(r2)*...*deg(rn)+1
 coefficients?
 ...


Yes, at most, except you don't need "+1". (For each radical ri, you will 
at most need to pick a power between 0 to deg(ri)-1 to index into the 
coefficients.)
Apr 20 2017
parent reply "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
On Thu, Apr 20, 2017 at 02:51:12PM +0200, Timon Gehr via Digitalmars-d wrote:

 On 20.04.2017 03:00, H. S. Teoh via Digitalmars-d wrote:

 On Thu, Apr 20, 2017 at 02:01:20AM +0200, Timon Gehr via Digitalmars-d wrote:
 [...]

 Yes, there is in fact a beautifully simple way to do better. :)
 ...


 
 Ahh, so *that's* what it's all about. I figured that's what I was
 missing. :-D  Thanks, I'll include this in QRat soon.




Update: QRat now supports ^^. :-) Integral exponents only, of course.  I
also implemented negative exponents, because QRat supports division and
the same algorithm can be easily reused for that purpose.

Interestingly enough, std.math has an overload of pow() that pretty much
does exactly the same thing, except that its sig constraints require a
built-in floating-point type. I'm half-tempted to submit a Phobos PR to
relax the sig constraints so that we could actually use it for QRat
without essentially duplicating the code.


[...]

 BTW, you are right that with std.bigint, computation using a linear
 number of additions is actually faster for my example (100000th
 Fibonacci number).  The asymptotic running time of the version with
 pow on QRats is lower though, so there ought to be a crossover point.
 (It is Θ(n^2) vs.  O(n^log₂(3)·log(n)). std.bigint does not implement
 anything that is asymptotically faster than Karatsuba.)


Yeah, probably there is a crossover point. But it might be quite large.
I suppose one could make a graph of the running times for increasing n,
and either find the crossover point that way or extrapolate using the
known curve shapes.

Having said that, I haven't scrutinized the performance characteristics
of QRat too carefully just yet -- there is probably room for
optimization.


[...]

 That would certainly be nice. Note that QRat is basically already
 there for different quadratic radicals, the only reason it does
 not work already is that we cannot use a QRat as the base field
 instead of ℚ (because ℚ is hardcoded).


 
 Oh?  I didn't try it myself, but if QRat itself passes
 isArithmeticType, I'd venture to say QRat!(n, QRat!m) ought to
 work... There are some hidden assumptions about properties of the
 rationals, though, but I surmise none that couldn't be replaced by
 prerequisites about the relative linear dependence of the mixed
 radicals over Q.  ...


 
 The issue is that gcd does not work on QRats. If QRat had two
 coefficients from an arbitrary (possibly ordered) field instead of
 encoding rationals explicitly, I think it would work.


You're right, without gcd it won't work.  The current implementation is
a bit overzealous on using gcd (cf. the division algorithm), mainly
because I'm concerned with integer overflow on native types. Probably
some of these uses can be dispensed with, where BigInt is involved. But
normalize() still needs gcd, otherwise sgn() and opCmp() may produce
wrong results.

One thought is that if there is a QInt base type (i.e., implementing
numbers of the form a+b√r, without the denominator), then we could
implement a gcd algorithm for it, and we'd be able to instantiate
QRat!(r, QInt).



 What I had in mind, though, was a more direct approach that perhaps
 may reduce the total number of operations, since if the code is
 aware that multiple radicals are involved, it could potentially
 factor out some commonalities to minimize recomputations.  ...


 
 This is probably the case.


Upon reviewing the algorithms I've come up with in the past, it appears
that QRat!(r, QInt) may in fact produce essentially the same code.



 Also, the current implementation of QRat fixes the radical at
 compile-time; I wanted to see if I could dynamically handle
 arbitrary radicals at runtime. It would have to be restricted by
 only allowing operations between two QRats of the same extension, of
 course, but if the code could handle arbitrary extensions
 dynamically, then that restriction could be lifted and we could
 (potentially, anyway) support arbitrary combinations of expressions
 involving radicals. That would be far more useful than QRat, for
 some of the things I'd like to use it for.  ...


 
 What applications do you have in mind? Computational geometry?


Yes. In particular, manipulating the coordinates of certain kinds of
polytopes. I currently do have code that can do this with
floating-point, but I'd like to be able to deal with exact coordinates
rather than floating-point approximations.


[...]

 All your conjectures are true, except the last one. (Galois theory
 is not an obstacle, because here, we only need to consider
 splitting fields of particularly simple polynomials that are
 solvable in radicals.)


 
 That's nice to know. But I suppose Galois theory *would* become an
 obstacle if I wanted to implement, for example, Q(x) for an
 arbitrary algebraic x?  ...


 
 All that the result about the quintic really says is that you will
 not, in general, be able to express x using field operations on
 radicals. It is still possible to compute the roots to arbitrary
 precision.


Oh I know that; I'm not really concerned with computing roots to
arbitrary precision here though, but more with implementing precise
arithmetic on expressions involving said roots.



 Computing the field operations in ℚ(x) will actually still be quite
 straightforward but you'd have to think about what to do with toString
 and opCmp. (Or more generally, you'd have to think about how to pick
 one of the roots of a given polynomial.)


Hmm, good point.  I suppose I haven't really thought through the
consequences of what might happen if I implemented a reciprocation
algorithm for an algebraic number k where the defining polynomial for k
may have multiple roots. At some point, assumptions about which root is
being used need to come into play, I suppose.  How to encode this in the
API and internally in the code is an interesting question.



 You can even mix radicals of different degrees.


 
 Yes, I've thought about that before. So it should be possible to
 represent Q(r1,r2,...rn) using deg(r1)*deg(r2)*...*deg(rn)+1
 coefficients?
 ...


 
 Yes, at most, except you don't need "+1". (For each radical ri, you
 will at most need to pick a power between 0 to deg(ri)-1 to index into
 the coefficients.)


[...]

The +1 is for the denominator, assuming integer coefficients.  Since
having 2^n rational coefficients is equivalent to having 2^n integer
coefficients (which are half the size of rational coefficients,
computer-representation-wise) + 1 denominator.

Though it's arguable whether this really saves that much once you get
out of the realm of native machine integer types into BigInt. It's
debatable whether that much is really saved in terms of space and CPU
time if you have two or more rational coefficients with denominators
that are relatively prime to each other, so that merging them all into a
single denominator may just cause all coefficients to explode in size
whereas separate rational coefficients could in fact be more compact.

But, at least in theory, if you're dealing with relatively small members
in the field, you could fit everything in native int types and having
2^n + 1 integer coefficients may perform better than having 2^n rational
coefficients.  I suspect the applicability of this is rather narrow,
however, because once you get past a small handful of radicals, the
coefficients (esp. intermediate coefficients computed during
reciprocation) will easily overflow native machine int types, thus
necessitating BigInt coefficients pretty quickly.

The last time I attempted an implementation with 3-4 separate radicals
many years ago, I found that even small starting coefficients (i.e., 1-2
digits) quickly exploded in internal algorithms due to repeated
multiplication, so that after just a small number of operations I was
already running into integer overflows. This was back when I was still
doing it in C/C++... I did attempt a re-implementation using libgmp, but
never finished.


T

-- 
"Real programmers can write assembly code in any language. :-)" -- Larry Wall
Apr 20 2017
parent reply Timon Gehr <timon.gehr gmx.ch> writes:
On 20.04.2017 20:29, H. S. Teoh via Digitalmars-d wrote:

 On Thu, Apr 20, 2017 at 02:51:12PM +0200, Timon Gehr via Digitalmars-d wrote:

 On 20.04.2017 03:00, H. S. Teoh via Digitalmars-d wrote:

 On Thu, Apr 20, 2017 at 02:01:20AM +0200, Timon Gehr via Digitalmars-d wrote:
 [...]

 Yes, there is in fact a beautifully simple way to do better. :)
 ...


 Ahh, so *that's* what it's all about. I figured that's what I was
 missing. :-D  Thanks, I'll include this in QRat soon.




 Update: QRat now supports ^^. :-) Integral exponents only, of course.  I
 also implemented negative exponents, because QRat supports division and
 the same algorithm can be easily reused for that purpose.
 ...


Nice! :)


 ...

 [...]

 BTW, you are right that with std.bigint, computation using a linear
 number of additions is actually faster for my example (100000th
 Fibonacci number).  The asymptotic running time of the version with
 pow on QRats is lower though, so there ought to be a crossover point.
 (It is Θ(n^2) vs.  O(n^log₂(3)·log(n)). std.bigint does not implement
 anything that is asymptotically faster than Karatsuba.)


 Yeah, probably there is a crossover point. But it might be quite large.
 I suppose one could make a graph of the running times for increasing n,
 and either find the crossover point that way or extrapolate using the
 known curve shapes.

 Having said that, I haven't scrutinized the performance characteristics
 of QRat too carefully just yet -- there is probably room for
 optimization.


Gcd is the problem. The following code which implements a strategy based 
on matrix multiplication instead of QRat multiplication is significantly 
faster than naive linear computation:

BigInt fib(long n){
     BigInt[2] 
a=[BigInt(0),BigInt(1)],b=[BigInt(1),BigInt(2)],c=[BigInt(1),BigInt(-1)];
      for(;n;n>>=1){
	    foreach(i;1-n&1..2){
             auto d=a[i]*a[1];
             a[i]=a[i]*b[1]+c[i]*a[1];
             b[i]=b[i]*b[1]-d;
             c[i]=c[i]*c[1]-d;
         }
     }
     return a[0];
}

If I change the gcd computation in QRat (line 233) from

auto g = gcd(abs(a), abs(b), c);

to

auto g = gcd(abs(a), c, abs(b));

I get performance that is a lot closer to the matrix version and also 
beats the linear computation. (This is because if one of the operands is 
1, gcd is cheap to compute.)



 ...

 Yes, at most, except you don't need "+1". (For each radical ri, you
 will at most need to pick a power between 0 to deg(ri)-1 to index into
 the coefficients.)


 [...]

 The +1 is for the denominator, assuming integer coefficients.  Since
 having 2^n rational coefficients is equivalent to having 2^n integer
 coefficients (which are half the size of rational coefficients,
 computer-representation-wise) + 1 denominator.
 ...


Ah, I see. Personally, I'm more in the one denominator per coefficient 
camp. :) I think having a designated ℚ type is cleaner, and it might 
even be more performant.
Apr 20 2017
parent reply Timon Gehr <timon.gehr gmx.ch> writes:
On 20.04.2017 21:11, Timon Gehr wrote:

 Update: QRat now supports ^^. :-) Integral exponents only, of course.  I
 also implemented negative exponents, because QRat supports division and
 the same algorithm can be easily reused for that purpose.
 ...


 Nice! :)


It does not work with BigInt-based QRats (T(1) does not work, as 1 does 
not implicitly convert to BigInt.)
Apr 20 2017
parent Timon Gehr <timon.gehr gmx.ch> writes:
On 20.04.2017 21:18, Timon Gehr wrote:

 On 20.04.2017 21:11, Timon Gehr wrote:

 Update: QRat now supports ^^. :-) Integral exponents only, of course.  I
 also implemented negative exponents, because QRat supports division and
 the same algorithm can be easily reused for that purpose.
 ...


 Nice! :)


 It does not work with BigInt-based QRats (T(1) does not work, as 1 does
 not implicitly convert to BigInt.)


I guess the best fix is to templatize the QRat constructor such that it 
accepts all argument types that can be used to construct the coefficients.
Apr 20 2017
prev sibling parent "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
On Thu, Apr 20, 2017 at 09:11:56PM +0200, Timon Gehr via Digitalmars-d wrote:

 On 20.04.2017 20:29, H. S. Teoh via Digitalmars-d wrote:


[...]

 Having said that, I haven't scrutinized the performance
 characteristics of QRat too carefully just yet -- there is probably
 room for optimization.
 


 
 Gcd is the problem. [...]


[...]

 If I change the gcd computation in QRat (line 233) from
 
 auto g = gcd(abs(a), abs(b), c);
 
 to
 
 auto g = gcd(abs(a), c, abs(b));
 
 I get performance that is a lot closer to the matrix version and also beats
 the linear computation. (This is because if one of the operands is 1, gcd is
 cheap to compute.)


Fixed, thanks!


[...]

 The +1 is for the denominator, assuming integer coefficients.  Since
 having 2^n rational coefficients is equivalent to having 2^n integer
 coefficients (which are half the size of rational coefficients,
 computer-representation-wise) + 1 denominator.  ...


 
 Ah, I see. Personally, I'm more in the one denominator per coefficient
 camp.  :) I think having a designated ℚ type is cleaner, and it might
 even be more performant.


It's hard to say without profiling it, I think.  Having one denominator
per coefficient does need more storage space, and you do have to
separately reduce each rational coefficient to lowest terms per
operation. I think some profiling is needed to know for sure.

Besides, Phobos is badly in need of a Rational type that's compatible
with both native int types and BigInt. Maybe I should adapt some of the
code in QRat for that purpose. :-D  (Since rationals, after all, are a
subset of QRats.)


On Thu, Apr 20, 2017 at 09:25:19PM +0200, Timon Gehr via Digitalmars-d wrote:

 On 20.04.2017 21:18, Timon Gehr wrote:

 On 20.04.2017 21:11, Timon Gehr wrote:

 Update: QRat now supports ^^. :-) Integral exponents only, of
 course.  I also implemented negative exponents, because QRat
 supports division and the same algorithm can be easily reused
 for that purpose.  ...


 
 Nice! :)


 
 It does not work with BigInt-based QRats (T(1) does not work, as 1
 does not implicitly convert to BigInt.)


 
 I guess the best fix is to templatize the QRat constructor such that
 it accepts all argument types that can be used to construct the
 coefficients.


Fixed, thanks!


T

-- 
If the comments and the code disagree, it's likely that *both* are wrong. --
Christopher
Apr 20 2017