diff --git a/cyclotomic.cabal b/cyclotomic.cabal --- a/cyclotomic.cabal +++ b/cyclotomic.cabal @@ -1,5 +1,6 @@ Name: cyclotomic -Version: 0.4.4.1 +Version: 0.5.0.0 +Stability: experimental Synopsis: A subfield of the complex numbers for exact calculation. Description: The cyclotomic numbers are a subset of the complex numbers that are represented exactly, enabling exact @@ -16,16 +17,33 @@ License-file: LICENSE Author: Scott N. Walck Maintainer: Scott N. Walck +Copyright: (c) Scott N. Walck 2012-2017 Category: Math Build-type: Simple -Cabal-version: >=1.6 -Tested-with: GHC == 7.8.4, GHC == 7.10.2 +Extra-source-files: test/Properties.hs +Cabal-version: >= 1.10 +Tested-with: GHC == 7.8.4, GHC == 7.10.2, GHC == 7.10.3, GHC == 8.0.1 +Bug-reports: https://github.com/walck/cyclotomic/issues + +Test-suite cyclotomic-tests + type: exitcode-stdio-1.0 + main-is: Properties.hs + build-depends: base, QuickCheck >= 2.4, cyclotomic, + test-framework, HUnit, test-framework-hunit, + test-framework-quickcheck2, + test-framework-smallcheck + default-language: Haskell2010 + Hs-source-dirs: test + Library - Exposed-modules: Data.Complex.Cyclotomic - Build-depends: base >= 4.2 && < 5, - containers >= 0.3, - arithmoi >= 0.4 - Hs-source-dirs: src + ghc-options: -Wall + Exposed-modules: Data.Complex.Cyclotomic, Data.Number.RealCyclotomic + Build-depends: base >= 4.2 && < 4.11, + containers >= 0.3, + arithmoi >= 0.4 + default-language: Haskell2010 + Hs-source-dirs: src + Source-repository head type: git - location: https://github.com/walck/cyclotomic + location: https://github.com/walck/cyclotomic.git diff --git a/src/Data/Complex/Cyclotomic.hs b/src/Data/Complex/Cyclotomic.hs --- a/src/Data/Complex/Cyclotomic.hs +++ b/src/Data/Complex/Cyclotomic.hs @@ -3,7 +3,7 @@ {- | Module : Data.Complex.Cyclotomic -Copyright : (c) Scott N. Walck 2012-2013 +Copyright : (c) Scott N. Walck 2012-2017 License : GPL-3 (see LICENSE) Maintainer : Scott N. Walck Stability : experimental @@ -86,6 +86,7 @@ , dft , dftInv , rootsQuadEq + , heron ) where @@ -473,13 +474,13 @@ isGaussianRat c = isRat (real c) && isRat (imag c) -- | Export as an inexact complex number. -toComplex :: Cyclotomic -> Complex Double +toComplex :: RealFloat a => Cyclotomic -> Complex a toComplex c = sum [fromRational r * en^p | (p,r) <- M.toList (coeffs c)] where en = exp (0 :+ 2*pi/n) n = fromIntegral (order c) -- | Export as an inexact real number if possible. -toReal :: Cyclotomic -> Maybe Double +toReal :: RealFloat a => Cyclotomic -> Maybe a toReal c | isReal c = Just $ realPart (toComplex c) | otherwise = Nothing @@ -561,3 +562,14 @@ aa = fromRational a bb = fromRational b sqrtDisc = sqrtRat (b*b - 4*a*c) + +-- | Heron's formula for the area of a triangle with +-- side lengths a, b, c. +heron :: Rational -- ^ a + -> Rational -- ^ b + -> Rational -- ^ c + -> Cyclotomic -- ^ area of triangle +heron a b c + = sqrtRat (s * (s-a) * (s-b) * (s-c)) + where + s = (a + b + c) / 2 diff --git a/src/Data/Number/RealCyclotomic.hs b/src/Data/Number/RealCyclotomic.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/RealCyclotomic.hs @@ -0,0 +1,149 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Safe #-} + +{- | +Module : Data.Number.RealCyclotomic +Copyright : (c) Scott N. Walck 2012-2017 +License : GPL-3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +The real cyclotomic numbers are a subset of the real numbers with +the following properties: + + 1. The real cyclotomic numbers are represented exactly, enabling exact + computations and equality comparisons. + + 2. The real cyclotomic numbers contain the rationals. + As a consequence, the real cyclotomic numbers are a dense subset of the + real numbers. + + 3. The real cyclotomic numbers contain the square roots of all nonnegative rational numbers. + + 4. The real cyclotomic numbers form a field: they are closed under addition, subtraction, + multiplication, and division. + + 5. The real cyclotomic numbers contain the sine and cosine of all rational + multiples of pi (equivalently, the sine and cosine of any rational number + of degrees or any rational number of revolutions). + + Floating point numbers do not do well with equality comparison: + +>(sqrt 2 + sqrt 3)^2 == 5 + 2 * sqrt 6 +> -> False + + "Data.Number.RealCyclotomic" represents these numbers exactly, allowing equality comparison: + +>(sqrtRat 2 + sqrtRat 3)^2 == 5 + 2 * sqrtRat 6 +> -> True + + 'RealCyclotomic's can be exported as inexact real numbers using the 'toReal' function: + +>sqrtRat 2 +> -> e(8) - e(8)^3 +>toReal $ sqrtRat 2 +> -> 1.414213562373095 + +This module is based on the module 'Data.Complex.Cyclotomic'. +Usually you would only import one of the modules 'Data.Number.RealCyclotomic' +or 'Data.Complex.Cyclotomic', depending on whether you wanted only +real numbers (this module) or complex numbers (the other). +Functions such as @sqrtRat@, @sinDeg@, @cosDeg@ are defined +in both modules, with different type signatures, so their +names will conflict if both modules are imported. +-} + +module Data.Number.RealCyclotomic + ( RealCyclotomic + , sqrtRat + , sinDeg + , cosDeg + , sinRev + , cosRev + , isRat + , toRat + , toReal + , goldenRatio + , heron + ) + where + +import qualified Data.Complex.Cyclotomic as Cyc +import Data.Complex + ( realPart + ) + +-- | A real cyclotomic number. +newtype RealCyclotomic = RealCyclotomic Cyc.Cyclotomic + deriving (Eq) + +-- | @abs@ and @signum@ are undefined. +instance Num RealCyclotomic where + RealCyclotomic x + RealCyclotomic y = RealCyclotomic (x + y) + RealCyclotomic x - RealCyclotomic y = RealCyclotomic (x - y) + RealCyclotomic x * RealCyclotomic y = RealCyclotomic (x * y) + negate (RealCyclotomic x) = RealCyclotomic (negate x) + fromInteger n = RealCyclotomic (fromInteger n) + abs = undefined + signum = undefined + +instance Fractional RealCyclotomic where + recip (RealCyclotomic x) = RealCyclotomic (recip x) + fromRational r = RealCyclotomic (fromRational r) + +instance Show RealCyclotomic where + show (RealCyclotomic x) = show x + +-- I need to do Ord first. +-- A Real instance would make realToFrac work. +-- instance Real RealCyclotomic where +-- toRational c = toRational (toReal c) + +-- | The square root of a 'Rational' number. +sqrtRat :: Rational -> RealCyclotomic +sqrtRat r + | r >= 0 = RealCyclotomic (Cyc.sqrtRat r) + | otherwise = error "sqrtRational needs a nonnegative argument" + +-- | Sine function with argument in degrees. +sinDeg :: Rational -> RealCyclotomic +sinDeg r = RealCyclotomic (Cyc.sinDeg r) + +-- | Cosine function with argument in degrees. +cosDeg :: Rational -> RealCyclotomic +cosDeg r = RealCyclotomic (Cyc.cosDeg r) + +-- | Sine function with argument in revolutions. +sinRev :: Rational -> RealCyclotomic +sinRev r = RealCyclotomic (Cyc.sinRev r) + +-- | Cosine function with argument in revolutions. +cosRev :: Rational -> RealCyclotomic +cosRev r = RealCyclotomic (Cyc.cosRev r) + +-- | Is the cyclotomic a rational? +isRat :: RealCyclotomic -> Bool +isRat (RealCyclotomic r) = Cyc.isRat r + +-- | Return an exact rational number if possible. +toRat :: RealCyclotomic -> Maybe Rational +toRat (RealCyclotomic r) = Cyc.toRat r + +-- | Export as an inexact real number. +toReal :: RealFloat a => RealCyclotomic -> a +toReal (RealCyclotomic r) = realPart (Cyc.toComplex r) + +-- | The golden ratio, @(1 + √5)/2@. +goldenRatio :: RealCyclotomic +goldenRatio = (1 + sqrtRat 5) / 2 + +-- | Heron's formula for the area of a triangle with +-- side lengths a, b, c. +heron :: Rational -- ^ a + -> Rational -- ^ b + -> Rational -- ^ c + -> RealCyclotomic -- ^ area of triangle +heron a b c + = sqrtRat (s * (s-a) * (s-b) * (s-c)) + where + s = (a + b + c) / 2 diff --git a/test/Properties.hs b/test/Properties.hs new file mode 100644 --- /dev/null +++ b/test/Properties.hs @@ -0,0 +1,193 @@ +{-# OPTIONS_GHC -Wall #-} + +module Main where + +import Data.Complex.Cyclotomic +import Test.Framework + ( defaultMain + , testGroup + ) +import Test.Framework.Providers.HUnit + ( testCase + ) +import Test.Framework.Providers.QuickCheck2 + ( testProperty + ) +import qualified Test.Framework.Providers.SmallCheck as S +import qualified Test.Framework.Providers.API as T +import Test.QuickCheck + ( Gen + , elements + , Arbitrary(..) + , shrinkRealFrac + ) +import Test.HUnit + ( (@?=) + , Assertion + ) +import Data.List + ( nub + ) +import Data.Ratio + ( (%) + ) + +main :: IO () +main = defaultMain tests + +tests :: [T.Test] +tests = [test1a + ,test2b + ,test3b + ,test4b + ,test5 + ,S.withDepth 10 (S.testProperty "SmallCheck prop_square_sqrtRat" prop_square_sqrtRat) + ,qc_square_sqrtRat + ,qc_Gauss + ,qc_dftInv_dft + ,qc_dft_dftInv + ,qc_sum_quadratic_roots + ] + +rationals :: [Rational] +rationals = 0 % 1 : [sign * k % j | n <- [0..], m <- [0..n-1], sign <- [1,-1] + , let k = m + 1, let j = n - m, gcd k j == 1] + +rationalList :: Integer -> [Rational] +rationalList m = nub [n % d | n <- [-m..m], d <- [1..m]] + +test1a :: T.Test +test1a = testGroup "polarRat" [polarRatTest p q | p <- [0..10], q <- [1..10]] + +polarRatAssertion :: Integer -> Integer -> Assertion +polarRatAssertion p q = polarRat 1 (p % q) @?= e q^p + +polarRatTest :: Integer -> Integer -> T.Test +polarRatTest p q = testCase ("polarRat 1 (" ++ show p ++ " % " ++ show q ++ ")") (polarRatAssertion p q) + +test2b :: T.Test +test2b = testGroup "sqrtRat r ^ 2 == r for the following values of r" + [testCase (show r) (sqrtRat r ^ (2::Int) @?= fromRational r) + | r <- take 100 rationals] + +test3b :: T.Test +test3b = testGroup "sqrtRat (r*r) == abs r for the following values of r" + [testCase (show r) (sqrtRat (r*r) @?= fromRational (abs r)) + | r <- take 100 rationals] + +test4b :: T.Test +test4b = testGroup "z * (1 / z) == 1 for the following values of z" + [testCase (show z) (z * (1 / z) @?= 1) + | n <- [1..10], m <- [1..10], let z = e n + e m, z /= 0] + +test5 :: T.Test +test5 = testGroup "Heron's formula" + [testCase "Try Heron" (heron 3 4 5 @?= 6)] + +---------------- +-- Properties -- +---------------- + +prop_square_sqrtRat :: Int -> Bool +prop_square_sqrtRat n = sqrtRat (fromIntegral n) ^ (2::Int) == fromIntegral n + +prop_Gauss :: Integer -> Bool +prop_Gauss n = let nn = 2 * abs n + 1 + in sum [e nn^(j*j `mod` nn) | j <- [1..(nn - 1) `div` 2]] + == if nn `mod` 4 == 1 + then (-1 + sqrtInteger nn) / 2 + else (-1 + i*sqrtInteger nn) / 2 + +prop_dftInv_dft :: [Rational] -> Bool +prop_dftInv_dft rs = dftInv (dft cs) == cs + where cs = map fromRational rs + +prop_dft_dftInv :: [Rational] -> Bool +prop_dft_dftInv rs = dft (dftInv cs) == cs + where cs = map fromRational rs + +prop_sum_quadratic_roots :: (Rational, Rational, Rational) -> Bool +prop_sum_quadratic_roots (a, b, c) + = case rootsQuadEq a b c of + Nothing -> a == 0 + Just (r1,r2) -> r1 + r2 == fromRational (-b / a) + +prop_sum_quadratic_roots_small :: (SmallRational, SmallRational, SmallRational) -> Bool +prop_sum_quadratic_roots_small (SmallRational a, SmallRational b, SmallRational c) + = case rootsQuadEq a b c of + Nothing -> a == 0 + Just (r1,r2) -> r1 + r2 == fromRational (-b / a) + +---------------------- +-- QuickCheck Tests -- +---------------------- + +qc_square_sqrtRat :: T.Test +qc_square_sqrtRat + = T.plusTestOptions (T.TestOptions + {T.topt_seed = Nothing + ,T.topt_maximum_generated_tests = Just 15 + ,T.topt_maximum_unsuitable_generated_tests = Nothing + ,T.topt_maximum_test_size = Just 15 + ,T.topt_maximum_test_depth = Nothing + ,T.topt_timeout = Nothing + }) + $ testProperty "QuickCheck prop_square_sqrtRat" prop_square_sqrtRat + +qc_Gauss :: T.Test +qc_Gauss + = T.plusTestOptions (T.TestOptions + {T.topt_seed = Nothing + ,T.topt_maximum_generated_tests = Nothing + ,T.topt_maximum_unsuitable_generated_tests = Nothing + ,T.topt_maximum_test_size = Nothing + ,T.topt_maximum_test_depth = Nothing + ,T.topt_timeout = Nothing + }) + $ testProperty "QuickCheck prop_Gauss" prop_Gauss + +qc_dftInv_dft :: T.Test +qc_dftInv_dft + = T.plusTestOptions (T.TestOptions + {T.topt_seed = Nothing + ,T.topt_maximum_generated_tests = Just 15 + ,T.topt_maximum_unsuitable_generated_tests = Nothing + ,T.topt_maximum_test_size = Just 30 + ,T.topt_maximum_test_depth = Nothing + ,T.topt_timeout = Nothing + }) + $ testProperty "QuickCheck prop_dftInv_dft" prop_dftInv_dft + +qc_dft_dftInv :: T.Test +qc_dft_dftInv + = T.plusTestOptions (T.TestOptions + {T.topt_seed = Nothing + ,T.topt_maximum_generated_tests = Just 15 + ,T.topt_maximum_unsuitable_generated_tests = Nothing + ,T.topt_maximum_test_size = Just 30 + ,T.topt_maximum_test_depth = Nothing + ,T.topt_timeout = Nothing + }) + $ testProperty "QuickCheck prop_dft_dftInv" prop_dft_dftInv + +qc_sum_quadratic_roots :: T.Test +qc_sum_quadratic_roots + = testProperty "QuickCheck prop_sum_quadratic_roots" prop_sum_quadratic_roots_small + +---------------------- +-- QuickCheck Stuff -- +---------------------- + +data SmallRational = SmallRational Rational + deriving (Show,Ord,Eq) + +smallRationalList :: [SmallRational] +smallRationalList = map SmallRational (rationalList 3) + +smallRationalGen :: Gen SmallRational +smallRationalGen = elements smallRationalList + +instance Arbitrary SmallRational where + arbitrary = smallRationalGen + shrink (SmallRational r) = map SmallRational (shrinkRealFrac r) +