Exact arithmetic with quadratic irrationals

[Timon Gehr via Digitalmars-d]

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