Problems with truncate big float values in Haskell
Hi, Firstly I desire one year brilliant for all on group. I have one simple problem but not solved for me for some mistake mine. I have one Rational number (3 + sqrt 5) and I need get the last three digits from integral part of that number. 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). I have attempted that on ghci for various values and get same result with different powers: let x = (truncate ((3 + sqrt 5) ^ 2000000)) `mod` 1000 and result is 216. let y = (truncate ((3 + sqrt 5) ^ 20000000)) `mod` 1000 and result is 216. let w = (truncate ((3 + sqrt 5) ^ 200000000)) `mod` 1000 and result is 216. let z = (truncate ((3 + sqrt 5) ^ 2000000000)) `mod` 1000 and result is 216. The result is same for all: 216, and is incorrect because any that expressions have one exclusive value. What I do for complete that task correctly and what I'm mistaken here?! Thank you, Josenildo Silva
On 2017-01-05 01:16 PM, Henson wrote:
I have one Rational number (3 + sqrt 5)
Irrational Number.
let x = (truncate ((3 + sqrt 5) ^ 2000000)) `mod` 1000 and result is 216.
You need to firstly read http://floating-point-gui.de/formats/fp/ , and secondly google for "goldberg floating point" and read "what every computer scientist should know about floating-point arithmetic".
On 2017-01-05 02:30 PM, Albert Y. C. Lai wrote:
On 2017-01-05 01:16 PM, Henson wrote:
I have one Rational number (3 + sqrt 5)
Irrational Number.
let x = (truncate ((3 + sqrt 5) ^ 2000000)) `mod` 1000 and result is 216.
You need to firstly read http://floating-point-gui.de/formats/fp/ , and secondly google for "goldberg floating point" and read "what every computer scientist should know about floating-point arithmetic".
And the conclusion at that point should be: A direct Double approach will not work. ("direct" means you don't do your own math first.) You will need to be either indirect (do your own math first) or forget Double. In fact I would do both indirect and forget Double. Let r = 3 + sqrt 5, s = 3 - sqrt 5. s is the evil twin of r. Define (as in math, not code) h n = r^n + s^n for natural n. (a) Prove: r and s are the roots of x^2 - 6x + 4 = 0. And so r^2 = 6r - 4, s^2 = 6s - 4. (How did I think up that equation? From (x - r)*(x - s).) (b) Prove: h satisfies the following recurrence h 0 = 2 h 1 = 6 h n = 6 * h (n-1) - 4 * h (n-2) for n>=2 (Sketch: Part (a) will be useful.) (c) Prove: floor(r^n) = h n - 1. (Sketch: 0 < s^n < 1, and the recurrence tells you that h n is a natural number.) (d) Give and/or code up an algorithm to compute the last three decimal digits of floor(r^n); it should use only O(n) time, O(lg n) space, and integer/natural arithmetic. In fact, the lg n space is only for counting from 0 to n, the algorithm is constant-space otherwise (for mod-1000 arithmetic). (Actually you could get it down to O(lg n) time too, but that's for another day.) (e) Do a similar thing to floor((the gold ratio)^n). How I came up with this idea: I recalled that the fibonacci numbers are related to (the golden ratio)^n (for example CLRS has a few lines on this). So I just had to adapt it to 3 + sqrt 5. * * * For a decade, I haven't had a good excuse to use my smug saying. So this is a good chance to say: The language of computing is not C, Node.js, or Haskell. The language of computing is mathematics.
On Thu, Jan 5, 2017 at 2:18 PM Henson
I have one Rational number (3 + sqrt 5) and I need get the last three digits from integral part of that number.
But that number is irrational, not rational. The Haskell sqrt function from the Prelude is a method of the Floating class, so it works by default on the Double type, which is a floating point number type with a fixed-size finite representation that approximates real numbers with varying precision over a limited range.
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).
I have attempted that on ghci for various values and get same result with different powers:
let x = (truncate ((3 + sqrt 5) ^ 2000000)) `mod` 1000 and result is 216.
let y = (truncate ((3 + sqrt 5) ^ 20000000)) `mod` 1000 and result is 216.
let w = (truncate ((3 + sqrt 5) ^ 200000000)) `mod` 1000 and result is 216.
let z = (truncate ((3 + sqrt 5) ^ 2000000000)) `mod` 1000 and result is 216.
The result is same for all:
216, and is incorrect because any that expressions have one exclusive value.
The Double type has a special value, Infinity, to represent (roughly) numbers too large to represent otherwise. You can get at it with a zero division: λ> 1/0 :: Double Infinity You can also get it with computations resulting in numbers too large to represent otherwise: λ> 10^308 :: Double 1.0000000000000006e308 λ> 10^309 :: Double Infinity The truncate method of the RealFrac class applied to this Infinity value can convert it to the smallest Integer value that is represented as Infinity when converted to Double: λ> truncate (1/0 :: Double) 179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137216 Indeed, the last three digits of the decimal representation of that number are 216. The immediately preceding Integer is indeed converted to a non-Infinity Double value: λ> fromInteger (pred (truncate (1/0 :: Double))) :: Double 1.7976931348623157e308 I believe this is all expected behavior. This article is very helpful and provides detailed explanations of many related topics: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html --- While playing around with this, I discovered a curious difference between Float and Double. The smallest Integer that does not yield Infinity on conversion to Float is not anywhere close to the value of the expression «truncate (1/0 :: Float) :: Integer»: λ> let biggestFloat :: Integer; biggestFloat = 340282356779733661637539395458142568447 λ> fromInteger biggestFloat :: Float 3.4028235e38 λ> fromInteger (succ biggestFloat) :: Float Infinity λ> truncate (1/0 :: Float) 340282366920938463463374607431768211456 λ> truncate (1/0 :: Float) - biggestFloat 10141204801825835211973625643009 With Double, «truncate (1/0 :: Double)» is indeed the immediate successor of the largest Integer not converted to Infinity: λ> let biggestDouble :: Integer; biggestDouble = 179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137215 λ> fromInteger biggestDouble :: Double 1.7976931348623157e308 λ> fromInteger (succ biggestDouble) :: Double Infinity λ> truncate (1/0 :: Double) - biggestDouble 1
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-Co... it seems that everything you need is there.
participants (4)
-
Albert Y. C. Lai -
Henson -
Manuel Gómez -
Richard A. O'Keefe