Re: [Haskell-cafe] Fractional sqrt
Henning Thielemann wrote: there is already an implementation of continued fractions for approximation of roots and transcendent functions by Jan Skibinski: http://darcs.haskell.org/numeric-quest/Fraction.hs I wrote: Wow, nice. Now - how was I supposed to have found that? http://www.haskell.org/haskellwiki/Libraries_and_tools/ It turns out that it is mentioned on the wiki, on the Mathematics and Physics subpage of Libraries and Tools. http://www.haskell.org/haskellwiki/Libraries_and_tools/Mathematics But the references were a bit unclear and out of date. I added a new page for the Numeric Quest library, and added it to the Mathematics category. (Should it be in any other category?) I also updated the references to Numeric Quest on the Mathematics and Physics page. Can someone with access to darcs.haskell.org please fix this library? darcs get currently does not seem to work for it. http://darcs.haskell.org/numeric-quest/ Also, these modules should each be moved into the hierarchy somewhere. Does anyone have suggestions? It may help to look at: http://www.haskell.org/haskellwiki/Numeric_Quest Thanks, Yitz ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
On Mon, 22 Jan 2007, Yitzchak Gale wrote: Henning Thielemann wrote: there is already an implementation of continued fractions for approximation of roots and transcendent functions by Jan Skibinski: http://darcs.haskell.org/numeric-quest/Fraction.hs I wrote: Wow, nice. Now - how was I supposed to have found that? Simon Peyton-Jones wrote: Perhaps one of you can add this to the libraries-and-tools page? http://www.haskell.org/haskellwiki/Libraries_and_tools/ I added a new page for the Numeric Quest library, and added it to the Mathematics category. (Should it be in any other category?) I also updated the references to Numeric Quest on the Mathematics and Physics page. Great work! However I would not throw away the links to the Internet archive because the pages stored there link to the original HTML pages which are lost. It's cumbersome to navigate through the archived pages. The Darcs repository at cvs.haskell.org is no official update by Jan Skibinski - I set up that repository to simplify further development. I added the keyword continued fractions to the Wiki page so the module can be discovered by other people. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
On Mon, Jan 22, 2007 at 03:26:38PM +0200, Yitzchak Gale wrote: Can someone with access to darcs.haskell.org please fix this library? darcs get currently does not seem to work for it. http://darcs.haskell.org/numeric-quest/ I've fixed the permissions, although applying patches in the future might break them again. Thanks Ian ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
I wrote: I added a new page for the Numeric Quest library... http://www.haskell.org/haskellwiki/Numeric_Quest I also updated the references to Numeric Quest on the Mathematics and Physics page. http://www.haskell.org/haskellwiki/Libraries_and_tools/Mathematics Henning Thielemann wrote I would not throw away the links to the Internet archive For modules that were previously mentioned on the Mathematics page and not re-hosted in your darcs repo, I left the links as they were. For modules that are in the repo, the original version can always be obtained from the repo. In addition, there is a link on the Numeric Quest page to the author's main Haskell page in the Internet Archive. From there, it is easy to navigate with one click to any of the author's original work still available in the Archive, including all of the modules that are also in the darcs repo. Perhaps I should add a note to that effect. I intentionally removed the direct links to the Archive from the Mathematics page for modules that are now in your repo. I estimate that people clicking from there are more likely to be lead astray than to be helped. Do you think it is important to have a separate link for each repo module to the Internet Archive? If so, they can be added to the Numeric Quest page. The Darcs repository at cvs.haskell.org is no official update by Jan Skibinski - I set up that repository to simplify further development. OK - I'll add a note to that effect. I added the keyword continued fractions to the Wiki page so the module can be discovered by other people. Good, thanks. -Yitz ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
On Sun, 21 Jan 2007, Yitzchak Gale wrote: You always get the best approximation for a Rational using continued fractions. For most calculations that is easier said than done, but for square root we are very lucky. Certainly no surprise - there is already an implementation of continued fractions for approximation of roots and transcendent functions by Jan Skibinski: http://darcs.haskell.org/numeric-quest/Fraction.hs ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
Henning Thielemann wrote: Certainly no surprise - there is already an implementation of continued fractions for approximation of roots and transcendent functions by Jan Skibinski: http://darcs.haskell.org/numeric-quest/Fraction.hs Wow, nice. Now - how was I supposed to have found that? It has been know at least since early Knuth that continued fractions are the best approximation to rationals (in a certain sense). It is too bad that Data.Ratio took the idea of simplest fraction from Scheme for the approxRational function. I wrote my own replacement using continued fractions. Now I see that I have reinvented several wheels. Thanks, Yitz ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
G'day all. I said: I've also extended the range for approxSmallSqrt here from (0,255) to (0,271). It is left as an exercise as to why this might be a good idea. (Hint: 272 is approximately 16.5*16.5.) The correct answer, for those playing at home, is it's because it WAS a good idea when I wrote the code for integer square roots. Mindlessly cutting and pasting it into code for fractions makes it not such a good idea. Cheers, Andrew Bromage ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
Sorry folks, it is just wrong to use Newton's method
for Rational.
Andrew Bromage wrote:
First off, note that for fractions, sqrt(p/q) = sqrt p / sqrt q.
Don't do that for Rational - you lose precious precision.
The whole idea for calculations in Rational is to find
a lucky denominator that is not too big but gives you
a good approximation. Newton's method will converge
very quickly - and it will also blow up your stack
very quickly as the denominator increases exponentially.
You always get the best approximation for a Rational
using continued fractions. For most calculations that
is easier said than done, but for square root we are
very lucky.
Is the continued fraction algorithm fast enough?
Well, on my old machine, the following code computes
sqrt 2 to within 10^(-1000) with no noticeable delay.
This can obviously be optimized in many ways,
but I am not sure it is worth the effort.
Regards,
Yitz
\begin{code}
-- Arbitrary precision square root for Rational
module SquareRoot where
import Data.Ratio
-- The square root of x, within error tolerance e
-- Warning: diverges if e=0 and x is not a perfect square
rSqrt :: Rational - Rational - Rational
rSqrt x _ | x 0 = error rSqrt of negative
rSqrt _ e | e 0 = error rSqrt negative tolerance
rSqrt x e = firstWithinE $ map (uncurry (%)) $ convergents $
sqrtCF (numerator x) (denominator x)
where
firstWithinE = head . dropWhile off
off y = abs (x - y^2 - e^2) 2 * e * y
-- The continued fraction for the square root of positive p%q
sqrtCF :: Integer - Integer - [Integer]
sqrtCF p q = cf 1 0 1
where
-- Continued fraction for (a*sqrt(p%q)+b)/c
cf a b c
| c == 0= []
| otherwise = x : cf (a' `div` g) (b' `div` g) (c' `div` g)
where
x = (iSqrt (p*a*a) q + b) `div` c
a' = a * c * q
e = c * x - b
b' = c * e * q
c' = p*a*a - q*e*e
g = foldl1 gcd [a', b', c']
-- The greatest integer less than or equal to sqrt (p%q) for positive p%q
iSqrt :: Integer - Integer - Integer
iSqrt p q = sq 0 (1 + p `div` q)
where
sq y z
| z - y 2 = y
| q * m * m = p = sq m z
| otherwise = sq y m
where m = (y + z) `div` 2
-- The following could go in a separate module
-- for general continued fractions
-- The convergents of a continued fraction
convergents :: [Integer] - [(Integer, Integer)]
convergents cf = drop 2 $ zip h k
where
h = 0 : 1 : zipWith3 affine cf (drop 1 h) h
k = 1 : 0 : zipWith3 affine cf (drop 1 k) k
affine x y z = x * y + z
/end{code}
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
On Thu, 18 Jan 2007, Novák Zoltán wrote: I have to admit that I do not fully understand the Haskell numerical tower... Now I'm using the Newton method: mysqrt :: (Fractional a) = a - a mysqrt x = (iterate (\y - (x / y + y) / 2.0 ) 1.0) !!2000 But I want a faster solution. (Not by replacing 2000 with 100... :) Newton method for sqrt is very fast. It converges quadratically, that is in each iteration the number of correct digits doubles. The problem is to find a good starting approximation. What about comparing powers of two and their squares against 'x' and use an appropriate power of two as starting value? Sure, this will require Real instance and won't work on complex numbers. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
Hi, Thanks for the answers. The best solution would be a general purpose arbitrary precision math library for Haskell. I found two: http://medialab.freaknet.org/bignum/ http://r6.ca/FewDigits/ I think both uses power series and have trigonometric functions too. (I only need sqrt, so probably I will not use them, just the simple Newton alg.) It would be great, if a fast implementation of arbitrary precision real arithmetic Class would be part of the Standard libraries. Because Haskell is so great for numerical computations already (functions, composing functions, higher order f. with typechecking), a good math library would make Matlab obsolete. :) Zoltan Novak ps. Call me Zoltan (that is my first name), just we write it in different order in Hungary. (Names are in hungarian notation here: the type description precedes the actual name) ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
Hi Zoltán, I only need sqrt, so probably I will... use... just the simple Newton alg. It is still not clear to me what type you want to work in. Is it Rational? In that case, you don't need the Newton algorithm. realToFrac . sqrt . realToFrac works fine, as you originally suggested. If that gives too much precision, you can use (`approxRational` epsilon) . sqrt . realToFrac Newton will not work in general for Fractional. It works for Floating, but there you already have the sqrt function. Are you using a type that has a Fractional instance, does not have a Floating instance, and Newton works and is needed? What type is that? a good math library would make Matlab obsolete. :) Yes, that would certainly be nice. This has been discussed at length several times on this list in the past. I wonder if there has been any further progress on it. Regards, Yitz ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
If it's arbitrary precision floating point that you want then sqrt should where it already is, as a member of Floating. (I find arbitrary precision real to be an oxymoron, the real numbers are the real numbers, they already have arbitrary precision.) For a real number module, you can use, e.g., David Lester's implementation, http://www.augustsson.net/Darcs/CReal/ -- Lennart On Jan 19, 2007, at 07:13 , Novák Zoltán wrote: Hi, Thanks for the answers. The best solution would be a general purpose arbitrary precision math library for Haskell. I found two: http://medialab.freaknet.org/bignum/ http://r6.ca/FewDigits/ I think both uses power series and have trigonometric functions too. (I only need sqrt, so probably I will not use them, just the simple Newton alg.) It would be great, if a fast implementation of arbitrary precision real arithmetic Class would be part of the Standard libraries. Because Haskell is so great for numerical computations already (functions, composing functions, higher order f. with typechecking), a good math library would make Matlab obsolete. :) Zoltan Novak ps. Call me Zoltan (that is my first name), just we write it in different order in Hungary. (Names are in hungarian notation here: the type description precedes the actual name) ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
G'day all.
Quoting Henning Thielemann [EMAIL PROTECTED]:
Newton method for sqrt is very fast. It converges quadratically, that is
in each iteration the number of correct digits doubles. The problem is to
find a good starting approximation.
Yup. So how might we go about doing this?
First off, note that for fractions, sqrt(p/q) = sqrt p / sqrt q.
Secondly, note that for floating point numbers,
sqrt(m * b^e) = sqrt m * b^(e/2).
So an approach to get an initial approximation to p/q might be:
1. Compute approximate square roots for p and q, then divide.
2. For each of p and q, express as a floating-point number
m * b^(2e), compute an approximate square root for m (which will
be less than 2b), then multiply by b^e.
We've now reduced the range of approximate square roots
that we want to numbers in the range 0 to 2b, where b is a parameter
that you can use to tune the algorithm.
OK, so let's pick b=16 out of the air. Our code might look something
like this:
sqrtApprox :: Rational - Rational
sqrtApprox x = sqrtApprox' (numerator x) / sqrtApprox' (denominator x)
sqrtApprox' :: Integer - Rational
sqrtApprox' n
| n 0 = error sqrtApprox'
| otherwise = approx n 1
where
approx n acc
| n 256 = (acc%1) * approxSmallSqrt (fromIntegral n)
| otherwise = approx (n `div` 256) (acc*16)
approxSmallSqrt :: Int - Rational
approxSmallSqrt n = {- detail omitted -}
The approxSmallSqrt function is a good candidate for the memoising
CAFs approach:
http://haskell.org/hawiki/MemoisingCafs
You can populate the memo table by using a fixed number of iterations
of the Newton-Raphson algorithm on a suitably sane initial estimate.
You can get a slightly better estimate by being a bit more careful about
rounding. Suppose that n = 2 + 256 * n', then a good estimate for
sqrt n is 16 * sqrt n'. But if n = 255 + 256 * n', the algorithm above
would also estimate sqrt n as 16 * sqrt n'. It's only slightly more
expensive (and occasionally better) to round up instead, and use
16 * sqrt (n'+1) as the estimate.
So you might want to do this instead:
sqrtApprox' :: Integer - Rational
sqrtApprox' n
| n 0 = error sqrtApprox'
| otherwise = approx n 1
where
approx n acc
| n 272
= (acc%1) * approxSmallSqrt (fromIntegral n)
| r 16
= approx q (acc*16)
| otherwise
= approx (q+1) (acc*16)
where (q,r) = n `divMod` 256
approxSmallSqrt :: Int - Rational
approxSmallSqrt n = {- detail omitted -}
I've also extended the range for approxSmallSqrt here from (0,255) to
(0,271). It is left as an exercise as to why this might be a good idea.
(Hint: 272 is approximately 16.5*16.5.)
Cheers,
Andrew Bromage
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
Hello, On Friday 19 January 2007 16:48, [EMAIL PROTECTED] wrote: ... sqrtApprox' :: Integer - Rational sqrtApprox' n | n 0 = error sqrtApprox' | otherwise = approx n 1 where approx n acc | n 256 = (acc%1) * approxSmallSqrt (fromIntegral n) | otherwise = approx (n `div` 256) (acc*16) ... In the Haskell report section 14.4 (and also e.g. in GHC.Float), we find -- Compute the (floor of the) log of i in base b. -- Simplest way would be just divide i by b until it's smaller then b, but that would -- be very slow! We are just slightly more clever. integerLogBase :: Integer - Integer - Int integerLogBase b i | i b = 0 | otherwise = doDiv (i `div` (b^l)) l where -- Try squaring the base first to cut down the number of divisions. l = 2 * integerLogBase (b*b) i doDiv :: Integer - Int - Int doDiv x y | x b = y | otherwise = doDiv (x `div` b) (y+1) Something like this could probably be used to find a suitable sqrt approximation for an Integer: sqrtApproxViaLog i = 2^(integerLogBase 2 i `div`2) Best regards Thorkil ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] Fractional sqrt
Hello, I would like to use the sqrt function on Fractional numbers. (mysqrt :: (Fractional a) = a-a) Half of the problem is solved by: Prelude :t (realToFrac.sqrt) (realToFrac.sqrt) :: (Fractional b, Real a, Floating a) = a - b For the other half I tried: Prelude :t (realToFrac.sqrt.realToFrac) (realToFrac.sqrt.realToFrac) :: (Fractional b, Real a) = a - b Prelude :t (realToFrac.sqrt.fromRational.toRational) (realToFrac.sqrt.fromRational.toRational) :: (Fractional b, Real a) = a - b Prelude :t (realToFrac.sqrt.realToFrac.fromRational.toRational) (realToFrac.sqrt.realToFrac.fromRational.toRational) :: (Fractional b, Real a) = a - b I have to admit that I do not fully understand the Haskell numerical tower... Now I'm using the Newton method: mysqrt :: (Fractional a) = a - a mysqrt x = (iterate (\y - (x / y + y) / 2.0 ) 1.0) !!2000 But I want a faster solution. (Not by replacing 2000 with 100... :) Zoltan ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
I don't see a much better way than using something like Newton- Raphson and testing for some kind of convergence. The Fractional class can contain many things; for instance it contains rational numbers. So your mysqrt function would have to be able to cope with returning arbitrary precision results. As a first step you should specify what mysqrt should return when it can't return the exact result. For instance, what would you like mysqrt (2%1) to return? -- Lennart On Jan 18, 2007, at 18:15 , Novák Zoltán wrote: Hello, I would like to use the sqrt function on Fractional numbers. (mysqrt :: (Fractional a) = a-a) Half of the problem is solved by: Prelude :t (realToFrac.sqrt) (realToFrac.sqrt) :: (Fractional b, Real a, Floating a) = a - b For the other half I tried: Prelude :t (realToFrac.sqrt.realToFrac) (realToFrac.sqrt.realToFrac) :: (Fractional b, Real a) = a - b Prelude :t (realToFrac.sqrt.fromRational.toRational) (realToFrac.sqrt.fromRational.toRational) :: (Fractional b, Real a) = a - b Prelude :t (realToFrac.sqrt.realToFrac.fromRational.toRational) (realToFrac.sqrt.realToFrac.fromRational.toRational) :: (Fractional b, Real a) = a - b I have to admit that I do not fully understand the Haskell numerical tower... Now I'm using the Newton method: mysqrt :: (Fractional a) = a - a mysqrt x = (iterate (\y - (x / y + y) / 2.0 ) 1.0) !!2000 But I want a faster solution. (Not by replacing 2000 with 100... :) Zoltan ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Fractional sqrt
Lennart Augustsson wrote: I don't see a much better way than using something like Newton- Raphson and testing for some kind of convergence. The Fractional class can contain many things; for instance it contains rational numbers. So your mysqrt function would have to be able to cope with returning arbitrary precision results. As a first step you should specify what mysqrt should return when it can't return the exact result. For instance, what would you like mysqrt (2%1) to return? Fractional also contains RealFloat a = Complex a. For that, you need to sqrt the magnitude, divide the phase by 2, and pick a branch. Which branch to pick depends on how you are using the calculation. Perhaps Zoltán has some other instance of Fractional that is not in the standard libraries. One can define an instance of Fractional for a numeric approximation of a p-adic field, where taking the square root means something like recursively taking square roots in the field of p elements and collecting the results. Or for an arbitrary finite field (if you are willing to accept an occasional _|_ for fromRational). Or for algebraic function fields. Etc. I really don't think there is anything that will work in general for Fractional. You have to look at each specific type. Regards, Yitz ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
