diff --git a/CHANGES.markdown b/CHANGES.markdown new file mode 100644 --- /dev/null +++ b/CHANGES.markdown @@ -0,0 +1,4 @@ +* 0.1 (19 April 2015) + + initial release, in conjunction with `diagrams-1.3` --- some + functionality split out from `diagrams-lib` diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.markdown b/README.markdown new file mode 100644 --- /dev/null +++ b/README.markdown @@ -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` diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/diagrams-solve.cabal b/diagrams-solve.cabal new file mode 100644 --- /dev/null +++ b/diagrams-solve.cabal @@ -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 diff --git a/src/Diagrams/Solve/Polynomial.hs b/src/Diagrams/Solve/Polynomial.hs new file mode 100644 --- /dev/null +++ b/src/Diagrams/Solve/Polynomial.hs @@ -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 diff --git a/src/Diagrams/Solve/Tridiagonal.hs b/src/Diagrams/Solve/Tridiagonal.hs new file mode 100644 --- /dev/null +++ b/src/Diagrams/Solve/Tridiagonal.hs @@ -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: +-- + +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 (). +-- 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"