Exact arithmetic with quadratic irrationals

[Timon Gehr via Digitalmars-d]

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.