[Haskell-cafe] Problems with truncate big float values in Haskell

[Richard A. O'Keefe]

On 6/01/17 7:16 AM, Henson wrote:
> I have one Rational number (3 + sqrt 5) and I need get the last three
> digits from integral part of that number.

3 + sqrt 5 is not a rational number.

Prelude> :type 3 + sqrt 5
3 + sqrt 5 :: Floating a => a

So it's not a value of type Rational either.

What's more, the value is approximately
5.236067977499789696409173668731276235440618359611525724270897
(thanks, bc(1)) so it doesn't HAVE three digits in its integral
part.

>
> The question is the number is arbitrary and that number always is power
> of n and n is between 2(Two) and 2000000000(Two billions).

Ah.  So what you are really asking about is the integral part of
(3 + sqrt 5)^n, which clearly lies somewhere strictly between 5^n
and 6^n.  If n can be 2x10**9, those are going to be very big
numbers.

It is useful to understand just what kind of number you have here.
You have (a + b * sqrt p) where p is a prime.
That is a quadratic surd
http://mathworld.wolfram.com/QuadraticSurd.html
or quadratic irrational number
https://en.wikipedia.org/wiki/Quadratic_irrational_number

Numbers (a + b sqrt c)/d where a and b are integers, d is a positive
integer, and c is a positive square-free number, are
real quadratic irrationals.  The set of such numbers with the
same value of c forms a field.

For fixed c, let a b e f be rationals.

(a + b sqrt c) ± (e + f sqrt c) = (a ± e) + (b ± f) sqrt c

(a + b sqrt c) × (e + f sqrt c) =
   (a × e + b × f × c) + (a × f + b × e) sqrt c

If a, b, e, f are integers, the sum difference and product
also have integer coefficients, which is nice.

So you can do

data S = S Integer Integer
c :: Integer
c = 5

unit :: S
unit = S 1 0

times :: S -> S -> S
times (S a b) (S e f) = S (a*e+b*f*c) (a*f+b*e)

power :: S -> Int -> S
power x n =
   case compare n 0 of
     LT -> error "negative exponent"
     EQ -> unit
     GT -> f x (n-1) x
           where f _ 0 y = y
                 f x n y = g x n
                 g x n y = if even n then g (times x x) (n`quot`2) y
                           else f x (n-1) (times x y)

This lets you compute power (S 3 1) n *precisely* for any
non-negative Int n.  What it doesn't let you do is get the
integer part.  Using a good enough approximation to sqrt 5,
I think the answers corresponding to n = 0..20 are
[1,5,27,143,751,935,607,903,991,335,47,943,471,55,447,463,991,95,607,263,152]

But there is another approach.
A number is a quadratic irrational if and only if it is
a regular periodic continued fraction.
sqrt 5 = 2 + 1/(4 + 1/(4 + 1/(4 + ...)))
so 3 + sqrt 5 = 5 + 1/(4 + 1/(4 + 1/(4 + ...)))
-- which you will note *is* a regular periodic continued fraction --
so now all you have to do is find out how to multiply
continued fractions.

Haskell being lazy makes it a nice choice for handling continued
fractions, but since these are *periodic* continued fractions we
could do this even in a strict language.

A web search will quickly turn up things like
http://perl.plover.com/yak/cftalk/INFO/gosper.txt
which admittedly isn't easy reading if you don't already
understand continued fractions.

On the other hand, *another* web search turns up
https://hackage.haskell.org/package/continued-fractions
https://hackage.haskell.org/package/cf
and looking for example at
https://hackage.haskell.org/package/continued-fractions-0.9.1.1/docs/Math-ContinuedFraction.html
it seems that everything you need is there.