galois-fft (empty) → 0.1.0

raw patch · 9 files changed

+420/−0 lines, 9 filesdep +QuickCheckdep +basedep +criterion

Dependencies added: QuickCheck, base, criterion, elliptic-curve, galois-fft, galois-field, pairing, poly, protolude, quickcheck-instances, tasty, tasty-discover, tasty-hunit, tasty-quickcheck, vector

Files

ChangeLog.md +5/−0
LICENSE +21/−0
README.md +66/−0
bench/Main.hs +16/−0
bench/Poly.hs +38/−0
galois-fft.cabal +93/−0
src/FFT.hs +140/−0
test/Main.hs +1/−0
test/TestFFT.hs +40/−0

+ ChangeLog.md view

@@ -0,0 +1,5 @@+# Change log for galois-fft++## 0.1.0++* Initial release.

+ LICENSE view

@@ -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.

+ README.md view

@@ -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.+```

+ bench/Main.hs view

@@ -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+      ]

+ bench/Poly.hs view

@@ -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+    ]

+ galois-fft.cabal view

@@ -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

+ src/FFT.hs view

@@ -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

+ test/Main.hs view

@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF tasty-discover -optF --tree-display #-}

+ test/TestFFT.hs view

@@ -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