I want to find out the first n that satisfies . It's a simple and easy stuff if I use another language, like c/c++, but I don't know how to implement it in Haskell.
#include <iostream>
long double term(int k) { return 1.0/(k*k+2.0*k); }
int main() {
long double total = 0.0;
for (int k=1;;k++) {
total += term(k);
if (total>=2.99/4.0) {
std::cout << k << std::endl;
break;
}
}
return 0;
}
I used dropWhile with an ordered list and take 1 to pick up the first one.
term k = 1.0/(k*k+2.0*k)
termSum n = sum $ take n $ map term [1..]
main = do
let [(n,val)] = take 1 $ dropWhile (\(a,b)->b <= 2.99/4.0) $ map (\n->(n,termSum n)) [1..]
print n
I know it's horrible. What is the best and intuitive way to write this?
Re: Thank you for the great answers! The one using fix function seems to be the fastest in my machine (Redhat 6.4 64bit / 80GB memory)
method#0 take 1 and dropWhile (my initial implementation)
threshold=0.74999 n=99999 time=52.167 sec
method#1 using fix function
threshold=0.74999 n=99999 time=0.005 sec
threshold=0.74999999 n=101554197 time=1.077 sec
threshold=0.7499999936263 n=134217004 time=1.407 sec
method#2 working backwards
threshold=0.74999 n=99999 time=0.026 sec
threshold=0.74999999 n=101554197 time=21.523 sec
threshold=0.7499999936263 n=134217004 time=25.247 sec
method#3 imperative way
threshold=0.74999 n=99999 time=0.008 sec
threshold=0.74999999 n=101554197 time=2.460 sec
threshold=0.7499999936263 n=134217004 time=3.254 sec
ReRe: I noticed that whatever implementation I used (fix, imperative way, or recursive way), if the threshold is larger than 0.7499999936264... it never ends.. in order for f(n) to be larger than 0.7499999936264, I thought we just needed to compute the terms up to 150,000,000 since ![f(n)=\frac_{3n^2+5n}^{4n^2+12n+8}]. I used Integer instead of Int, but it did not help either. Is there any reason why it does not finish if I set the threshold larger than 0.7499999936264 ...?
I like to work backwards in these sort of situations:
main = print k where
k = 1 + length (takeWhile (< (2.99/4)) partialSums)
partialSums = scanl1 (+) terms
terms = [ 1.0/(k*k+2.0*k) | k <- [1..] ]
How this works:
terms
is an infinite list, but since Haskell is lazy, we'll only calculate as much of each term as we demand:
λ terms = [ 1.0/(k*k+2.0*k) | k <- [1..] ] :: [Double]
λ take 5 terms
[0.3333333333333333,0.125,6.666666666666667e-2,4.1666666666666664e-2,2.857142857142857e-2]
λ :p terms
terms = 0.3333333333333333 : 0.125 : 6.666666666666667e-2 :
4.1666666666666664e-2 : 2.857142857142857e-2 : (_t5::[Double])
partialSums
is another infinite list, based on the contents of terms
(using scanl1
). It lets us amortize the work that you do computing termSum
:
λ partialSums = scanl1 (+) terms
λ take 5 partialSums
[0.3333333333333333,0.4583333333333333,0.525,0.5666666666666667,0.5952380952380952]
The takeWhile (< (2.99/4))
then determines how many terms of partialSums
we need to generate and thereby how many terms of
terms
we need to generate:
λ length (takeWhile (< (2.99/4)) partialSums)
398
If we check, we can see that the sum of the first 398 terms
are less than 2.99 / 4
, but the 399th bumps it over:
λ sum (take 398 terms) < 2.99/4
True
λ sum (take 399 terms) < 2.99/4
False
Or, equivalently that the 397th partial sum (0-based index) is less than the target and that the 398th is not:
λ partialSums !! 397 < 2.99/4
True
λ partialSums !! 398 < 2.99/4
False
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With