diff --git a/ChangeLog.md b/ChangeLog.md
new file mode 100644
--- /dev/null
+++ b/ChangeLog.md
@@ -0,0 +1,5 @@
+# Change log for galois-fft
+
+## 0.1.0
+
+* Initial release.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2019 Adjoint
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,66 @@
+<p align="center">
+<a href="https://www.adjoint.io">
+  <img width="250" src="./.assets/adjoint.png" alt="Adjoint Logo" />
+</a>
+</p>
+
+# galois-fft
+
+Fast Fourier Transforms over finite fields. Provides functionality for
+polynomial evaluation, polynomial interpolation, and computation of Lagrange
+polynomials.
+
+In a finite field F with 2^m elements. We can define a discrete Fourier
+transform by selecting 2^m - 1 roots of unity ω ∈ F.
+
+## Example
+
+
+```haskell
+import Protolude
+
+import Data.Curve.Weierstrass.BN254 (Fr)
+import Data.Pairing.BN254           (getRootOfUnity)
+
+import FFT
+
+k :: Int
+k = 5
+
+polySize :: Int
+polySize = 2^k
+
+leftCoeffs, rightCoeffs :: [Fr]
+leftCoeffs = map fromIntegral [1..polySize]
+rightCoeffs = map fromIntegral (reverse [1..polySize])
+
+main :: IO ()
+main = do
+  print $ interpolate getRootOfUnity leftCoeffs
+  print $ fftMult getRootOfUnity leftCoeffs rightCoeffs
+  pure ()
+```
+
+## License
+
+```
+Copyright (c) 2018-2019 Adjoint Inc.
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE
+OR OTHER DEALINGS IN THE SOFTWARE.
+```
diff --git a/bench/Main.hs b/bench/Main.hs
new file mode 100644
--- /dev/null
+++ b/bench/Main.hs
@@ -0,0 +1,16 @@
+{-# LANGUAGE NoImplicitPrelude #-}
+
+-- To get the benchmarking data, run "stack bench".
+
+module Main where
+
+import Protolude
+
+import Criterion.Main
+
+import Poly
+
+main :: IO ()
+main = defaultMain
+      [ bgroup "Polynomial operations" Poly.benchmarks
+      ]
diff --git a/bench/Poly.hs b/bench/Poly.hs
new file mode 100644
--- /dev/null
+++ b/bench/Poly.hs
@@ -0,0 +1,38 @@
+{-# LANGUAGE NoImplicitPrelude #-}
+
+module Poly (benchmarks) where
+
+import Protolude
+
+import Criterion.Main
+import Data.Curve.Weierstrass.BN254 (Fr)
+import Data.Pairing.BN254           (getRootOfUnity)
+
+import FFT
+
+k :: Int
+k = 5
+
+polySize :: Int
+polySize = 2^k
+
+leftCoeffs, rightCoeffs :: [Fr]
+leftCoeffs = map fromIntegral [1..polySize]
+rightCoeffs = map fromIntegral (reverse [1..polySize])
+
+points :: [(Fr,Fr)]
+points
+  = map (\i -> (getRootOfUnity k ^ i, fromIntegral i))
+        [1..polySize]
+
+fftPoints :: [Fr]
+fftPoints = map snd points
+
+benchmarks :: [Benchmark]
+benchmarks
+  = [ bench "FFT-based multiplication"
+            $ nf (uncurry $ fftMult getRootOfUnity)
+                 (leftCoeffs, rightCoeffs)
+    , bench "FFT-based interpolation"
+            $ nf (interpolate getRootOfUnity) fftPoints
+    ]
diff --git a/galois-fft.cabal b/galois-fft.cabal
new file mode 100644
--- /dev/null
+++ b/galois-fft.cabal
@@ -0,0 +1,93 @@
+cabal-version: 1.12
+
+-- This file has been generated from package.yaml by hpack version 0.31.2.
+--
+-- see: https://github.com/sol/hpack
+--
+-- hash: 501e6cfc24248ef770ef9809bf5fcddd3908419a34a6bfd35be74ee2883c24e4
+
+name:           galois-fft
+version:        0.1.0
+synopsis:       FFTs over finite fields
+description:    Finite field polynomial arithmetic based on fast Fourier transforms
+category:       Cryptography
+homepage:       https://github.com/adjoint-io/galois-fft#readme
+bug-reports:    https://github.com/adjoint-io/galois-fft/issues
+maintainer:     Adjoint Inc (info@adjoint.io)
+license:        MIT
+license-file:   LICENSE
+build-type:     Simple
+extra-source-files:
+    README.md
+    ChangeLog.md
+
+source-repository head
+  type: git
+  location: https://github.com/adjoint-io/galois-fft
+
+library
+  exposed-modules:
+      FFT
+  other-modules:
+      Paths_galois_fft
+  hs-source-dirs:
+      src
+  default-extensions: LambdaCase OverloadedStrings NoImplicitPrelude
+  ghc-options: -freverse-errors -O2 -Wall
+  build-depends:
+      base >=4.10 && <5
+    , elliptic-curve >=0.3 && <0.4
+    , galois-field >=1 && <2
+    , poly >=0.3.2
+    , protolude >=0.2 && <0.3
+    , vector >=0.12 && <0.13
+  default-language: Haskell2010
+
+test-suite fft-tests
+  type: exitcode-stdio-1.0
+  main-is: Main.hs
+  other-modules:
+      TestFFT
+      Paths_galois_fft
+  hs-source-dirs:
+      test
+  default-extensions: LambdaCase OverloadedStrings NoImplicitPrelude
+  ghc-options: -freverse-errors -O2 -Wall -main-is Main
+  build-depends:
+      QuickCheck
+    , base >=4.10 && <5
+    , elliptic-curve >=0.3 && <0.4
+    , galois-fft
+    , galois-field >=1 && <2
+    , pairing
+    , poly >=0.3.2
+    , protolude >=0.2 && <0.3
+    , quickcheck-instances
+    , tasty
+    , tasty-discover
+    , tasty-hunit
+    , tasty-quickcheck
+    , vector >=0.12 && <0.13
+  default-language: Haskell2010
+
+benchmark fft-benchmarks
+  type: exitcode-stdio-1.0
+  main-is: Main.hs
+  other-modules:
+      Poly
+      Paths_galois_fft
+  hs-source-dirs:
+      bench
+  default-extensions: LambdaCase OverloadedStrings NoImplicitPrelude
+  ghc-options: -freverse-errors -O2 -Wall -main-is Main
+  build-depends:
+      base >=4.10 && <5
+    , criterion
+    , elliptic-curve >=0.3 && <0.4
+    , galois-fft
+    , galois-field >=1 && <2
+    , pairing
+    , poly >=0.3.2
+    , protolude >=0.2 && <0.3
+    , vector >=0.12 && <0.13
+  default-language: Haskell2010
diff --git a/src/FFT.hs b/src/FFT.hs
new file mode 100644
--- /dev/null
+++ b/src/FFT.hs
@@ -0,0 +1,140 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TupleSections #-}
+
+-- | Discrete Fourier transforms for polynomial interpolation
+module FFT
+  ( CoeffVec,
+    dftNaive,
+    fft,
+    fftMult,
+    fftTargetPoly,
+    interpolate,
+    inverseDft,
+  )
+where
+
+import Data.Field.Galois (GaloisField, pow)
+import qualified Data.List as List
+import Data.Poly (VPoly, monomial, toPoly)
+import Data.Vector (fromList)
+import Protolude
+
+-- | Polynomial represented as a coefficient vector, little-endian
+type CoeffVec f = [f]
+
+-- | Discrete Fourier transform. Can be interpreted as some polynomial
+-- evaluated at certain roots of unity. (In our case the length of
+-- these lists will be a power of two.)
+type DFT f = [f]
+
+-- | Evaluate a polynomial given by its coefficient vector
+evalPoly :: Num f => CoeffVec f -> f -> f
+evalPoly coeffs x = foldr (\c rest -> c + x * rest) 0 coeffs
+
+-- | Naive discrete Fourier transformation performed by evaluating the
+-- polynomial at the appropriate roots of unity.
+dftNaive ::
+  Num f =>
+  -- | principal 2^k-th root of unity
+  f ->
+  -- | polynomial coefficients, length should be 2^k for
+  -- some k
+  CoeffVec f ->
+  DFT f
+dftNaive omega_n as = map (\i -> evalPoly as (omega_n ^ i)) [0 .. length as - 1]
+
+-- | Split a list into a list containing the odd-numbered and one with
+-- the even-numbered elements.
+split :: [a] -> ([a], [a])
+split = foldr (\a (r1, r2) -> (a : r2, r1)) ([], [])
+
+-- | Calculate ceiling of log base 2 of an integer.
+log2 :: Int -> Int
+log2 x = floorLog + correction
+  where
+    floorLog = finiteBitSize x - 1 - countLeadingZeros x
+    correction =
+      if countTrailingZeros x < floorLog
+        then 1
+        else 0
+
+-- | Fast Fourier transformation.
+fft ::
+  GaloisField k =>
+  -- | function that gives for input n the principal (2^n)-th root of unity
+  (Int -> k) ->
+  -- | length should be n
+  CoeffVec k ->
+  DFT k
+fft omega_n as =
+  case length as of
+    1 -> as
+    n ->
+      let (as0, as1) = split as
+          y0 = fft omega_n as0
+          y1 = fft omega_n as1
+          omegas = map (pow (omega_n (log2 n))) [0 .. n]
+       in combine y0 y1 omegas
+  where
+    combine y0 y1 omegas =
+      (\xs -> map fst xs ++ map snd xs)
+        $ map (\(yk0, yk1, currentOmega) -> (yk0 + currentOmega * yk1, yk0 - currentOmega * yk1))
+        $ List.zip3 y0 y1 omegas
+
+-- | Inverse discrete Fourier transformation, uses FFT.
+inverseDft :: GaloisField k => (Int -> k) -> DFT k -> CoeffVec k
+inverseDft primRootsUnity dft =
+  let n = fromIntegral . length $ dft
+   in map (/ n) $
+        fft (recip . primRootsUnity) dft
+
+-- | Append minimal amount of zeroes until the list has a length which
+-- is a power of two.
+padToNearestPowerOfTwo :: Num f => [f] -> [f]
+padToNearestPowerOfTwo [] = []
+padToNearestPowerOfTwo xs = padToNearestPowerOfTwoOf (length xs) xs
+
+-- | Given n, append zeroes until the list has length 2^n.
+padToNearestPowerOfTwoOf ::
+  Num f =>
+  -- | n
+  Int ->
+  -- | list which should have length <= 2^n
+  [f] ->
+  -- | list which will have length 2^n
+  [f]
+padToNearestPowerOfTwoOf i xs = xs ++ replicate padLength 0
+  where
+    padLength = nearestPowerOfTwo - length xs
+    nearestPowerOfTwo = bit $ log2 i
+
+-- | Create a polynomial that goes through the given values.
+interpolate :: GaloisField k => (Int -> k) -> [k] -> VPoly k
+interpolate primRoots pts = toPoly . fromList $ inverseDft primRoots (padToNearestPowerOfTwo pts)
+
+-- | Multiply polynomials using FFT
+fftMult :: GaloisField k => (Int -> k) -> CoeffVec k -> CoeffVec k -> CoeffVec k
+fftMult primRoots l r = inverseDft primRoots $ zipWith (*) dftL dftR
+  where
+    n = 2 * max (length l) (length r)
+    paddedDft x = fft primRoots (padToNearestPowerOfTwoOf n x)
+    dftL = paddedDft l
+    dftR = paddedDft r
+
+-- XXX make this actually go fast
+-- polyWithZeroesAt
+--    :: Fractional f
+--    => (Int -> f)
+--    -> [f]
+--    -> CoeffVec f
+-- polyWithZeroesAt primRoots
+--   = foldl' (fftMult primRoots) [1]
+--     . map (\xcoord -> [-xcoord, 1])
+
+-- XXX make this actually use FFT mult
+fftTargetPoly :: GaloisField k => (Int -> k) -> Int -> VPoly k
+fftTargetPoly primRoots numRoots =
+  foldl' (*) (monomial 0 1) ((\i -> toPoly . fromList $ [- pow omega i, 1]) <$> [0 .. 2 ^ k - 1 :: Integer])
+  where
+    k = log2 numRoots
+    omega = primRoots k
diff --git a/test/Main.hs b/test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Main.hs
@@ -0,0 +1,1 @@
+{-# OPTIONS_GHC -F -pgmF tasty-discover -optF --tree-display #-}
diff --git a/test/TestFFT.hs b/test/TestFFT.hs
new file mode 100644
--- /dev/null
+++ b/test/TestFFT.hs
@@ -0,0 +1,40 @@
+{-# LANGUAGE DataKinds #-}
+
+module TestFFT where
+
+import Data.Curve.Weierstrass.BN254 (Fr)
+import Data.Field.Galois
+import Data.Pairing.BN254 (getRootOfUnity)
+import FFT
+import Protolude
+import Test.Tasty.HUnit
+import Test.Tasty.QuickCheck
+
+data TestPoly f = TestPoly (CoeffVec f)
+  deriving (Show)
+
+testExp :: Int
+testExp = 8
+
+instance Arbitrary f => Arbitrary (TestPoly f) where
+  arbitrary = TestPoly <$> vectorOf (2 ^ testExp) arbitrary
+
+omega :: Fr
+omega = getRootOfUnity testExp
+
+unit_omega_is_primitive_root :: Assertion
+unit_omega_is_primitive_root =
+  assertBool "omega is not a correct primitive root of unity" $
+    isRootOfUnity (toU' omega :: RootsOfUnity 256 Fr)
+
+-- We are choosing 256 because there are n=2^8=256 roots of unity
+-- in this test case
+
+prop_fftCorrect :: TestPoly Fr -> Bool
+prop_fftCorrect (TestPoly coeffs) =
+  dftNaive omega coeffs == fft getRootOfUnity coeffs
+
+prop_inverseDftCorrect :: TestPoly Fr -> Bool
+prop_inverseDftCorrect (TestPoly coeffs) =
+  inverseDft getRootOfUnity (fft getRootOfUnity coeffs) == coeffs
+    && fft getRootOfUnity (inverseDft getRootOfUnity coeffs) == coeffs