diagrams-solve 0.1.3.1 → 0.2

raw patch · 5 files changed

+185/−177 lines, 5 filesdep −deepseqdep −tasty-hunitPVP ok

version bump matches the API change (PVP)

Dependencies removed: deepseq, tasty-hunit

API changes (from Hackage documentation)

Files

• CHANGES.markdown +10/−0
• diagrams-solve.cabal +2/−4
• src/Diagrams/Solve/Polynomial.hs +101/−105
• src/Diagrams/Solve/Tridiagonal.hs +42/−47
• tests/Test.hs +30/−21

@@ -1,3 +1,13 @@+* 0.2 (2 May 2026)++  - Fixes for various pattern-match warnings+  - Remove some unnecessary dependencies (`deepseq`, `tasty-hunit`)+  - The types of `solveTriDiagonal` and `solveCyclicTridiagonal` have+    changed to take arguments of type `NonEmpty a` instead of `[a]`.+    Previously, they simply crashed when given empty lists as+    arguments.+  - Test with GHC 9.14+ * 0.1.3.1 (19 Feb 2025)    Test with up through GHC 9.12

@@ -1,5 +1,5 @@ name:                diagrams-solve-version:             0.1.3.1+version:             0.2 synopsis:            Pure Haskell solver routines used by diagrams description:         Pure Haskell solver routines used by the diagrams                      project.  Currently includes finding real roots@@ -15,7 +15,7 @@ build-type:          Simple extra-source-files:  README.markdown, CHANGES.markdown cabal-version:       >=1.10-Tested-with:         GHC ==8.0.2 || ==8.2.2 || ==8.4.4 || ==8.6.5 || ==8.8.4 || ==8.10.7 || ==9.0.2 || ==9.2.8 || ==9.4.8 || ==9.6.5 || ==9.8.2 || ==9.10.1 || ==9.12.1+Tested-with:         GHC ==8.0.2 || ==8.2.2 || ==8.4.4 || ==8.6.5 || ==8.8.4 || ==8.10.7 || ==9.0.2 || ==9.2.8 || ==9.4.8 || ==9.6.5 || ==9.8.2 || ==9.10.1 || ==9.12.1 || ==9.14.1 Source-repository head   type:     git   location: http://github.com/diagrams/diagrams-solve.git@@ -34,8 +34,6 @@   hs-source-dirs: tests   default-language:    Haskell2010   build-depends:       base >= 4.2 && < 5.0,-                       deepseq >= 1.3 && < 1.6,                        diagrams-solve,                        tasty >= 0.10 && < 1.6,-                       tasty-hunit >= 0.9.2 && < 0.11,                        tasty-quickcheck >= 0.8 && < 0.12

@@ -1,4 +1,3 @@------------------------------------------------------------------------------ -- | -- Module      :  Diagrams.Solve.Polynomial -- Copyright   :  (c) 2011-2015 diagrams-solve team (see LICENSE)@@ -6,26 +5,22 @@ -- Maintainer  :  diagrams-discuss@googlegroups.com -- -- Exact solving of low-degree (n <= 4) polynomials.----------------------------------------------------------------------------------module Diagrams.Solve.Polynomial-       ( quadForm-       , cubForm-       , quartForm-       , cubForm'-       , quartForm'-       ) where--import           Data.List (maximumBy)-import           Data.Ord  (comparing)+module Diagrams.Solve.Polynomial (+  quadForm,+  cubForm,+  quartForm,+  cubForm',+  quartForm',+) where -import           Prelude   hiding ((^))-import qualified Prelude   as P ((^))+import Data.Maybe (fromMaybe, listToMaybe)+import Prelude hiding ((^))+import qualified Prelude as P ((^))  -- | The fundamental circle constant, /i.e./ ratio between a circle's --   circumference and radius. tau :: Floating a => a-tau = 2*pi+tau = 2 * pi  -- | A specialization of (^) to Integer --   c.f. http://comments.gmane.org/gmane.comp.lang.haskell.libraries/21164@@ -33,7 +28,7 @@ --   power's Integral type... was a genuinely bad idea." - Edward Kmett -- --   Note there are rewrite rules in GHC.Real to expand small exponents.-(^) :: (Num a) => a -> Integer -> a+(^) :: Num a => a -> Integer -> a (^) = (P.^)  -- | Utility function used to avoid singularities@@ -48,35 +43,30 @@ -- | The quadratic formula. quadForm :: (Floating d, Ord d) => d -> d -> d -> [d] quadForm a b c--    -- There are infinitely many solutions in this case,-    -- so arbitrarily return 0+  -- There are infinitely many solutions in this case,+  -- so arbitrarily return 0   | a == 0 && b == 0 && c == 0 = [0]--    -- c /= 0+  -- c /= 0   | a == 0 && b == 0 = []--    -- linear-  | a == 0    = [-c/b]--    -- no real solutions-  | d < 0     = []--    -- ax^2 + c = 0-  | b == 0    = [sqrt (-c/a), -sqrt (-c/a)]--    -- multiplicity 2 solution-  | d == 0    = [-b/(2*a)]--    -- see http://www.mpi-hd.mpg.de/astrophysik/HEA/internal/Numerical_Recipes/f5-6.pdf-  | otherwise = [q/a, c/q]- where d = b^2 - 4*a*c-       q = -1/2*(b + signum b * sqrt d)+  -- linear+  | a == 0 = [-c / b]+  -- no real solutions+  | d < 0 = []+  -- ax^2 + c = 0+  | b == 0 = [sqrt (-c / a), -sqrt (-c / a)]+  -- multiplicity 2 solution+  | d == 0 = [-b / (2 * a)]+  -- see http://www.mpi-hd.mpg.de/astrophysik/HEA/internal/Numerical_Recipes/f5-6.pdf+  | otherwise = [q / a, c / q]+ where+  d = b ^ 2 - 4 * a * c+  q = -1 / 2 * (b + signum b * sqrt d) {-# INLINE quadForm #-}  _quadForm_prop :: Double -> Double -> Double -> Bool _quadForm_prop a b c = all (aboutZero' 1e-10 . eval) (quadForm a b c)-  where eval x = a*x^2 + b*x + c+ where+  eval x = a * x ^ 2 + b * x + c  ------------------------------------------------------------ -- Cubic formula@@ -88,34 +78,32 @@ --   list of all real roots. First argument is tolerance. cubForm' :: (Floating d, Ord d) => d -> d -> d -> d -> d -> [d] cubForm' toler a b c d-  | aboutZero' toler a      = quadForm b c d--    -- three real roots, use trig method to avoid complex numbers-  | delta >  0              = map trig [0,1,2]--    -- one real root of multiplicity 3-  | delta == 0 && disc == 0 = [ -b/(3*a) ]--    -- two real roots, one of multiplicity 2-  | delta == 0 && disc /= 0 = [ (b*c - 9*a*d)/(2*disc)-                              , (9*a^2*d - 4*a*b*c + b^3)/(a * disc)-                              ]--    -- one real root (and two complex)-  | otherwise               = [-b/(3*a) - cc/(3*a) + disc/(3*a*cc)]-- where delta  = 18*a*b*c*d - 4*b^3*d + b^2*c^2 - 4*a*c^3 - 27*a^2*d^2-       disc   = 3*a*c - b^2-       qq     = sqrt(-27*(a^2)*delta)-       qq'    = if abs (xx + qq) > abs (xx - qq) then qq else -qq-       cc     = cubert (1/2*(qq' + xx))-       xx     = 2*b^3 - 9*a*b*c + 27*a^2*d-       p      = disc/(3*a^2)-       q      = xx/(27*a^3)-       phi = 1/3*acos(3*q/(2*p)*sqrt(-3/p))-       trig k = 2 * sqrt(-p/3) * cos(phi - k*tau/3) - b/(3*a)-       cubert x | x < 0     = -((-x)**(1/3))-                | otherwise = x**(1/3)+  | aboutZero' toler a = quadForm b c d+  -- three real roots, use trig method to avoid complex numbers+  | delta > 0 = map trig [0, 1, 2]+  -- one real root of multiplicity 3+  | delta == 0 && disc == 0 = [-b / (3 * a)]+  -- two real roots, one of multiplicity 2+  | delta == 0 && disc /= 0 =+      [ (b * c - 9 * a * d) / (2 * disc)+      , (9 * a ^ 2 * d - 4 * a * b * c + b ^ 3) / (a * disc)+      ]+  -- one real root (and two complex)+  | otherwise = [-b / (3 * a) - cc / (3 * a) + disc / (3 * a * cc)]+ where+  delta = 18 * a * b * c * d - 4 * b ^ 3 * d + b ^ 2 * c ^ 2 - 4 * a * c ^ 3 - 27 * a ^ 2 * d ^ 2+  disc = 3 * a * c - b ^ 2+  qq = sqrt (-27 * (a ^ 2) * delta)+  qq' = if abs (xx + qq) > abs (xx - qq) then qq else -qq+  cc = cubert (1 / 2 * (qq' + xx))+  xx = 2 * b ^ 3 - 9 * a * b * c + 27 * a ^ 2 * d+  p = disc / (3 * a ^ 2)+  q = xx / (27 * a ^ 3)+  phi = 1 / 3 * acos (3 * q / (2 * p) * sqrt (-3 / p))+  trig k = 2 * sqrt (-p / 3) * cos (phi - k * tau / 3) - b / (3 * a)+  cubert x+    | x < 0 = -((-x) ** (1 / 3))+    | otherwise = x ** (1 / 3) {-# INLINE cubForm' #-}  -- | Solve the cubic equation ax^3 + bx^2 + cx + d = 0, returning a@@ -127,18 +115,20 @@  _cubForm_prop :: Double -> Double -> Double -> Double -> Bool _cubForm_prop a b c d = all (aboutZero' 1e-5 . eval) (cubForm a b c d)-  where eval x = a*x^3 + b*x^2 + c*x + d-           -- Basically, however large you set the tolerance it seems-           -- that quickcheck can always come up with examples where-           -- the returned solutions evaluate to something near zero-           -- but larger than the tolerance (but it takes it more-           -- tries the larger you set the tolerance). Wonder if this-           -- is an inherent limitation or (more likely) a problem-           -- with numerical stability.  If this turns out to be an-           -- issue in practice we could, say, use the solutions-           -- generated here as very good guesses to a numerical-           -- solver which can give us a more precise answer?+ where+  eval x = a * x ^ 3 + b * x ^ 2 + c * x + d +-- Basically, however large you set the tolerance it seems+-- that quickcheck can always come up with examples where+-- the returned solutions evaluate to something near zero+-- but larger than the tolerance (but it takes it more+-- tries the larger you set the tolerance). Wonder if this+-- is an inherent limitation or (more likely) a problem+-- with numerical stability.  If this turns out to be an+-- issue in practice we could, say, use the solutions+-- generated here as very good guesses to a numerical+-- solver which can give us a more precise answer?+ ------------------------------------------------------------ -- Quartic formula ------------------------------------------------------------@@ -155,33 +145,37 @@   -- x(ax^3+bx^2+cx+d)   | aboutZero' toler c0 = 0 : cubForm c4 c3 c2 c1   -- substitute solutions of y back to x-  | otherwise = map (\x->x-(a/4)) roots-    where-      -- eliminate c4: x^4+ax^3+bx^2+cx+d-      [a,b,c,d] = map (/c4) [c3,c2,c1,c0]-      -- eliminate cubic term via x = y - a/4-      -- reduced quartic: y^4 + py^2 + qy + r = 0-      p = b - 3/8*a^2-      q = 1/8*a^3-a*b/2+c-      r = (-3/256)*a^4+a^2*b/16-a*c/4+d+  | otherwise = map (\x -> x - (a / 4)) roots+ where+  -- eliminate c4: x^4+ax^3+bx^2+cx+d+  a = c3 / c4+  b = c2 / c4+  c = c1 / c4+  d = c0 / c4+  -- eliminate cubic term via x = y - a/4+  -- reduced quartic: y^4 + py^2 + qy + r = 0+  p = b - 3 / 8 * a ^ 2+  q = 1 / 8 * a ^ 3 - a * b / 2 + c+  r = (-3 / 256) * a ^ 4 + a ^ 2 * b / 16 - a * c / 4 + d -      -- | roots of the reduced quartic-      roots | aboutZero' toler r =-                0 : cubForm 1 0 p q   -- no constant term: y(y^3 + py + q) = 0-            | u < -toler || v < -toler = []     -- no real solutions due to square root-            | otherwise      = s1++s2 -- solutions of the quadratics+  -- \| roots of the reduced quartic+  roots+    | aboutZero' toler r =+        0 : cubForm 1 0 p q -- no constant term: y(y^3 + py + q) = 0+    | u < -toler || v < -toler = [] -- no real solutions due to square root+    | otherwise = s1 ++ s2 -- solutions of the quadratics -      -- solve the resolvent cubic - only one solution is needed-      z:_ = cubForm 1 (-p/2) (-r) (p*r/2 - q^2/8)+  -- solve the resolvent cubic - only one solution is needed+  z = fromMaybe 0 . listToMaybe $ cubForm 1 (-p / 2) (-r) (p * r / 2 - q ^ 2 / 8) -      -- solve the two quadratic equations-      -- y^2 ± v*y-(±u-z)-      u = z^2 - r-      v = 2*z - p-      u' = if aboutZero' toler u then 0 else sqrt u-      v' = if aboutZero' toler v then 0 else sqrt v-      s1 = quadForm 1 (if q<0 then -v' else v') (z-u')-      s2 = quadForm 1 (if q<0 then v' else -v') (z+u')+  -- solve the two quadratic equations+  -- y^2 ± v*y-(±u-z)+  u = z ^ 2 - r+  v = 2 * z - p+  u' = if aboutZero' toler u then 0 else sqrt u+  v' = if aboutZero' toler v then 0 else sqrt v+  s1 = quadForm 1 (if q < 0 then -v' else v') (z - u')+  s2 = quadForm 1 (if q < 0 then v' else -v') (z + u') {-# INLINE quartForm' #-}  -- | Solve the quartic equation c4 x^4 + c3 x^3 + c2 x^2 + c1 x + c0 = 0, returning a@@ -193,5 +187,7 @@  _quartForm_prop :: Double -> Double -> Double -> Double -> Double -> Bool _quartForm_prop a b c d e = all (aboutZero' 1e-5 . eval) (quartForm a b c d e)-  where eval x = a*x^4 + b*x^3 + c*x^2 + d*x + e-           -- Same note about tolerance as for cubic+ where+  eval x = a * x ^ 4 + b * x ^ 3 + c * x ^ 2 + d * x + e++-- Same note about tolerance as for cubic

@@ -1,5 +1,5 @@ {-# OPTIONS_GHC -fno-warn-name-shadowing #-}------------------------------------------------------------------------------+ -- | -- Module      :  Diagrams.Solve.Tridiagonal -- Copyright   :  (c) 2011-2015 diagrams-solve team (see LICENSE)@@ -7,67 +7,62 @@ -- Maintainer  :  diagrams-discuss@googlegroups.com -- -- Solving of tridiagonal and cyclic tridiagonal linear systems.----------------------------------------------------------------------------------module Diagrams.Solve.Tridiagonal-       ( solveTriDiagonal-       , solveCyclicTriDiagonal-       ) where+module Diagrams.Solve.Tridiagonal (+  solveTriDiagonal,+  solveCyclicTriDiagonal,+) where +import Data.List.NonEmpty (NonEmpty (..), (<|))+import qualified Data.List.NonEmpty as NE+ -- | @solveTriDiagonal as bs cs ds@ solves a system of the form @A*X = ds@ --   where 'A' is an 'n' by 'n' matrix with 'bs' as the main diagonal --   and 'as' the diagonal below and 'cs' the diagonal above.  See: --   <http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm>--solveTriDiagonal :: Fractional a => [a] -> [a] -> [a] -> [a] -> [a]-solveTriDiagonal as (b0:bs) (c0:cs) (d0:ds) = h cs' ds'-  where-    cs' = c0 / b0 : f cs' as bs cs-    f _ [_] _ _ = []-    f (c':cs') (a:as) (b:bs) (c:cs) = c / (b - c' * a) : f cs' as bs cs-    f _ _ _ _ = error "solveTriDiagonal.f: impossible!"--    ds' = d0 / b0 : g ds' as bs cs' ds-    g _ [] _ _ _ = []-    g (d':ds') (a:as) (b:bs) (c':cs') (d:ds) = (d - d' * a)/(b - c' * a) : g ds' as bs cs' ds-    g _ _ _ _ _ = error "solveTriDiagonal.g: impossible!"+solveTriDiagonal :: Fractional a => [a] -> NonEmpty a -> NonEmpty a -> NonEmpty a -> NonEmpty a+solveTriDiagonal as (b0 :| bs) (c0 :| cs) (d0 :| ds) = h cs' ds'+ where+  cs' = c0 / b0 : f cs' as bs cs+  f _ [_] _ _ = []+  f (c' : cs') (a : as) (b : bs) (c : cs) = c / (b - c' * a) : f cs' as bs cs+  f _ _ _ _ = error "solveTriDiagonal.f: impossible!" -    h _ [d] = [d]-    h (c:cs) (d:ds) = let xs@(x:_) = h cs ds in d - c * x : xs-    h _ _ = error "solveTriDiagonal.h: impossible!"+  ds' = d0 / b0 : g ds' as bs cs' ds+  g _ [] _ _ _ = []+  g (d' : ds') (a : as) (b : bs) (c' : cs') (d : ds) = (d - d' * a) / (b - c' * a) : g ds' as bs cs' ds+  g _ _ _ _ _ = error "solveTriDiagonal.g: impossible!" -solveTriDiagonal _ _ _ _ = error "arguments 2,3,4 to solveTriDiagonal must be nonempty"+  h _ [d] = d :| []+  h (c : cs) (d : ds) = let xs@(x :| _) = h cs ds in d - c * x <| xs+  h _ _ = error "solveTriDiagonal.h: impossible!"  -- Helper that applies the passed function only to the last element of a list modifyLast :: (a -> a) -> [a] -> [a]-modifyLast _ []     = []-modifyLast f [a]    = [f a]-modifyLast f (a:as) = a : modifyLast f as+modifyLast _ [] = []+modifyLast f [a] = [f a]+modifyLast f (a : as) = a : modifyLast f as --- Helper that builds a list of length n of the form: '[s,m,m,...,m,m,e]'-sparseVector :: Int -> a -> a -> a -> [a]-sparseVector n s m e-    | n < 1     = []-    | otherwise = s : h (n - 1)-  where-    h 1 = [e]-    h n = m : h (n - 1)+-- Helper that builds a non-empty list of the form+-- '[s,m,m,...,m,m,e]', with the same length as the given list+sparseVector :: NonEmpty x -> a -> a -> a -> NonEmpty a+sparseVector (_ :| ds) s m e = s :| h ds+ where+  h [] = []+  h [_] = [e]+  h (_ : ds) = m : h ds  -- | Solves a system similar to the tri-diagonal system using a special case --   of the Sherman-Morrison formula (<http://en.wikipedia.org/wiki/Sherman-Morrison_formula>). --   This code is based on /Numerical Recpies in C/'s @cyclic@ function in section 2.7.-solveCyclicTriDiagonal :: Fractional a => [a] -> [a] -> [a] -> [a] -> a -> a -> [a]-solveCyclicTriDiagonal as (b0:bs) cs ds alpha beta = zipWith ((+) . (fact *)) zs xs-  where-    l = length ds-    gamma = -b0-    us = sparseVector l gamma 0 alpha--    bs' = (b0 - gamma) : modifyLast (subtract (alpha*beta/gamma)) bs+solveCyclicTriDiagonal :: Fractional a => [a] -> NonEmpty a -> NonEmpty a -> NonEmpty a -> a -> a -> NonEmpty a+solveCyclicTriDiagonal as (b0 :| bs) cs ds alpha beta = NE.zipWith ((+) . (fact *)) zs xs+ where+  gamma = -b0+  us = sparseVector ds gamma 0 alpha -    xs@(x:_) = solveTriDiagonal as bs' cs ds-    zs@(z:_) = solveTriDiagonal as bs' cs us+  bs' = (b0 - gamma) :| modifyLast (subtract (alpha * beta / gamma)) bs -    fact = -(x + beta * last xs / gamma) / (1.0 + z + beta * last zs / gamma)+  xs@(x :| _) = solveTriDiagonal as bs' cs ds+  zs@(z :| _) = solveTriDiagonal as bs' cs us -solveCyclicTriDiagonal _ _ _ _ _ _ = error "second argument to solveCyclicTriDiagonal must be nonempty"+  fact = -((x + beta * NE.last xs / gamma) / (1.0 + z + beta * NE.last zs / gamma))

@@ -3,32 +3,41 @@ import Data.List (sort) import Diagrams.Solve.Polynomial -import Test.Tasty (defaultMain, testGroup, TestTree)-import           Test.Tasty.QuickCheck+import Test.Tasty (TestTree, defaultMain, testGroup)+import Test.Tasty.QuickCheck  tests :: TestTree-tests = testGroup "Solve" [-         testProperty "solutions found satisfy quadratic equation" $-         \a b c -> let sat x =  a * x * x + b * x + c =~ 0 in all sat (quadForm a b c)--- could verify number of solutions, but we would just duplicate the function definition-        , testProperty "solutions found satisfy cubic equation" $-         \a b c d -> let sat x =  a * x * x * x + b * x * x + c * x + d =~ (0 :: Double) in all sat (cubForm a b c d)---- some specific examples and regression tests-        , testGroup "Solve specific examples" [-            testProperty "1 * x^3 + -886.7970773009183 * x^2 + 262148.4783430062 * x + -264000817.775054 = 0" $-                let [r] = cubForm 1 (-886.7970773009183) 262148.4783430062 (-264000817.775054) in-                r =~ 915.4538593912--          , testProperty "1 * u^4 + -240 * u^3 + 25449 * u^2 + -1325880 * u + 26471900.25 = 0" $-                let [r1, r2] = sort $ quartForm 1 (-240) 25449 (-1325880) 26471900.25 in-                r1 =~ 50.6451 && r2 =~ 69.3549-            ]+tests =+  testGroup+    "Solve"+    [ testProperty "solutions found satisfy quadratic equation" $+        \a b c -> let sat x = a * x * x + b * x + c =~ (0 :: Double) in all sat (quadForm a b c)+    , -- could verify number of solutions, but we would just duplicate the function definition+      testProperty "solutions found satisfy cubic equation" $+        \a b c d -> let sat x = a * x * x * x + b * x * x + c * x + d =~ (0 :: Double) in all sat (cubForm a b c d)+    , -- some specific examples and regression tests+      testGroup+        "Solve specific examples"+        [ testProperty "1 * x^3 + -886.7970773009183 * x^2 + 262148.4783430062 * x + -264000817.775054 = 0" $+            let rs = cubForm 1 (-886.7970773009183) 262148.4783430062 (-264000817.775054)+             in rs =~ [915.4538593912 :: Double]+        , testProperty "1 * u^4 + -240 * u^3 + 25449 * u^2 + -1325880 * u + 26471900.25 = 0" $+            let rs = sort $ quartForm 1 (-240) 25449 (-1325880) 26471900.25+             in rs =~ [50.6451, 69.3549 :: Double]         ]+    ] -(=~) :: Double -> Double -> Bool-(=~) a b = abs (a - b) < 0.001+class Eqish a where+  (=~) :: a -> a -> Bool infix 4 =~++instance Eqish Double where+  (=~) a b = abs (a - b) < 0.001++instance Eqish a => Eqish [a] where+  [] =~ [] = True+  (a : as) =~ (b : bs) = a =~ b && as =~ bs+  _ =~ _ = False  main :: IO () main = defaultMain tests