complex: Add initial complex support.
authorNick Bowler <nbowler@draconx.ca>
Sat, 20 Feb 2010 23:26:02 +0000 (20 18:26 -0500)
committerNick Bowler <nbowler@draconx.ca>
Sat, 20 Feb 2010 23:26:02 +0000 (20 18:26 -0500)
Data/Floating/CMath/Complex.hs [new file with mode: 0644]
Data/Floating/Types.hs
altfloat.cabal
complex.c [new file with mode: 0644]

diff --git a/Data/Floating/CMath/Complex.hs b/Data/Floating/CMath/Complex.hs
new file mode 100644 (file)
index 0000000..411b3a7
--- /dev/null
@@ -0,0 +1,161 @@
+{-
+ - Copyright (C) 2010 Nick Bowler.
+ -
+ - License BSD2:  2-clause BSD license.  See LICENSE for full terms.
+ - This is free software: you are free to change and redistribute it.
+ - There is NO WARRANTY, to the extent permitted by law.
+ -}
+
+-- | Bindings to the standard C complex library.  The FFI does not provide a
+-- mechanism to call complex-valued functions, so we create a Storable instance
+-- for Complex CDouble which exploits the fact that, in C, @double _Complex@
+-- has the same representation and alignment requirements as a @double[2]@ with
+-- the first element being the real part and the second being the imaginary
+-- part.  A similar instance is created for Complex CFloat.
+--
+-- Bindings are not provided for the cimag, creal and conj functions as they
+-- wouldn't be useful in a Haskell program.
+{-# LANGUAGE ForeignFunctionInterface #-}
+module Data.Floating.CMath.Complex (
+    -- * Trigonometric functions
+    c_cacos, c_casin, c_catan, c_ccos, c_csin, c_ctan,
+
+    -- * Hyperbolic functions
+    c_cacosh, c_casinh, c_catanh, c_ccosh, c_csinh, c_ctanh,
+
+    -- * Exponential and logarithmic functions
+    c_cexp, c_clog,
+
+    -- * Power and absolute-value functions
+    c_cabs, c_csqrt, c_cpow,
+
+    -- * Manipulation functions
+    c_carg, c_cproj
+) where
+
+import Data.Floating.CMath.Instances
+import Data.Floating.Classes
+import Data.Floating.Types
+
+import Foreign
+import Foreign.C
+
+unwrap :: (Storable a, PrimFloat a) => (Ptr (Complex a) -> IO ())
+    -> Complex a -> Complex a
+unwrap f x = unsafePerformIO . with x $ \p -> f p >> peek p
+
+-- 7.3.5 Trigonometric functions
+foreign import ccall unsafe "cacos_wrap"
+    c_cacos_wrap :: Ptr (Complex CDouble) -> IO ()
+foreign import ccall unsafe "casin_wrap"
+    c_casin_wrap :: Ptr (Complex CDouble) -> IO ()
+foreign import ccall unsafe "catan_wrap"
+    c_catan_wrap :: Ptr (Complex CDouble) -> IO ()
+foreign import ccall unsafe "ccos_wrap"
+    c_ccos_wrap :: Ptr (Complex CDouble) -> IO ()
+foreign import ccall unsafe "csin_wrap"
+    c_csin_wrap :: Ptr (Complex CDouble) -> IO ()
+foreign import ccall unsafe "ctan_wrap"
+    c_ctan_wrap :: Ptr (Complex CDouble) -> IO ()
+
+-- 7.3.6 Hyperbolic functions
+foreign import ccall unsafe "cacosh_wrap"
+    c_cacosh_wrap :: Ptr (Complex CDouble) -> IO ()
+foreign import ccall unsafe "casinh_wrap"
+    c_casinh_wrap :: Ptr (Complex CDouble) -> IO ()
+foreign import ccall unsafe "catanh_wrap"
+    c_catanh_wrap :: Ptr (Complex CDouble) -> IO ()
+foreign import ccall unsafe "ccosh_wrap"
+    c_ccosh_wrap :: Ptr (Complex CDouble) -> IO ()
+foreign import ccall unsafe "csinh_wrap"
+    c_csinh_wrap :: Ptr (Complex CDouble) -> IO ()
+foreign import ccall unsafe "ctanh_wrap"
+    c_ctanh_wrap :: Ptr (Complex CDouble) -> IO ()
+
+-- 7.3.7 Exponential and logarithmic functions
+foreign import ccall unsafe "cexp_wrap"
+    c_cexp_wrap :: Ptr (Complex CDouble) -> IO ()
+foreign import ccall unsafe "clog_wrap"
+    c_clog_wrap :: Ptr (Complex CDouble) -> IO ()
+
+-- 7.3.8 Power and asbolute-value functions
+foreign import ccall unsafe "cabs_wrap"
+    c_cabs_wrap :: Ptr (Complex CDouble) -> IO CDouble
+foreign import ccall unsafe "csqrt_wrap"
+    c_csqrt_wrap :: Ptr (Complex CDouble) -> IO ()
+foreign import ccall unsafe "cpow_wrap"
+    c_cpow_wrap :: Ptr (Complex CDouble) -> Ptr (Complex CDouble) -> IO ()
+
+-- 7.3.9 Manipulation functions
+foreign import ccall unsafe "carg_wrap"
+    c_carg_wrap :: Ptr (Complex CDouble) -> IO CDouble
+foreign import ccall unsafe "cproj_wrap"
+    c_cproj_wrap :: Ptr (Complex CDouble) -> IO ()
+
+c_cacos :: Complex CDouble -> Complex CDouble
+c_cacos = unwrap c_cacos_wrap
+
+c_casin :: Complex CDouble -> Complex CDouble
+c_casin = unwrap c_casin_wrap
+
+c_catan :: Complex CDouble -> Complex CDouble
+c_catan = unwrap c_catan_wrap
+
+c_ccos :: Complex CDouble -> Complex CDouble
+c_ccos = unwrap c_ccos_wrap
+
+c_csin :: Complex CDouble -> Complex CDouble
+c_csin = unwrap c_csin_wrap
+
+c_ctan :: Complex CDouble -> Complex CDouble
+c_ctan = unwrap c_ctan_wrap
+
+c_cacosh :: Complex CDouble -> Complex CDouble
+c_cacosh = unwrap c_cacosh_wrap
+
+c_casinh :: Complex CDouble -> Complex CDouble
+c_casinh = unwrap c_casinh_wrap
+
+c_catanh :: Complex CDouble -> Complex CDouble
+c_catanh = unwrap c_catanh_wrap
+
+c_ccosh :: Complex CDouble -> Complex CDouble
+c_ccosh = unwrap c_ccosh_wrap
+
+c_csinh :: Complex CDouble -> Complex CDouble
+c_csinh = unwrap c_csinh_wrap
+
+c_ctanh :: Complex CDouble -> Complex CDouble
+c_ctanh = unwrap c_ctanh_wrap
+
+c_cexp :: Complex CDouble -> Complex CDouble
+c_cexp = unwrap c_cexp_wrap
+
+c_clog :: Complex CDouble -> Complex CDouble
+c_clog = unwrap c_clog_wrap
+
+c_cabs :: Complex CDouble -> CDouble
+c_cabs x = unsafePerformIO . with x $ c_cabs_wrap
+
+c_cpow :: Complex CDouble -> Complex CDouble -> Complex CDouble
+c_cpow x y = unsafePerformIO . with x $ \px -> with y $ \py -> do
+    c_cpow_wrap px py
+    peek px
+
+c_csqrt :: Complex CDouble -> Complex CDouble
+c_csqrt = unwrap c_csqrt_wrap
+
+c_carg :: Complex CDouble -> CDouble
+c_carg x = unsafePerformIO . with x $ c_carg_wrap
+
+c_cproj :: Complex CDouble -> Complex CDouble
+c_cproj = unwrap c_cproj_wrap
+
+instance (Storable a, PrimFloat a) => Storable (Complex a) where
+    sizeOf    (real :+ _) = 2 * sizeOf real
+    alignment (real :+ _) = alignment real
+    peek ptr  = do
+        [real, imag] <- peekArray 2 (castPtr ptr)
+        return $! real :+ imag
+    poke ptr (real :+ imag) = do
+        pokeArray (castPtr ptr) [real, imag]
index d0274a7..8c84401 100644 (file)
@@ -11,7 +11,7 @@
 {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, OverlappingInstances #-}
 {-# LANGUAGE MagicHash #-}
 module Data.Floating.Types (
-    Double(..), Float(..), FloatConvert(..)
+    Double(..), Float(..), FloatConvert(..), Complex(..)
 ) where
 
 import Prelude hiding (Double, Float)
@@ -20,6 +20,7 @@ import GHC.Integer
 import GHC.Prim
 import Foreign.C
 
+import Data.Floating.Classes
 import Data.Ratio
 
 import Unsafe.Coerce
@@ -32,6 +33,11 @@ data Double = D# Double#
 -- | The Float type. 
 data Float  = F# Float#
 
+infix 6 :+
+-- | Complex numbers.
+data (PrimFloat a) => Complex a = !a :+ !a
+    deriving Eq
+
 -- | Coercion to floating point types.
 class FloatConvert a b where
     -- | Convert to a floating point type.  Conversions from integers and real
index 085e198..9c9609b 100644 (file)
@@ -72,9 +72,10 @@ Library
         Build-Depends: integer
 
     Include-Dirs: .
-    C-Sources: cfloat.c c99-compat.c
+    C-Sources: cfloat.c complex.c c99-compat.c
     Exposed-Modules:
         Data.Floating.CMath,
+        Data.Floating.CMath.Complex,
         Data.Floating.Classes,
         Data.Floating.Types,
         Data.Floating.Types.Double,
diff --git a/complex.c b/complex.c
new file mode 100644 (file)
index 0000000..30f6e7d
--- /dev/null
+++ b/complex.c
@@ -0,0 +1,41 @@
+#include <complex.h>
+
+#define complex_wrap1(T, func) void func ## _wrap(T complex *a) { \
+       *a = func(*a); \
+}
+
+complex_wrap1(double, cacos)
+complex_wrap1(double, casin)
+complex_wrap1(double, catan)
+complex_wrap1(double, ccos)
+complex_wrap1(double, csin)
+complex_wrap1(double, ctan)
+
+complex_wrap1(double, cacosh)
+complex_wrap1(double, casinh)
+complex_wrap1(double, catanh)
+complex_wrap1(double, ccosh)
+complex_wrap1(double, csinh)
+complex_wrap1(double, ctanh)
+
+double cabs_wrap(const double complex *a)
+{
+       return cabs(*a);
+}
+
+complex_wrap1(double, cexp)
+complex_wrap1(double, clog)
+
+complex_wrap1(double, csqrt)
+
+void cpow_wrap(double complex *a, const double complex *b)
+{
+       *a = cpow(*a, *b);
+}
+
+complex_wrap1(double, cproj)
+
+double carg_wrap(const double complex *a)
+{
+       return carg(*a);
+}