[Haskell-cafe] Problems with truncate big float values in Haskell
Richard A. O'Keefe
ok at cs.otago.ac.nz
Fri Jan 6 01:22:56 UTC 2017
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.
More information about the Haskell-Cafe
mailing list