diagrams-solve (empty) → 0.1

raw patch · 7 files changed

+345/−0 lines, 7 filesdep +basesetup-changed

Dependencies added: base

Files

CHANGES.markdown +4/−0
LICENSE +30/−0
README.markdown +18/−0
Setup.hs +2/−0
diagrams-solve.cabal +28/−0
src/Diagrams/Solve/Polynomial.hs +190/−0
src/Diagrams/Solve/Tridiagonal.hs +73/−0

+ CHANGES.markdown view

@@ -0,0 +1,4 @@+* 0.1 (19 April 2015)++  initial release, in conjunction with `diagrams-1.3` --- some+  functionality split out from `diagrams-lib`

+ LICENSE view

@@ -0,0 +1,30 @@+Copyright (c) 2015, various++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of various nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

+ README.markdown view

@@ -0,0 +1,18 @@+[![Build Status](https://travis-ci.org/diagrams/diagrams-solve.png?branch=master)](https://travis-ci.org/diagrams/diagrams-solve)++Miscellaneous pure-Haskell solver routines used in+[diagrams](http://projects.haskell.org/diagrams/), a Haskell embedded+domain-specific language for compositional, declarative drawing.++This is split out into a separate package with no dependencies on the+rest of diagrams in case it is useful to others, but no particular+guarantees are made as to the suitability or correctness of the code+(though we are certainly open to bug reports).++Currently the package contains:++  - functions to find real roots of quadratic, cubic, and quartic+    polynomials, in `Diagrams.Solve.Polynomial`++  - functions to solve tridiagonal and cyclic tridiagonal systems of+    linear equations, in `Diagrams.Solve.Tridiagonal`

+ Setup.hs view

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

+ diagrams-solve.cabal view

@@ -0,0 +1,28 @@+name:                diagrams-solve+version:             0.1+synopsis:            Pure Haskell solver routines used by diagrams+description:         Pure Haskell solver routines used by the diagrams+                     project.  Currently includes finding real roots+                     of low-degree (n < 5) polynomials, and solving+                     tridiagonal and cyclic tridiagonal linear+                     systems.+homepage:            http://projects.haskell.org/diagrams+license:             BSD3+license-file:        LICENSE+author:              various+maintainer:          diagrams-discuss@googlegroups.com+category:            Math+build-type:          Simple+extra-source-files:  README.markdown, CHANGES.markdown+cabal-version:       >=1.10+Tested-with:         GHC == 7.4.2, GHC == 7.6.3, GHC == 7.8.4, GHC == 7.10.1+Source-repository head+  type:     git+  location: http://github.com/diagrams/diagrams-solve.git++library+  exposed-modules:     Diagrams.Solve.Polynomial,+                       Diagrams.Solve.Tridiagonal+  build-depends:       base >=4.5 && < 4.9+  hs-source-dirs:      src+  default-language:    Haskell2010

+ src/Diagrams/Solve/Polynomial.hs view

@@ -0,0 +1,190 @@+-----------------------------------------------------------------------------+-- |+-- Module      :  Diagrams.Solve.Polynomial+-- Copyright   :  (c) 2011-2015 diagrams-solve team (see LICENSE)+-- License     :  BSD-style (see LICENSE)+-- 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)++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++-- | A specialization of (^) to Integer+--   c.f. http://comments.gmane.org/gmane.comp.lang.haskell.libraries/21164+--   for discussion. "The choice in (^) and (^^) to overload on the+--   power's Integral type... was a genuinely bad idea." - Edward Kmett+(^) :: (Num a) => a -> Integer -> a+(^) = (P.^)++-- | Utility function used to avoid singularities+aboutZero' :: (Ord a, Num a) => a -> a -> Bool+aboutZero' toler x = abs x < toler++------------------------------------------------------------+-- Quadratic formula+------------------------------------------------------------++-- | 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+  | a == 0 && b == 0 && c == 0 = [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)++_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++------------------------------------------------------------+-- Cubic formula+------------------------------------------------------------++-- See http://en.wikipedia.org/wiki/Cubic_formula#General_formula_of_roots++-- | Solve the cubic equation ax^3 + bx^2 + cx + d = 0, returning a+--   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'    | aboutZero' toler disc = maximumBy (comparing (abs . (+xx))) [qq, -qq]+              | otherwise = 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)++-- | Solve the cubic equation ax^3 + bx^2 + cx + d = 0, returning a+--   list of all real roots within 1e-10 tolerance+--   (although currently it's closer to 1e-5)+cubForm :: (Floating d, Ord d) => d -> d -> d -> d -> [d]+cubForm = cubForm' 1e-10++_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?++------------------------------------------------------------+-- Quartic formula+------------------------------------------------------------++-- Based on http://tog.acm.org/resources/GraphicsGems/gems/Roots3b/and4.c+-- as of 5/12/14, with help from http://en.wikipedia.org/wiki/Quartic_function++-- | Solve the quartic equation c4 x^4 + c3 x^3 + c2 x^2 + c1 x + c0 = 0, returning a+--   list of all real roots. First argument is tolerance.+quartForm' :: (Floating d, Ord d) => d -> d -> d -> d -> d -> d -> [d]+quartForm' toler c4 c3 c2 c1 c0+  -- obvious cubic+  | aboutZero' toler c4 = cubForm c3 c2 c1 c0+  -- 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++      -- | 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 < 0 || v < 0 = []     -- 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 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 quartic equation c4 x^4 + c3 x^3 + c2 x^2 + c1 x + c0 = 0, returning a+--   list of all real roots within 1e-10 tolerance+--   (although currently it's closer to 1e-5)+quartForm :: (Floating d, Ord d) => d -> d -> d -> d -> d -> [d]+quartForm = quartForm' 1e-10++_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

+ src/Diagrams/Solve/Tridiagonal.hs view

@@ -0,0 +1,73 @@+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Diagrams.Solve.Tridiagonal+-- Copyright   :  (c) 2011-2015 diagrams-solve team (see LICENSE)+-- License     :  BSD-style (see LICENSE)+-- Maintainer  :  diagrams-discuss@googlegroups.com+--+-- Solving of tridiagonal and cyclic tridiagonal linear systems.+--+-----------------------------------------------------------------------------+module Diagrams.Solve.Tridiagonal+       ( solveTriDiagonal+       , solveCyclicTriDiagonal+       ) where++-- | @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!"++    h _ [d] = [d]+    h (c:cs) (d:ds) = let xs@(x:_) = h cs ds in d - c * x : xs+    h _ _ = error "solveTriDiagonal.h: impossible!"++solveTriDiagonal _ _ _ _ = error "arguments 2,3,4 to solveTriDiagonal must be nonempty"++-- 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++-- 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)++-- | 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++    xs@(x:_) = solveTriDiagonal as bs' cs ds+    zs@(z:_) = solveTriDiagonal as bs' cs us++    fact = -(x + beta * last xs / gamma) / (1.0 + z + beta * last zs / gamma)++solveCyclicTriDiagonal _ _ _ _ _ _ = error "second argument to solveCyclicTriDiagonal must be nonempty"