diff --git a/CHANGES.markdown b/CHANGES.markdown --- a/CHANGES.markdown +++ b/CHANGES.markdown @@ -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 diff --git a/diagrams-solve.cabal b/diagrams-solve.cabal --- a/diagrams-solve.cabal +++ b/diagrams-solve.cabal @@ -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 diff --git a/src/Diagrams/Solve/Polynomial.hs b/src/Diagrams/Solve/Polynomial.hs --- a/src/Diagrams/Solve/Polynomial.hs +++ b/src/Diagrams/Solve/Polynomial.hs @@ -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 diff --git a/src/Diagrams/Solve/Tridiagonal.hs b/src/Diagrams/Solve/Tridiagonal.hs --- a/src/Diagrams/Solve/Tridiagonal.hs +++ b/src/Diagrams/Solve/Tridiagonal.hs @@ -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: -- - -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 (). -- 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)) diff --git a/tests/Test.hs b/tests/Test.hs --- a/tests/Test.hs +++ b/tests/Test.hs @@ -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