fixedprec 0.2.1.1 → 0.2.2.0

raw patch · 3 files changed

+44/−2 lines, 3 files

Files

• ChangeLog +3/−0
• Data/Number/FixedPrec.hs +40/−1
• fixedprec.cabal +1/−1

@@ -1,5 +1,8 @@ ChangeLog +v0.2.2.0 2014/03/04+	(2014/03/04) PS1 - added high-precision quadratic equation solver.+ v0.2.1.1 2014/02/05 	(2014/02/04) PS1 - improved efficiency of square root. 	(2014/02/04) PS1 - added profiling options.

@@ -24,6 +24,7 @@      -- * Other operations   fractional,+  solve_quadratic,   log_double   ) where @@ -88,6 +89,28 @@       a = 2^(b `div` 2)       b = hibit (fromIntegral n) +-- | Find the ceiling of the larger solution of a quadratic+-- equation. Specifically, given the polynomial /p/(/x/) = /x/² + /bx/+-- + /c/, where /b/ and /c/ are integers, find the smallest integer+-- /x/ ≥ -/b/\/2 satisfying /p/(/x/) ≥ 0, if /b/^2 - 4/c/ ≥ 0.+-- +-- This is done using integer arithmetic, so there are no rounding+-- errors. It generalizes 'intsqrt'.+intquad :: Integer -> Integer -> Maybe Integer+intquad b c+  | b^2 - 4*c < 0 = Nothing+  | x1^2 + b*x1 + c >= 0 = Just x1+  | otherwise = Just (iterate x0)+  where+    iterate x+      | px <= 0 && px + 2*x+1+b > 0 = x+1+      | otherwise = iterate ((x^2 - c) `div` (2*x + b))+      where+        px = x^2 + b*x + c+    x1 = -(b `div` 2)+    x0 = x1 + 2^(h `div` 2)+    h = hibit (b^2 - 4*c)+ -- ---------------------------------------------------------------------- -- ** Other general-purpose functions             @@ -100,7 +123,7 @@ -- In other words, let /n/ = ⌊log[sub /b/] /x/⌋ and  -- /r/ = /x/ /b/[sup −/n/]. This can be more efficient than 'floor' -- ('logBase' /b/ /x/) depending on the type; moreover, it also works--- for exact types such as 'Rational' and 'EReal'.+-- for exact types such as 'Rational' and 'QRootTwo'. floorlog :: (Fractional b, Ord b) => b -> b -> (Integer, b) floorlog b x      | x <= 0            = error "floorlog: argument not positive"@@ -274,6 +297,22 @@ fractional a@(F x) = F (x `mod` one) where   p = getprec a   one = (decshiftL p 1)++-- | Solve the quadratic equation /x/^2 + /bx/ + /c/ = 0 with maximal+-- possible precision, using a numerically stable method. Return the+-- pair (/x1/, /x2/) of solutions with /x1/ <= /x2/, or 'Nothing' if no+-- solution exists.+-- +-- This is far more precise, and probably more efficient, than naively+-- using the quadratic formula.+solve_quadratic :: (Precision e) => FixedPrec e -> FixedPrec e -> Maybe (FixedPrec e, FixedPrec e)+solve_quadratic b c = do+  let p = getprec b + 3+      b' = floor (b * 10^p)+      c' = floor (c * 100^p)+  x2' <- intquad b' c'+  let x1' = -b' - x2'+  return (fromInteger x1' / 10^p, fromInteger x2' / 10^p)  -- ---------------------------------------------------------------------- -- ** Power series

@@ -1,6 +1,6 @@ name:           fixedprec -- Don't forget to update the ChangeLog-version:        0.2.1.1+version:        0.2.2.0 license:        BSD3 cabal-version:  >= 1.8 build-type:	Simple