diff --git a/JacobiRoots.hs b/JacobiRoots.hs deleted file mode 100644 # file too large to diff: JacobiRoots.hs diff --git a/Makefile b/Makefile new file mode 100644 --- /dev/null +++ b/Makefile @@ -0,0 +1,7 @@ +all : src/JacobiRootsBinary.hs + +src/JacobiRootsBinary.hs : gen/JacobiRootsRaw300.hs gen/Convert.hs + cd gen && runhaskell -Wall -Werror Convert.hs + +clean : + rm -f src/JacobiRootsBinary.hs diff --git a/compute_roots.nb b/compute_roots.nb --- a/compute_roots.nb +++ b/compute_roots.nb @@ -10,10 +10,10 @@ NotebookFileLineBreakTest NotebookFileLineBreakTest NotebookDataPosition[ 157, 7] -NotebookDataLength[ 12371, 294] -NotebookOptionsPosition[ 12000, 276] -NotebookOutlinePosition[ 12357, 292] -CellTagsIndexPosition[ 12314, 289] +NotebookDataLength[ 14346, 313] +NotebookOptionsPosition[ 13970, 295] +NotebookOutlinePosition[ 14328, 311] +CellTagsIndexPosition[ 14285, 308] WindowFrame->Normal*) (* Beginning of Notebook Content *) @@ -80,7 +80,7 @@ RowBox[{ RowBox[{ RowBox[{"funName", "[", "name_", "]"}], ":=", - RowBox[{"\"\\"", "<>", "name", "<>", "\"\\""}]}], + RowBox[{"\"\\"", "<>", "name", "<>", "\"\\""}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ @@ -91,26 +91,11 @@ RowBox[{"ToString", "[", "roots", "]"}], ",", RowBox[{"{", RowBox[{ - RowBox[{"\"\<{{\>\"", "\[Rule]", "\"\< [ [\>\""}], ",", + RowBox[{"\"\<{{\>\"", "\[Rule]", "\"\< [ []\\n , [\>\""}], ",", RowBox[{"\"\<{\>\"", "\[Rule]", "\"\<[\>\""}], ",", RowBox[{"\"\<},\>\"", "\[Rule]", "\"\<]\\n ,\>\""}], ",", RowBox[{"\"\<}}\>\"", "\[Rule]", "\"\<]\\n ]\\n\>\""}]}], "}"}]}], - "]"}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", - RowBox[{"(*", - RowBox[{ - RowBox[{ - RowBox[{"specialize", "[", - RowBox[{"name_", ",", "type_"}], "]"}], ":=", - RowBox[{"\"\<{-# SPECIALIZE \>\"", "<>", - RowBox[{"funName", "[", "name", "]"}], "<>", "\"\< :: [[\>\"", "<>", - "type", "<>", "\"\<]] #-}\\n\>\""}]}], ";"}], - "*)"}]}], "\[IndentingNewLine]", - RowBox[{ - RowBox[{ - RowBox[{ - RowBox[{"specialize", "[", - RowBox[{"name_", ",", "type_"}], "]"}], ":=", "\"\<\>\""}], ";"}], - "\[IndentingNewLine]"}], "\[IndentingNewLine]", + "]"}]}], ";"}], "\[IndentingNewLine]"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{ @@ -119,7 +104,7 @@ RowBox[{"\"\<-- | roots of shifted Jacobi polynomials with alpha=\>\"", "<>", RowBox[{"ToString", "[", "\[Alpha]", "]"}], "<>", "\"\<, beta=\>\"", "<>", RowBox[{"ToString", "[", "\[Beta]", "]"}], "<>", "\"\<\\n\>\""}]}], - ";"}], "\[IndentingNewLine]", "\n", + ";"}], "\n", "\[IndentingNewLine]", RowBox[{"(*", RowBox[{ RowBox[{"typeSignature", "[", "name_", "]"}], ":=", @@ -142,10 +127,6 @@ RowBox[{ RowBox[{"doc", "[", RowBox[{"name", ",", "\[Alpha]", ",", "\[Beta]"}], "]"}], "<>", - RowBox[{"specialize", "[", - RowBox[{"name", ",", "\"\\""}], "]"}], "<>", - RowBox[{"specialize", "[", - RowBox[{"name", ",", "\"\\""}], "]"}], "<>", RowBox[{"typeSignature", "[", "name", "]"}], "<>", RowBox[{"funName", "[", "name", "]"}], "<>", "\"\< =\\n\>\"", "<>", RowBox[{"toHaskell", "[", @@ -155,10 +136,13 @@ RowBox[{ RowBox[{ RowBox[{"header", " ", "=", " ", - RowBox[{"\"\<{-# OPTIONS_GHC -Wall #-}\\n\\nmodule JacobiRoots ( \>\"", "<>", + RowBox[{ + "\"\<\\\n{-# OPTIONS_GHC -Wall #-}\n\n-- | Module with raw data, \ +generated from Mathematica\\n\\nmodule JacobiRootsRaw\>\"", "<>", + RowBox[{"ToString", "[", "order", "]"}], "<>", "\"\< ( \>\"", "<>", RowBox[{"funName", "[", "\"\\"", "]"}], "<>", "\"\<, \>\"", "<>", RowBox[{"funName", "[", "\"\\"", "]"}], "<>", - "\"\< ) where\\n\>\""}]}], ";"}], + "\"\< ) where\n\>\""}]}], ";"}], "\[IndentingNewLine]"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"t0", "=", @@ -203,8 +187,9 @@ RowBox[{ RowBox[{ RowBox[{"NotebookDirectory", "[", "]"}], "<>", - "\"\\""}], ",", "module", ",", "\"\\""}], - "]"}], ";"}]}], "Input", + "\"\\"", "<>", + RowBox[{"ToString", "[", "order", "]"}], "<>", "\"\<.hs\>\""}], ",", + "module", ",", "\"\\""}], "]"}], ";"}]}], "Input", CellChangeTimes->{{3.587987975710743*^9, 3.587988018967742*^9}, { 3.587988140782954*^9, 3.58798815728058*^9}, {3.587988288124055*^9, 3.587988349210535*^9}, {3.587988379483639*^9, 3.587988428435318*^9}, { @@ -250,9 +235,23 @@ 3.5880149096937437`*^9}, {3.588015249321007*^9, 3.588015250099921*^9}, { 3.588015742672408*^9, 3.588015751994816*^9}, {3.58801578267706*^9, 3.588015964903317*^9}, {3.5880160012766333`*^9, 3.588016029194995*^9}, { - 3.5880163126420937`*^9, 3.588016312947206*^9}}], + 3.5880163126420937`*^9, 3.588016312947206*^9}, {3.5880492477726107`*^9, + 3.588049265665798*^9}, 3.5880498267091303`*^9, {3.588050037648439*^9, + 3.588050063191792*^9}, 3.5880501433990107`*^9, {3.5880502595720873`*^9, + 3.5880502620409193`*^9}, {3.588050305322728*^9, 3.588050308864871*^9}, { + 3.58805340225362*^9, 3.5880534078772917`*^9}, {3.588053522769102*^9, + 3.5880535309152517`*^9}, {3.5880537119087257`*^9, 3.588053715411994*^9}, + 3.588057517471484*^9, {3.588057587804565*^9, 3.588057627662592*^9}, { + 3.588057658967712*^9, 3.588057703410576*^9}, {3.5880577363170643`*^9, + 3.5880577586845617`*^9}, 3.5880585557306213`*^9, {3.588064533314286*^9, + 3.58806455128642*^9}, {3.588064695121482*^9, 3.5880646961299763`*^9}, { + 3.588064790677484*^9, 3.588064790978847*^9}, {3.588064939633843*^9, + 3.588064948600685*^9}, {3.588065420861298*^9, 3.58806542302002*^9}, { + 3.5880669915461607`*^9, 3.58806700833703*^9}, {3.5880670393550777`*^9, + 3.588067048968004*^9}, {3.5880671241013193`*^9, 3.588067131159873*^9}, { + 3.5880672290497227`*^9, 3.588067229613203*^9}}], -Cell[BoxData["\<\"Radau roots computed in 99.155180 seconds\"\>"], "Output", +Cell[BoxData["\<\"Radau roots computed in 102.655302 seconds\"\>"], "Output", CellChangeTimes->{ 3.5880141771904573`*^9, 3.588014379606248*^9, 3.588014422221437*^9, { 3.588014530934038*^9, 3.588014551780951*^9}, 3.588014651040941*^9, { @@ -260,9 +259,19 @@ 3.588014818472938*^9, 3.5880148466863537`*^9}, 3.588014890897044*^9, 3.588014959037346*^9, 3.588015596428967*^9, 3.588015916598996*^9, { 3.5880159517687273`*^9, 3.588015965603405*^9}, 3.588016008640575*^9, - 3.588016058468951*^9, 3.5880164175632877`*^9, 3.588017582282074*^9}], + 3.588016058468951*^9, 3.5880164175632877`*^9, 3.588017582282074*^9, + 3.588050038856845*^9, 3.588050073854802*^9, 3.588050145238899*^9, + 3.588050263231791*^9, 3.588050310958125*^9, 3.5880529748624268`*^9, { + 3.588053402797201*^9, 3.588053408258072*^9}, {3.588053525294703*^9, + 3.588053531424819*^9}, 3.588053715749647*^9, 3.5880554757821198`*^9, + 3.588057597999413*^9, 3.588057629157517*^9, {3.588057685296475*^9, + 3.5880577076867027`*^9}, 3.588057760554085*^9, 3.5880578072788057`*^9, + 3.5880589719227962`*^9, 3.588064553163354*^9, 3.588064701928281*^9, + 3.5880647923135643`*^9, {3.588064940167183*^9, 3.588064949406036*^9}, + 3.58806583592706*^9, {3.5880670401113577`*^9, 3.588067051889412*^9}, + 3.588067133405834*^9, 3.588067335364573*^9}], -Cell[BoxData["\<\"Legendre roots computed in 362.251589 seconds\"\>"], \ +Cell[BoxData["\<\"Legendre roots computed in 376.058679 seconds\"\>"], \ "Output", CellChangeTimes->{ 3.5880141771904573`*^9, 3.588014379606248*^9, 3.588014422221437*^9, { @@ -271,10 +280,20 @@ 3.588014818472938*^9, 3.5880148466863537`*^9}, 3.588014890897044*^9, 3.588014959037346*^9, 3.588015596428967*^9, 3.588015916598996*^9, { 3.5880159517687273`*^9, 3.588015965603405*^9}, 3.588016008640575*^9, - 3.588016058468951*^9, 3.5880164175632877`*^9, 3.58801794453237*^9}] + 3.588016058468951*^9, 3.5880164175632877`*^9, 3.588017582282074*^9, + 3.588050038856845*^9, 3.588050073854802*^9, 3.588050145238899*^9, + 3.588050263231791*^9, 3.588050310958125*^9, 3.5880529748624268`*^9, { + 3.588053402797201*^9, 3.588053408258072*^9}, {3.588053525294703*^9, + 3.588053531424819*^9}, 3.588053715749647*^9, 3.5880554757821198`*^9, + 3.588057597999413*^9, 3.588057629157517*^9, {3.588057685296475*^9, + 3.5880577076867027`*^9}, 3.588057760554085*^9, 3.5880578072788057`*^9, + 3.5880589719227962`*^9, 3.588064553163354*^9, 3.588064701928281*^9, + 3.5880647923135643`*^9, {3.588064940167183*^9, 3.588064949406036*^9}, + 3.58806583592706*^9, {3.5880670401113577`*^9, 3.588067051889412*^9}, + 3.588067133405834*^9, 3.588067711421949*^9}] }, Open ]] }, -WindowSize->{807, 874}, +WindowSize->{1031, 975}, WindowMargins->{{11, Automatic}, {Automatic, 10}}, FrontEndVersion->"8.0 for Mac OS X x86 (32-bit, 64-bit Kernel) (February 23, \ 2011)", @@ -292,9 +311,9 @@ (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ -Cell[579, 22, 10188, 230, 658, "Input"], -Cell[10770, 254, 603, 8, 27, "Output"], -Cell[11376, 264, 608, 9, 27, "Output"] +Cell[579, 22, 10754, 229, 553, "Input"], +Cell[11336, 253, 1305, 18, 27, "Output"], +Cell[12644, 273, 1310, 19, 27, "Output"] }, Open ]] } ] diff --git a/gen/Convert.hs b/gen/Convert.hs new file mode 100644 --- /dev/null +++ b/gen/Convert.hs @@ -0,0 +1,27 @@ +{-# OPTIONS_GHC -Wall #-} + +module Main ( main ) where + +import Data.Binary ( encode ) + +import JacobiRootsRaw300 + +main :: IO () +main = do + writeFile "../src/JacobiRootsBinary.Hs" blah + +blah :: String +blah = unlines + [ "{-# OPTIONS_GHC -Wall #-}" + , "{-# Language OverloadedStrings #-}" + , "" + , "module JacobiRootsBinary ( allShiftedLegendreRootsBinary, allShiftedRadauRootsBinary ) where" + , "" + , "import Data.ByteString.Lazy ( ByteString )" + , "" + , "allShiftedLegendreRootsBinary :: ByteString" + , "allShiftedLegendreRootsBinary = " ++ show (encode allShiftedLegendreRootsRaw) + , "" + , "allShiftedRadauRootsBinary :: ByteString" + , "allShiftedRadauRootsBinary = " ++ show (encode allShiftedRadauRootsRaw) + ] diff --git a/jacobi-roots.cabal b/jacobi-roots.cabal --- a/jacobi-roots.cabal +++ b/jacobi-roots.cabal @@ -1,19 +1,50 @@ name: jacobi-roots -version: 0.1.0.0 -synopsis: roots of two shifted jacobi polynomials to double precision -description: This package provides pre-computed roots to the shifted Radau and Legendre polynomials. Roots to order of a few hundred are provided. +version: 0.2.0.1 license: GPL-3 license-file: LICENSE author: Greg Horn maintainer: gregmainland@gmail.com copyright: (c) 2013 Greg Horn category: Math +homepage: http://github.com/ghorn/jacobi-roots +bug-reports: http://github.com/ghorn/jacobi-roots/issues build-type: Simple -extra-source-files: compute_roots.nb cabal-version: >=1.10 +synopsis: Roots of two shifted Jacobi polynomials (Legendre and Radau) to double precision +description: + This package provides roots of shifted Radau and Legendre polynomials to double precision. + . + The roots are precomputed using a Mathematica script, and degrees of up to a few hundred are provided. + . + See the "JacobiRoots" module for documentation + . + To build this from scratch, first run compute_roots.nb in Mathematica to generate gen/JacobiRootsRaw300.hs. Then type "make" to run gen/Convert.hs, which creates src/JacobiRootsBinary.hs. Then the normal "cabal configure; cabal build" +extra-source-files: compute_roots.nb + gen/Convert.hs +-- gen/JacobiRootsRaw300.hs + Makefile + library + hs-source-dirs: src exposed-modules: JacobiRoots - build-depends: base >=4.6 && <4.7 + other-modules: JacobiRootsBinary + build-depends: base >= 4 && < 5, + vector, + bytestring, + binary default-language: Haskell2010 ghc-options: -O2 + +source-repository head + type: git + location: git://github.com/ghorn/jacobi-roots.git + +test-suite doctests + type: exitcode-stdio-1.0 + main-is: doctests.hs + build-depends: base >= 4 && < 5, + doctest >= 0.8 + default-language: Haskell2010 + ghc-options: -threaded + hs-source-dirs: tests diff --git a/src/JacobiRoots.hs b/src/JacobiRoots.hs new file mode 100644 --- /dev/null +++ b/src/JacobiRoots.hs @@ -0,0 +1,75 @@ +{-# OPTIONS_GHC -Wall #-} + +-- | Get the roots of shifted Legendre and Radau polynomials +-- +-- >>> shiftedLegendreRoots 3 +-- Just (fromList [0.11270166537925831,0.5,0.8872983346207417]) +-- >>> shiftedRadauRoots 2 +-- Just (fromList [0.1550510257216822,0.6449489742783178]) +-- +-- The roots are pre-computed and only a finite number are provided +-- +-- >>> (V.length allShiftedLegendreRoots, V.length allShiftedRadauRoots) +-- (301,301) +-- >>> shiftedLegendreRoots 1000000000000000 +-- Nothing +-- +-- There are N roots in the Nth Jacobi polynomial +-- +-- >>> import Data.Maybe ( fromJust ) +-- >>> V.length (fromJust (shiftedLegendreRoots 5)) +-- 5 +-- >>> all (\k -> k == V.length (fromJust (shiftedLegendreRoots k))) [0..(V.length allShiftedLegendreRoots - 1)] +-- True +-- >>> all (\k -> k == V.length (fromJust (shiftedRadauRoots k))) [0..(V.length allShiftedRadauRoots - 1)] +-- True + +module JacobiRoots + ( shiftedLegendreRoots + , shiftedRadauRoots + , allShiftedLegendreRoots + , allShiftedRadauRoots + ) where + +import Data.Binary ( decode ) +import qualified Data.Vector as V +import Data.Vector ( (!?) ) + +import JacobiRootsBinary ( allShiftedLegendreRootsBinary, allShiftedRadauRootsBinary ) + +-- | get the roots of the Nth shifted Legendre polynomial +-- +-- @ +-- 'shiftedLegendreRoots' == ('allShiftedLegendreRoots' '!?') +-- @ +-- +-- >>> mapM_ (print . shiftedLegendreRoots) [0..3] +-- Just (fromList []) +-- Just (fromList [0.5]) +-- Just (fromList [0.2113248654051871,0.7886751345948129]) +-- Just (fromList [0.11270166537925831,0.5,0.8872983346207417]) +shiftedLegendreRoots :: Int -> Maybe (V.Vector Double) +shiftedLegendreRoots = (allShiftedLegendreRoots !?) + +-- | get the roots of the Nth shifted Radau polynomial +-- +-- @ +-- 'shiftedRadauRoots' == ('allShiftedRadauRoots' '!?') +-- @ +-- +-- >>> mapM_ (print . shiftedRadauRoots) [0..3] +-- Just (fromList []) +-- Just (fromList [0.3333333333333333]) +-- Just (fromList [0.1550510257216822,0.6449489742783178]) +-- Just (fromList [8.858795951270394e-2,0.4094668644407347,0.787659461760847]) +shiftedRadauRoots :: Int -> Maybe (V.Vector Double) +shiftedRadauRoots = (allShiftedRadauRoots !?) + + +-- | roots of shifted Jacobi polynomials with alpha=0, beta=0 +allShiftedLegendreRoots :: V.Vector (V.Vector Double) +allShiftedLegendreRoots = V.fromList $ map V.fromList $ decode allShiftedLegendreRootsBinary + +-- | roots of shifted Jacobi polynomials with alpha=1, beta=0 +allShiftedRadauRoots :: V.Vector (V.Vector Double) +allShiftedRadauRoots = V.fromList $ map V.fromList $ decode allShiftedRadauRootsBinary diff --git a/src/JacobiRootsBinary.hs b/src/JacobiRootsBinary.hs new file mode 100644 # file too large to diff: src/JacobiRootsBinary.hs diff --git a/tests/doctests.hs b/tests/doctests.hs new file mode 100644 --- /dev/null +++ b/tests/doctests.hs @@ -0,0 +1,8 @@ +{-# OPTIONS_GHC -Wall #-} + +module Main where + +import Test.DocTest + +main :: IO () +main = doctest ["-isrc","src/JacobiRoots.hs"]