diff --git a/CHANGES.md b/CHANGES.md new file mode 100644 --- /dev/null +++ b/CHANGES.md @@ -0,0 +1,6 @@ +Changelog for lol-cpp project +================================ + +0.0.0.1 +---- + * Initial split from lol package. Now requires GHC >= 8.0. diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP.hs b/Crypto/Lol/Cyclotomic/Tensor/CPP.hs new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP.hs @@ -0,0 +1,454 @@ +{-| +Module : Crypto.Lol.Cyclotomic.Tensor.CPP +Description : Wrapper for a C++ implementation of the 'Tensor' interface. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX + +Wrapper for a C++ implementation of the 'Tensor' interface. +-} + +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE InstanceSigs #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE RebindableSyntax #-} +{-# LANGUAGE RoleAnnotations #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UndecidableInstances #-} + +module Crypto.Lol.Cyclotomic.Tensor.CPP (CT) where + +import Algebra.Additive as Additive (C) +import Algebra.Module as Module (C) +import Algebra.ZeroTestable as ZeroTestable (C) + +import Control.Applicative hiding ((*>)) +import Control.Arrow ((***)) +import Control.DeepSeq +import Control.Monad.Except +import Control.Monad.Identity (Identity (..), runIdentity) +import Control.Monad.Random +import Control.Monad.Trans as T (lift) + +import Data.Coerce +import Data.Constraint hiding ((***)) +import Data.Int +import Data.Maybe +import Data.Traversable as T +import Data.Vector.Generic as V (fromList, toList, unzip) +import Data.Vector.Storable as SV (Vector, convert, foldl', + fromList, generate, + length, map, replicate, + replicateM, thaw, thaw, toList, + unsafeFreeze, + unsafeWith, zipWith, (!)) +import Data.Vector.Storable.Mutable as SM hiding (replicate) + +import Foreign.Marshal.Utils (with) +import Foreign.Ptr + +import Crypto.Lol.CRTrans +import Crypto.Lol.Cyclotomic.Tensor +import Crypto.Lol.Cyclotomic.Tensor.CPP.Backend +import Crypto.Lol.Cyclotomic.Tensor.CPP.Extension +import Crypto.Lol.Cyclotomic.Tensor.CPP.Instances () +import Crypto.Lol.GaussRandom +import Crypto.Lol.Prelude as LP hiding + (replicate, + unzip, zip) +import Crypto.Lol.Types.FiniteField +import Crypto.Lol.Types.IZipVector +import Crypto.Lol.Types.Proto +import Crypto.Lol.Utils.ShowType + +import Data.Foldable as F + +import System.IO.Unsafe (unsafePerformIO) + +-- | Newtype wrapper around a Vector. +newtype CT' (m :: Factored) r = CT' { unCT :: Vector r } + deriving (Show, Eq, NFData) + +-- the first argument, though phantom, affects representation +type role CT' representational nominal + +-- GADT wrapper that distinguishes between Unbox and unrestricted +-- element types + +-- | An implementation of 'Tensor' backed by C++ code. +data CT (m :: Factored) r where + CT :: Storable r => CT' m r -> CT m r + ZV :: IZipVector m r -> CT m r + +deriving instance Show r => Show (CT m r) + +instance Show (ArgType CT) where + show _ = "CT" + +instance Eq r => Eq (CT m r) where + (ZV x) == (ZV y) = x == y + (CT x) == (CT y) = x == y + x@(CT _) == y = x == toCT y + y == x@(CT _) = x == toCT y + +instance (Protoable (IZipVector m r), Fact m, Storable r) => Protoable (CT m r) where + type ProtoType (CT m r) = ProtoType (IZipVector m r) + + toProto x@(CT _) = toProto $ toZV x + toProto (ZV x) = toProto x + + fromProto x = toCT <$> ZV <$> fromProto x + +toCT :: (Storable r) => CT m r -> CT m r +toCT v@(CT _) = v +toCT (ZV v) = CT $ zvToCT' v + +toZV :: (Fact m) => CT m r -> CT m r +toZV (CT (CT' v)) = ZV $ fromMaybe (error "toZV: internal error") $ + iZipVector $ convert v +toZV v@(ZV _) = v + +zvToCT' :: forall m r . (Storable r) => IZipVector m r -> CT' m r +zvToCT' v = coerce (convert $ unIZipVector v :: Vector r) + +wrap :: (Storable s, Storable r) => (CT' l s -> CT' m r) -> (CT l s -> CT m r) +{-# INLINABLE wrap #-} +wrap f (CT v) = CT $ f v +wrap f (ZV v) = CT $ f $ zvToCT' v + +wrapM :: (Storable s, Storable r, Monad mon) => (CT' l s -> mon (CT' m r)) + -> (CT l s -> mon (CT m r)) +{-# INLINABLE wrapM #-} +wrapM f (CT v) = CT <$> f v +wrapM f (ZV v) = CT <$> f (zvToCT' v) + +-- convert an CT' *twace* signature to Tagged one +type family Tw (r :: *) :: * where + Tw (CT' m' r -> CT' m r) = Tagged '(m,m') (Vector r -> Vector r) + Tw (Maybe (CT' m' r -> CT' m r)) = TaggedT '(m,m') Maybe (Vector r -> Vector r) + +type family Em r where + Em (CT' m r -> CT' m' r) = Tagged '(m,m') (Vector r -> Vector r) + Em (Maybe (CT' m r -> CT' m' r)) = TaggedT '(m,m') Maybe (Vector r -> Vector r) + + +---------- NUMERIC PRELUDE INSTANCES ---------- + +-- CJP: Additive, Ring are not necessary when we use zipWithT +-- EAC: This has performance implications for the CT backend, +-- which used a (very fast) C function for (*) and (+) +instance (Additive r, Storable r, Fact m) + => Additive.C (CT m r) where + (CT (CT' a)) + (CT (CT' b)) = CT $ CT' $ SV.zipWith (+) a b + a + b = toCT a + toCT b + negate (CT (CT' a)) = CT $ CT' $ SV.map negate a -- EAC: This probably should be converted to C code + negate a = negate $ toCT a + + zero = CT $ repl zero + +{- +instance (Fact m, Ring r, Storable r, Dispatch r) + => Ring.C (CT m r) where + (CT a@(CT' _)) * (CT b@(CT' _)) = CT $ (untag $ cZipDispatch dmul) a b + + fromInteger = CT . repl . fromInteger +-} + +instance (ZeroTestable r, Storable r) + => ZeroTestable.C (CT m r) where + --{-# INLINABLE isZero #-} + isZero (CT (CT' a)) = SV.foldl' (\ b x -> b && isZero x) True a + isZero (ZV v) = isZero v + +instance (GFCtx fp d, Fact m, Additive (CT m fp)) + => Module.C (GF fp d) (CT m fp) where + + r *> v = case v of + CT (CT' arr) -> CT $ CT' $ SV.fromList $ unCoeffs $ r *> Coeffs $ SV.toList arr + ZV zv -> ZV $ fromJust $ iZipVector $ V.fromList $ unCoeffs $ r *> Coeffs $ V.toList $ unIZipVector zv + +---------- Category-theoretic instances ---------- + +instance Fact m => Functor (CT m) where + -- Functor instance is implied by Applicative laws + fmap f x = pure f <*> x + +instance Fact m => Applicative (CT m) where + pure = ZV . pure + + (ZV f) <*> (ZV a) = ZV (f <*> a) + f@(ZV _) <*> v@(CT _) = f <*> toZV v + +instance Fact m => Foldable (CT m) where + -- Foldable instance is implied by Traversable + foldMap = foldMapDefault + +instance Fact m => Traversable (CT m) where + traverse f r@(CT _) = T.traverse f $ toZV r + traverse f (ZV v) = ZV <$> T.traverse f v + +instance Tensor CT where + + type TElt CT r = (Storable r, Dispatch r) + + entailIndexT = tag $ Sub Dict + entailEqT = tag $ Sub Dict + entailZTT = tag $ Sub Dict + -- entailRingT = tag $ Sub Dict + entailNFDataT = tag $ Sub Dict + entailRandomT = tag $ Sub Dict + entailShowT = tag $ Sub Dict + entailModuleT = tag $ Sub Dict + + scalarPow = CT . scalarPow' -- Vector code + + l = wrap $ basicDispatch dl + lInv = wrap $ basicDispatch dlinv + + mulGPow = wrap $ basicDispatch dmulgpow + mulGDec = wrap $ basicDispatch dmulgdec + + divGPow = wrapM $ dispatchGInv dginvpow + divGDec = wrapM $ dispatchGInv dginvdec + + crtFuncs = (,,,,) <$> + return (CT . repl) <*> + (wrap . untag (cZipDispatch dmul) <$> gCRT) <*> + (wrap . untag (cZipDispatch dmul) <$> gInvCRT) <*> + (wrap <$> untagT ctCRT) <*> + (wrap <$> untagT ctCRTInv) + + twacePowDec = wrap $ runIdentity $ coerceTw twacePowDec' + embedPow = wrap $ runIdentity $ coerceEm embedPow' + embedDec = wrap $ runIdentity $ coerceEm embedDec' + + tGaussianDec v = CT <$> cDispatchGaussian v + --tGaussianDec v = CT <$> coerceT' (gaussianDec v) + + -- we do not wrap this function because (currently) it can only be called on lifted types + gSqNormDec (CT v) = untag gSqNormDec' v + gSqNormDec (ZV v) = gSqNormDec (CT $ zvToCT' v) + + crtExtFuncs = (,) <$> (wrap <$> coerceTw twaceCRT') + <*> (wrap <$> coerceEm embedCRT') + + coeffs = wrapM $ coerceCoeffs coeffs' + + powBasisPow = (CT <$>) <$> coerceBasis powBasisPow' + + crtSetDec = (CT <$>) <$> coerceBasis crtSetDec' + + fmapT f = wrap $ coerce (SV.map f) + + zipWithT f v1' v2' = + let (CT (CT' v1)) = toCT v1' + (CT (CT' v2)) = toCT v2' + in CT $ CT' $ SV.zipWith f v1 v2 + + unzipT v = + let (CT (CT' x)) = toCT v + in (CT . CT') *** (CT . CT') $ unzip x + + {-# INLINABLE entailIndexT #-} + {-# INLINABLE entailEqT #-} + {-# INLINABLE entailZTT #-} + {-# INLINABLE entailNFDataT #-} + {-# INLINABLE entailRandomT #-} + {-# INLINABLE entailShowT #-} + {-# INLINABLE scalarPow #-} + {-# INLINABLE l #-} + {-# INLINABLE lInv #-} + {-# INLINABLE mulGPow #-} + {-# INLINABLE mulGDec #-} + {-# INLINABLE divGPow #-} + {-# INLINABLE divGDec #-} + {-# INLINABLE crtFuncs #-} + {-# INLINABLE twacePowDec #-} + {-# INLINABLE embedPow #-} + {-# INLINABLE embedDec #-} + {-# INLINABLE tGaussianDec #-} + {-# INLINABLE gSqNormDec #-} + {-# INLINE crtExtFuncs #-} + {-# INLINABLE coeffs #-} + {-# INLINABLE powBasisPow #-} + {-# INLINABLE crtSetDec #-} + {-# INLINABLE fmapT #-} + {-# INLINE zipWithT #-} + {-# INLINE unzipT #-} + +coerceTw :: (Functor mon) => TaggedT '(m, m') mon (Vector r -> Vector r) -> mon (CT' m' r -> CT' m r) +coerceTw = (coerce <$>) . untagT + +coerceEm :: (Functor mon) => TaggedT '(m, m') mon (Vector r -> Vector r) -> mon (CT' m r -> CT' m' r) +coerceEm = (coerce <$>) . untagT + +-- | Useful coersion for defining @coeffs@ in the @Tensor@ +-- interface. Using 'coerce' alone is insufficient for type inference. +coerceCoeffs :: Tagged '(m,m') (Vector r -> [Vector r]) -> CT' m' r -> [CT' m r] +coerceCoeffs = coerce + +-- | Useful coersion for defining @powBasisPow@ and @crtSetDec@ in the @Tensor@ +-- interface. Using 'coerce' alone is insufficient for type inference. +coerceBasis :: Tagged '(m,m') [Vector r] -> Tagged m [CT' m' r] +coerceBasis = coerce + +dispatchGInv :: forall m r . (Storable r, Fact m) + => (Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO Int16) + -> CT' m r -> Maybe (CT' m r) +dispatchGInv f = + let factors = proxy (marshalFactors <$> ppsFact) (Proxy::Proxy m) + totm = proxy (fromIntegral <$> totientFact) (Proxy::Proxy m) + numFacts = fromIntegral $ SV.length factors + in \(CT' x) -> unsafePerformIO $ do + yout <- SV.thaw x + ret <- SM.unsafeWith yout (\pout -> + SV.unsafeWith factors (\pfac -> + f pout totm pfac numFacts)) + if ret /= 0 + then Just . CT' <$> unsafeFreeze yout + else return Nothing + +withBasicArgs :: forall m r . (Fact m, Storable r) + => (Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ()) + -> CT' m r -> IO (CT' m r) +withBasicArgs f = + let factors = proxy (marshalFactors <$> ppsFact) (Proxy::Proxy m) + totm = proxy (fromIntegral <$> totientFact) (Proxy::Proxy m) + numFacts = fromIntegral $ SV.length factors + in \(CT' x) -> do + yout <- SV.thaw x + SM.unsafeWith yout (\pout -> + SV.unsafeWith factors (\pfac -> + f pout totm pfac numFacts)) + CT' <$> unsafeFreeze yout + +basicDispatch :: (Storable r, Fact m) + => (Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ()) + -> CT' m r -> CT' m r +basicDispatch f = unsafePerformIO . withBasicArgs f + +gSqNormDec' :: (Storable r, Fact m, Dispatch r) + => Tagged m (CT' m r -> r) +gSqNormDec' = return $ (!0) . unCT . unsafePerformIO . withBasicArgs dnorm + +ctCRT :: (Storable r, CRTrans mon r, Dispatch r, Fact m) + => TaggedT m mon (CT' m r -> CT' m r) +ctCRT = do + ru' <- ru + return $ \x -> unsafePerformIO $ + withPtrArray ru' (flip withBasicArgs x . dcrt) + +-- CTensor CRT^(-1) functions take inverse rus +ctCRTInv :: (Storable r, CRTrans mon r, Dispatch r, Fact m) + => TaggedT m mon (CT' m r -> CT' m r) +ctCRTInv = do + mhatInv <- snd <$> crtInfo + ruinv' <- ruInv + return $ \x -> unsafePerformIO $ + withPtrArray ruinv' (\ruptr -> with mhatInv (flip withBasicArgs x . dcrtinv ruptr)) + +cZipDispatch :: (Storable r, Fact m) + => (Ptr r -> Ptr r -> Int64 -> IO ()) + -> Tagged m (CT' m r -> CT' m r -> CT' m r) +cZipDispatch f = do -- in Tagged m + totm <- fromIntegral <$> totientFact + return $ coerce $ \a b -> unsafePerformIO $ do + yout <- SV.thaw a + SM.unsafeWith yout (\pout -> + SV.unsafeWith b (\pin -> + f pout pin totm)) + unsafeFreeze yout + +cDispatchGaussian :: forall m r var rnd . + (Storable r, Transcendental r, Dispatch r, Ord r, + Fact m, ToRational var, Random r, MonadRandom rnd) + => var -> rnd (CT' m r) +cDispatchGaussian var = flip proxyT (Proxy::Proxy m) $ do -- in TaggedT m rnd + -- get rus for (Complex r) + -- takes ru (not ruInv) to match RT + ruinv' <- mapTaggedT (return . fromMaybe (error "complexGaussianRoots")) ru + totm <- pureT totientFact + mval <- pureT valueFact + rad <- pureT radicalFact + yin <- T.lift $ realGaussians (var * fromIntegral (mval `div` rad)) totm + return $ unsafePerformIO $ + withPtrArray ruinv' (\ruptr -> withBasicArgs (dgaussdec ruptr) (CT' yin)) + +instance (Storable r, Random r, Fact m) => Random (CT' m r) where + --{-# INLINABLE random #-} + random = runRand $ replM (liftRand random) + + randomR = error "randomR nonsensical for CT'" + +instance (Storable r, Random (CT' m r)) => Random (CT m r) where + --{-# INLINABLE random #-} + random = runRand $ CT <$> liftRand random + + randomR = error "randomR nonsensical for CT" + +instance (NFData r) => NFData (CT m r) where + rnf (CT v) = rnf v + rnf (ZV v) = rnf v + +repl :: forall m r . (Fact m, Storable r) => r -> CT' m r +repl = let n = proxy totientFact (Proxy::Proxy m) + in coerce . SV.replicate n + +replM :: forall m r mon . (Fact m, Storable r, Monad mon) + => mon r -> mon (CT' m r) +replM = let n = proxy totientFact (Proxy::Proxy m) + in fmap coerce . SV.replicateM n + +scalarPow' :: forall m r . (Fact m, Additive r, Storable r) => r -> CT' m r +-- constant-term coefficient is first entry wrt powerful basis +scalarPow' = + let n = proxy totientFact (Proxy::Proxy m) + in \r -> CT' $ generate n (\i -> if i == 0 then r else zero) + +ru, ruInv :: (CRTrans mon r, Fact m, Storable r) + => TaggedT m mon [Vector r] +ru = do + mval <- pureT valueFact + wPow <- fst <$> crtInfo + LP.map + (\(p,e) -> do + let pp = p^e + pow = mval `div` pp + generate pp (wPow . (*pow))) <$> + pureT ppsFact + +ruInv = do + mval <- pureT valueFact + wPow <- fst <$> crtInfo + LP.map + (\(p,e) -> do + let pp = p^e + pow = mval `div` pp + generate pp (\i -> wPow $ -i*pow)) <$> + pureT ppsFact + +wrapVector :: forall mon m r . (Monad mon, Fact m, Ring r, Storable r) + => TaggedT m mon (Kron r) -> mon (CT' m r) +wrapVector v = do + vmat <- proxyT v (Proxy::Proxy m) + let n = proxy totientFact (Proxy::Proxy m) + return $ CT' $ generate n (flip (indexK vmat) 0) + +gCRT, gInvCRT :: (Storable r, CRTrans mon r, Fact m) + => mon (CT' m r) +gCRT = wrapVector gCRTK +gInvCRT = wrapVector gInvCRTK diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/Backend.hs b/Crypto/Lol/Cyclotomic/Tensor/CPP/Backend.hs new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/Backend.hs @@ -0,0 +1,337 @@ +{-| +Module : Crypto.Lol.Cyclotomic.Tensor.CPP.Backend +Description : Transforms Haskell types into C counterparts. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX + +This module contains the functions to transform Haskell types into their +C counterpart, and to transform polymorphic Haskell functions into C funtion +calls in a type-safe way. +-} + +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UndecidableInstances #-} + +{-# OPTIONS_GHC -fno-warn-unticked-promoted-constructors #-} + +module Crypto.Lol.Cyclotomic.Tensor.CPP.Backend +( Dispatch +, dcrt, dcrtinv +, dgaussdec +, dl, dlinv +, dnorm +, dmulgpow, dmulgdec +, dginvpow, dginvdec +, dmul +, marshalFactors +, CPP +, withArray, withPtrArray +) where + +import Crypto.Lol.Prelude as LP (Complex, PP, Proxy (..), Tagged, + map, mapM_, proxy, tag) +import Crypto.Lol.Reflects +import Crypto.Lol.Types.Unsafe.RRq +import Crypto.Lol.Types.Unsafe.ZqBasic + +import Data.Int +import Data.Vector.Storable as SV (Vector, fromList, + unsafeToForeignPtr0) +import Data.Vector.Storable.Internal (getPtr) + +import Foreign.ForeignPtr (touchForeignPtr) +import Foreign.Marshal.Array (withArray) +import Foreign.Marshal.Utils (with) +import Foreign.Ptr (Ptr, castPtr, plusPtr) +import Foreign.Storable (Storable (..)) + +import GHC.TypeLits -- for error message + +-- | Convert a list of prime powers to a suitable C representation. +marshalFactors :: [PP] -> Vector CPP +marshalFactors = SV.fromList . LP.map (\(p,e) -> (fromIntegral p, fromIntegral e)) + +-- http://stackoverflow.com/questions/6517387/vector-vector-foo-ptr-ptr-foo-io-a-io-a +-- | Evaluates a C function that takes an "a** ptr" on a list of Vectors. +withPtrArray :: (Storable a) => [Vector a] -> (Ptr (Ptr a) -> IO b) -> IO b +withPtrArray v f = do + let vs = LP.map SV.unsafeToForeignPtr0 v + ptrV = LP.map (\(fp,_) -> getPtr fp) vs + res <- withArray ptrV f + LP.mapM_ (\(fp,_) -> touchForeignPtr fp) vs + return res + +-- Note: These types need to be the same, otherwise something goes wrong on the C end... +-- | C representation of a prime power. +type CPP = (Int16, Int16) + +instance (Storable a, Storable b) + => Storable (a,b) where + sizeOf _ = sizeOf (undefined :: a) + sizeOf (undefined :: b) + alignment _ = max (alignment (undefined :: a)) (alignment (undefined :: b)) + peek p = do + a <- peek (castPtr p :: Ptr a) + b <- peek (castPtr (plusPtr p (sizeOf a)) :: Ptr b) + return (a,b) + poke p (a,b) = do + poke (castPtr p :: Ptr a) a + poke (castPtr (plusPtr p (sizeOf a)) :: Ptr b) b + +data ZqB64D -- for type safety purposes +data ComplexD +data DoubleD +data Int64D +data RRqD + +type family CTypeOf x where + CTypeOf (a,b) = EqCType a b (CTypeOf a) (CTypeOf b) + CTypeOf (ZqBasic (q :: k) Int64) = ZqB64D + CTypeOf Double = DoubleD + CTypeOf Int64 = Int64D + CTypeOf (Complex Double) = ComplexD + CTypeOf (RRq (q :: k) Double) = RRqD + + -- EAC: See #12237 and #11990 + CTypeOf (ZqBasic (q :: k) i) = TypeError (Text "Unsupported C type: " :<>: ShowType (ZqBasic q i) :$$: Text "Use Int64 as the base ring") + CTypeOf (Complex i) = TypeError (Text "Unsupported C type: " :<>: ShowType (Complex i) :$$: Text "Use Double as the base ring") + CTypeOf (RRq (q :: k) i) = TypeError (Text "Unsupported C type: " :<>: ShowType (RRq q i) :$$: Text "Use Double as the base ring") + CTypeOf a = TypeError (Text "Unsupported C type: " :<>: ShowType a) + +type family EqCType a b c d where + EqCType a b ZqB64D ZqB64D = ZqB64D + EqCType a b RRqD RRqD = RRqD + EqCType a b ComplexD ComplexD = ComplexD + EqCType a b c c = TypeError (Text "Cannot call C code on a tuple of type " :<>: ShowType a) + EqCType a b c d = TypeError (Text "You are trying to use CTensor on a tuple," :<>: + Text " but the tuple contains two different C types: " :$$: + ShowType a :<>: Text " and " :<>: ShowType b) + +-- returns the modulus as a nested list of moduli +class (Tuple a) => ZqTuple a where + type ModPairs a + getModuli :: Tagged a (ModPairs a) + +instance (Reflects q Int64) => ZqTuple (ZqBasic q Int64) where + type ModPairs (ZqBasic q Int64) = Int64 + getModuli = tag $ proxy value (Proxy::Proxy q) + +instance (Reflects q r, RealFrac r) => ZqTuple (RRq q r) where + type ModPairs (RRq q r) = Int64 + getModuli = tag $ round (proxy value (Proxy::Proxy q) :: r) + +instance (ZqTuple a, ZqTuple b) => ZqTuple (a, b) where + type ModPairs (a,b) = (ModPairs a, ModPairs b) + getModuli = + let as = proxy getModuli (Proxy::Proxy a) + bs = proxy getModuli (Proxy :: Proxy b) + in tag (as,bs) + +-- counts components in a nested tuple +class Tuple a where + numComponents :: Tagged a Int16 + +instance {-# Overlappable #-} Tuple a where + numComponents = tag 1 + +instance (Tuple a, Tuple b) => Tuple (a,b) where + numComponents = tag $ proxy numComponents (Proxy::Proxy a) + proxy numComponents (Proxy::Proxy b) + +-- | Single-argument synonym for @Dispatch'@. +type Dispatch r = (Dispatch' (CTypeOf r) r) + +-- | Class to safely match Haskell types with the appropriate C function. +class (repr ~ CTypeOf r) => Dispatch' repr r where + -- | Equivalent to 'Tensor's @crt@. + dcrt :: Ptr (Ptr r) -> Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO () + -- | Equivalent to 'Tensor's @crtInv@. + dcrtinv :: Ptr (Ptr r) -> Ptr r -> Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO () + -- | Equivalent to 'Tensor's @tGaussianDec@. + dgaussdec :: Ptr (Ptr (Complex r)) -> Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO () + -- | Equivalent to 'Tensor's @l@. + dl :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO () + -- | Equivalent to 'Tensor's @lInv@. + dlinv :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO () + -- | Equivalent to 'Tensor's @gSqNormDec@. + dnorm :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO () + -- | Equivalent to 'Tensor's @mulGPow@. + dmulgpow :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO () + -- | Equivalent to 'Tensor's @mulGDec@. + dmulgdec :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO () + -- | Equivalent to 'Tensor's @divGPow@. + dginvpow :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO Int16 + -- | Equivalent to 'Tensor's @divGDec@. + dginvdec :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO Int16 + -- | Equivalent to @zipWith (*)@ + dmul :: Ptr r -> Ptr r -> Int64 -> IO () + +instance (ZqTuple r, Storable (ModPairs r), CTypeOf r ~ RRqD) + => Dispatch' RRqD r where + dcrt = error "cannot call CT CRT on type RRq" + dcrtinv = error "cannot call CT CRTInv on type RRq" + dl = error "cannot call CT L on type RRq (though you probably should be able to)" + dlinv = error "cannot call CT LInv on type RRq (though you probably should be able to)" + dnorm = error "cannto call CT normSq on type RRq" + dmulgpow = error "cannot call CT mulGPow on type RRq" + dmulgdec = error "cannot call CT mulGDec on type RRq" + dginvpow = error "cannot call CT divGPow on type RRq" + dginvdec = error "cannot call CT divGDec on type RRq" + dmul = error "cannot call CT mul on type RRq" + dgaussdec = error "cannot call CT gaussianDec on type RRq" + +instance (ZqTuple r, Storable (ModPairs r), CTypeOf r ~ ZqB64D) + => Dispatch' ZqB64D r where + dcrt ruptr pout totm pfac numFacts = + let qs = proxy getModuli (Proxy::Proxy r) + numPairs = proxy numComponents (Proxy::Proxy r) + in with qs $ \qsptr -> + tensorCRTRq numPairs (castPtr pout) totm pfac numFacts (castPtr ruptr) (castPtr qsptr) + dcrtinv ruptr minv pout totm pfac numFacts = + let qs = proxy getModuli (Proxy::Proxy r) + numPairs = proxy numComponents (Proxy::Proxy r) + in with qs $ \qsptr -> + tensorCRTInvRq numPairs (castPtr pout) totm pfac numFacts (castPtr ruptr) (castPtr minv) (castPtr qsptr) + dl pout totm pfac numFacts = + let qs = proxy getModuli (Proxy::Proxy r) + numPairs = proxy numComponents (Proxy::Proxy r) + in with qs $ \qsptr -> + tensorLRq numPairs (castPtr pout) totm pfac numFacts (castPtr qsptr) + dlinv pout totm pfac numFacts = + let qs = proxy getModuli (Proxy::Proxy r) + numPairs = proxy numComponents (Proxy::Proxy r) + in with qs $ \qsptr -> + tensorLInvRq numPairs (castPtr pout) totm pfac numFacts (castPtr qsptr) + dnorm = error "cannot call CT normSq on type ZqBasic" + dmulgpow pout totm pfac numFacts = + let qs = proxy getModuli (Proxy::Proxy r) + numPairs = proxy numComponents (Proxy::Proxy r) + in with qs $ \qsptr -> + tensorGPowRq numPairs (castPtr pout) totm pfac numFacts (castPtr qsptr) + dmulgdec pout totm pfac numFacts = + let qs = proxy getModuli (Proxy::Proxy r) + numPairs = proxy numComponents (Proxy::Proxy r) + in with qs $ \qsptr -> + tensorGDecRq numPairs (castPtr pout) totm pfac numFacts (castPtr qsptr) + dginvpow pout totm pfac numFacts = + let qs = proxy getModuli (Proxy::Proxy r) + numPairs = proxy numComponents (Proxy::Proxy r) + in with qs $ \qsptr -> + tensorGInvPowRq numPairs (castPtr pout) totm pfac numFacts (castPtr qsptr) + dginvdec pout totm pfac numFacts = + let qs = proxy getModuli (Proxy::Proxy r) + numPairs = proxy numComponents (Proxy::Proxy r) + in with qs $ \qsptr -> + tensorGInvDecRq numPairs (castPtr pout) totm pfac numFacts (castPtr qsptr) + dmul aout bout totm = + let qs = proxy getModuli (Proxy::Proxy r) + numPairs = proxy numComponents (Proxy::Proxy r) + in with qs $ \qsptr -> + mulRq numPairs (castPtr aout) (castPtr bout) totm (castPtr qsptr) + dgaussdec = error "cannot call CT gaussianDec on type ZqBasic" + +-- products of Complex correspond to CRTExt of a Zq product +instance (Tuple r, CTypeOf r ~ ComplexD) => Dispatch' ComplexD r where + dcrt ruptr pout totm pfac numFacts = + tensorCRTC (proxy numComponents (Proxy::Proxy r)) (castPtr pout) totm pfac numFacts (castPtr ruptr) + dcrtinv ruptr minv pout totm pfac numFacts = + tensorCRTInvC (proxy numComponents (Proxy::Proxy r)) (castPtr pout) totm pfac numFacts (castPtr ruptr) (castPtr minv) + dl pout = + tensorLC (proxy numComponents (Proxy::Proxy r)) (castPtr pout) + dlinv pout = + tensorLInvC (proxy numComponents (Proxy::Proxy r)) (castPtr pout) + dnorm = error "cannot call CT normSq on type Complex Double" + dmulgpow pout = + tensorGPowC (proxy numComponents (Proxy::Proxy r)) (castPtr pout) + dmulgdec pout = + tensorGDecC (proxy numComponents (Proxy::Proxy r)) (castPtr pout) + dginvpow pout = + tensorGInvPowC (proxy numComponents (Proxy::Proxy r)) (castPtr pout) + dginvdec pout = + tensorGInvDecC (proxy numComponents (Proxy::Proxy r)) (castPtr pout) + dmul aout bout = + mulC (proxy numComponents (Proxy::Proxy r)) (castPtr aout) (castPtr bout) + dgaussdec = error "cannot call CT gaussianDec on type Comple Double" + +-- no support for products of Double +instance Dispatch' DoubleD Double where + dcrt = error "cannot call CT Crt on type Double" + dcrtinv = error "cannot call CT CrtInv on type Double" + dl pout = + tensorLDouble 1 (castPtr pout) + dlinv pout = + tensorLInvDouble 1 (castPtr pout) + dnorm pout = tensorNormSqD 1 (castPtr pout) + dmulgpow = error "cannot call CT mulGPow on type Double" + dmulgdec = error "cannot call CT mulGDec on type Double" + dginvpow = error "cannot call CT divGPow on type Double" + dginvdec = error "cannot call CT divGDec on type Double" + dmul = error "cannot call CT (*) on type Double" + dgaussdec ruptr pout totm pfac numFacts = + tensorGaussianDec 1 (castPtr pout) totm pfac numFacts (castPtr ruptr) + +-- no support for products of Z +instance Dispatch' Int64D Int64 where + dcrt = error "cannot call CT Crt on type Int64" + dcrtinv = error "cannot call CT CrtInv on type Int64" + dl pout = + tensorLR 1 (castPtr pout) + dlinv pout = + tensorLInvR 1 (castPtr pout) + dnorm pout = + tensorNormSqR 1 (castPtr pout) + dmulgpow pout = + tensorGPowR 1 (castPtr pout) + dmulgdec pout = + tensorGDecR 1 (castPtr pout) + dginvpow pout = + tensorGInvPowR 1 (castPtr pout) + dginvdec pout = + tensorGInvDecR 1 (castPtr pout) + dmul = error "cannot call CT (*) on type Int64" + dgaussdec = error "cannot call CT gaussianDec on type Int64" + +foreign import ccall unsafe "tensorLR" tensorLR :: Int16 -> Ptr Int64 -> Int64 -> Ptr CPP -> Int16 -> IO () +foreign import ccall unsafe "tensorLInvR" tensorLInvR :: Int16 -> Ptr Int64 -> Int64 -> Ptr CPP -> Int16 -> IO () +foreign import ccall unsafe "tensorLRq" tensorLRq :: Int16 -> Ptr (ZqBasic q Int64) -> Int64 -> Ptr CPP -> Int16 -> Ptr Int64 -> IO () +foreign import ccall unsafe "tensorLInvRq" tensorLInvRq :: Int16 -> Ptr (ZqBasic q Int64) -> Int64 -> Ptr CPP -> Int16 -> Ptr Int64 -> IO () +foreign import ccall unsafe "tensorLDouble" tensorLDouble :: Int16 -> Ptr Double -> Int64 -> Ptr CPP -> Int16 -> IO () +foreign import ccall unsafe "tensorLInvDouble" tensorLInvDouble :: Int16 -> Ptr Double -> Int64 -> Ptr CPP -> Int16 -> IO () +foreign import ccall unsafe "tensorLC" tensorLC :: Int16 -> Ptr (Complex Double) -> Int64 -> Ptr CPP -> Int16 -> IO () +foreign import ccall unsafe "tensorLInvC" tensorLInvC :: Int16 -> Ptr (Complex Double) -> Int64 -> Ptr CPP -> Int16 -> IO () + +foreign import ccall unsafe "tensorNormSqR" tensorNormSqR :: Int16 -> Ptr Int64 -> Int64 -> Ptr CPP -> Int16 -> IO () +foreign import ccall unsafe "tensorNormSqD" tensorNormSqD :: Int16 -> Ptr Double -> Int64 -> Ptr CPP -> Int16 -> IO () + +foreign import ccall unsafe "tensorGPowR" tensorGPowR :: Int16 -> Ptr Int64 -> Int64 -> Ptr CPP -> Int16 -> IO () +foreign import ccall unsafe "tensorGPowRq" tensorGPowRq :: Int16 -> Ptr (ZqBasic q Int64) -> Int64 -> Ptr CPP -> Int16 -> Ptr Int64 -> IO () +foreign import ccall unsafe "tensorGPowC" tensorGPowC :: Int16 -> Ptr (Complex Double) -> Int64 -> Ptr CPP -> Int16 -> IO () +foreign import ccall unsafe "tensorGDecR" tensorGDecR :: Int16 -> Ptr Int64 -> Int64 -> Ptr CPP -> Int16 -> IO () +foreign import ccall unsafe "tensorGDecRq" tensorGDecRq :: Int16 -> Ptr (ZqBasic q Int64) -> Int64 -> Ptr CPP -> Int16 -> Ptr Int64 -> IO () +foreign import ccall unsafe "tensorGDecC" tensorGDecC :: Int16 -> Ptr (Complex Double) -> Int64 -> Ptr CPP -> Int16 -> IO () +foreign import ccall unsafe "tensorGInvPowR" tensorGInvPowR :: Int16 -> Ptr Int64 -> Int64 -> Ptr CPP -> Int16 -> IO Int16 +foreign import ccall unsafe "tensorGInvPowRq" tensorGInvPowRq :: Int16 -> Ptr (ZqBasic q Int64) -> Int64 -> Ptr CPP -> Int16 -> Ptr Int64 -> IO Int16 +foreign import ccall unsafe "tensorGInvPowC" tensorGInvPowC :: Int16 -> Ptr (Complex Double) -> Int64 -> Ptr CPP -> Int16 -> IO Int16 +foreign import ccall unsafe "tensorGInvDecR" tensorGInvDecR :: Int16 -> Ptr Int64 -> Int64 -> Ptr CPP -> Int16 -> IO Int16 +foreign import ccall unsafe "tensorGInvDecRq" tensorGInvDecRq :: Int16 -> Ptr (ZqBasic q Int64) -> Int64 -> Ptr CPP -> Int16 -> Ptr Int64 -> IO Int16 +foreign import ccall unsafe "tensorGInvDecC" tensorGInvDecC :: Int16 -> Ptr (Complex Double) -> Int64 -> Ptr CPP -> Int16 -> IO Int16 + +foreign import ccall unsafe "tensorCRTRq" tensorCRTRq :: Int16 -> Ptr (ZqBasic q Int64) -> Int64 -> Ptr CPP -> Int16 -> Ptr (Ptr (ZqBasic q Int64)) -> Ptr Int64 -> IO () +foreign import ccall unsafe "tensorCRTC" tensorCRTC :: Int16 -> Ptr (Complex Double) -> Int64 -> Ptr CPP -> Int16 -> Ptr (Ptr (Complex Double)) -> IO () +foreign import ccall unsafe "tensorCRTInvRq" tensorCRTInvRq :: Int16 -> Ptr (ZqBasic q Int64) -> Int64 -> Ptr CPP -> Int16 -> Ptr (Ptr (ZqBasic q Int64)) -> Ptr (ZqBasic q Int64) -> Ptr Int64 -> IO () +foreign import ccall unsafe "tensorCRTInvC" tensorCRTInvC :: Int16 -> Ptr (Complex Double) -> Int64 -> Ptr CPP -> Int16 -> Ptr (Ptr (Complex Double)) -> Ptr (Complex Double) -> IO () + +foreign import ccall unsafe "tensorGaussianDec" tensorGaussianDec :: Int16 -> Ptr Double -> Int64 -> Ptr CPP -> Int16 -> Ptr (Ptr (Complex Double)) -> IO () + +foreign import ccall unsafe "mulRq" mulRq :: Int16 -> Ptr (ZqBasic q Int64) -> Ptr (ZqBasic q Int64) -> Int64 -> Ptr Int64 -> IO () +foreign import ccall unsafe "mulC" mulC :: Int16 -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int64 -> IO () diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/Extension.hs b/Crypto/Lol/Cyclotomic/Tensor/CPP/Extension.hs new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/Extension.hs @@ -0,0 +1,164 @@ +{-| +Module : Crypto.Lol.Cyclotomic.Tensor.CPP.Extension +Description : Embedding/twacing in various bases for CPP. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX + +CPP Tensor-specific functions for embedding/twacing in various bases. +-} + +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RebindableSyntax #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} + +module Crypto.Lol.Cyclotomic.Tensor.CPP.Extension +( embedPow', embedDec', embedCRT' +, twacePowDec', twaceCRT' +, coeffs', powBasisPow' +, crtSetDec' +, backpermute' +) where + +import Crypto.Lol.CRTrans +import Crypto.Lol.Cyclotomic.Tensor as T +import Crypto.Lol.Prelude as LP hiding (lift, null) +import Crypto.Lol.Types.FiniteField +import Crypto.Lol.Types.ZmStar + + +import Control.Applicative hiding (empty) +import Control.Monad.Trans (lift) + +import Data.Maybe +import Data.Reflection (reify) +import qualified Data.Vector as V +import Data.Vector.Storable as SV +import qualified Data.Vector.Unboxed as U + + +-- | /O(n)/ Yield the vector obtained by replacing each element @i@ of the +-- index vector by @xs'!'i@. This is equivalent to @'map' (xs'!') is@ but is +-- often much more efficient. +-- +-- > backpermute <0,3,2,3,1,0> = +backpermute' :: (Storable a) => + U.Vector Int -- ^ @is@ index vector (of length @n@) + -> Vector a -- ^ @xs@ value vector + -> Vector a +{-# INLINABLE backpermute' #-} +backpermute' is v = generate (U.length is) (\i -> v ! (is U.! i)) + +embedPow', embedDec' :: (Additive r, Storable r, m `Divides` m') + => Tagged '(m, m') (Vector r -> Vector r) +{-# INLINABLE embedPow' #-} +{-# INLINABLE embedDec' #-} +-- | Embeds an vector in the powerful basis of the the mth cyclotomic ring +-- to an vector in the powerful basis of the m'th cyclotomic ring when @m | m'@ +embedPow' = (\indices arr -> generate (U.length indices) $ \idx -> + let (j0,j1) = indices U.! idx + in if j0 == 0 + then arr ! j1 + else zero) <$> baseIndicesPow +-- | Embeds an vector in the decoding basis of the the mth cyclotomic ring +-- to an vector in the decoding basis of the m'th cyclotomic ring when @m | m'@ +embedDec' = (\indices arr -> generate (U.length indices) + (\idx -> maybe LP.zero + (\(sh,b) -> if b then negate (arr ! sh) else arr ! sh) + (indices U.! idx))) <$> baseIndicesDec + +-- | Embeds an vector in the CRT basis of the the mth cyclotomic ring +-- to an vector in the CRT basis of the m'th cyclotomic ring when @m | m'@ +embedCRT' :: forall mon m m' r . (CRTrans mon r, Storable r, m `Divides` m') + => TaggedT '(m, m') mon (Vector r -> Vector r) +embedCRT' = + (lift (proxyT crtInfo (Proxy::Proxy m') :: mon (CRTInfo r))) >> + (pureT $ backpermute' <$> baseIndicesCRT) + +-- | maps a vector in the powerful/decoding basis, representing an +-- O_m' element, to a vector of arrays representing O_m elements in +-- the same type of basis +coeffs' :: (Storable r, m `Divides` m') + => Tagged '(m, m') (Vector r -> [Vector r]) +coeffs' = flip (\x -> V.toList . V.map (`backpermute'` x)) + <$> extIndicesCoeffs + +-- | The "tweaked trace" function in either the powerful or decoding +-- basis of the m'th cyclotomic ring to the mth cyclotomic ring when +-- @m | m'@. +twacePowDec' :: forall m m' r . (Storable r, m `Divides` m') + => Tagged '(m, m') (Vector r -> Vector r) +{-# INLINABLE twacePowDec' #-} +twacePowDec' = backpermute' <$> extIndicesPowDec + +kronToVec :: forall mon m r . (Monad mon, Fact m, Ring r, Storable r) + => TaggedT m mon (Kron r) -> TaggedT m mon (Vector r) +kronToVec v = do + vmat <- v + let n = proxy totientFact (Proxy::Proxy m) + return $ generate n (flip (indexK vmat) 0) + +twaceCRT' :: forall mon m m' r . + (Storable r, CRTrans mon r, m `Divides` m') + => TaggedT '(m, m') mon (Vector r -> Vector r) +{-# INLINE twaceCRT' #-} +twaceCRT' = tagT $ do + g' <- proxyT (kronToVec gCRTK) (Proxy::Proxy m') + gInv <- proxyT (kronToVec gInvCRTK) (Proxy::Proxy m) + embed <- proxyT embedCRT' (Proxy::Proxy '(m,m')) + indices <- pure $ proxy extIndicesCRT (Proxy::Proxy '(m,m')) + (_, m'hatinv) <- proxyT crtInfo (Proxy::Proxy m') + let phi = proxy totientFact (Proxy::Proxy m) + phi' = proxy totientFact (Proxy::Proxy m') + mhat = fromIntegral $ proxy valueHatFact (Proxy::Proxy m) + hatRatioInv = m'hatinv * mhat + reltot = phi' `div` phi + -- tweak = mhat * g' / (m'hat * g) + tweak = SV.map (* hatRatioInv) $ SV.zipWith (*) (embed gInv) g' + return $ \ arr -> -- take true trace after mul-by-tweak + let v = backpermute' indices (SV.zipWith (*) tweak arr) + in generate phi $ \i -> foldl1' (+) $ SV.unsafeSlice (i*reltot) reltot v + +-- | The powerful extension basis, wrt the powerful basis. +-- Outputs a list of vectors in O_m' that are an O_m basis for O_m' +powBasisPow' :: forall m m' r . (m `Divides` m', Ring r, SV.Storable r) + => Tagged '(m, m') [SV.Vector r] +powBasisPow' = do + (_, phi, phi', _) <- indexInfo + idxs <- baseIndicesPow + return $ LP.map (\k -> generate phi' $ \j -> + let (j0,j1) = idxs U.! j + in if j0==k && j1==0 then one else zero) + [0..phi' `div` phi - 1] + +-- | A list of vectors representing the mod-p CRT set of the +-- extension O_m'/O_m +crtSetDec' :: forall m m' fp . + (m `Divides` m', PrimeField fp, Coprime (PToF (CharOf fp)) m', SV.Storable fp) + => Tagged '(m, m') [SV.Vector fp] +crtSetDec' = + let m'p = Proxy :: Proxy m' + p = proxy valuePrime (Proxy::Proxy (CharOf fp)) + phi = proxy totientFact m'p + d = proxy (order p) m'p + h :: Int = proxy valueHatFact m'p + hinv = recip $ fromIntegral h + in reify d $ \(_::Proxy d) -> do + let twCRTs' :: Kron (GF fp d) + = fromMaybe (error "internal error: crtSetDec': twCRTs") $ proxyT twCRTs m'p + zmsToIdx = proxy T.zmsToIndexFact m'p + elt j i = indexK twCRTs' j (zmsToIdx i) + trace' = trace :: GF fp d -> fp -- to avoid recomputing powTraces + cosets <- partitionCosets p + return $ LP.map (\is -> generate phi + (\j -> hinv * trace' + (LP.sum $ LP.map (elt j) is))) cosets diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/Instances.hs b/Crypto/Lol/Cyclotomic/Tensor/CPP/Instances.hs new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/Instances.hs @@ -0,0 +1,30 @@ +{-| +Module : Crypto.Lol.Cyclotomic.Tensor.CPP.Instances +Description : CPP Tensor-specific instances. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX + +CPP Tensor-specific instances. +-} + +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE StandaloneDeriving #-} + +module Crypto.Lol.Cyclotomic.Tensor.CPP.Instances where + +-- EAC: Do not import Crypto.Lol.Types, because it exports an IrreduciblePoly +-- instance which screw with GHC. Probably #10338. +import Crypto.Lol.Types.Unsafe.Complex +import Crypto.Lol.Types.Unsafe.RRq +import Crypto.Lol.Types.Unsafe.ZqBasic + +import Foreign.Storable + +deriving instance (Storable a) => Storable (Complex a) +deriving instance (Storable r) => Storable (RRq q r) +deriving instance (Storable i) => Storable (ZqBasic q i) diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/common.cpp b/Crypto/Lol/Cyclotomic/Tensor/CPP/common.cpp new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/common.cpp @@ -0,0 +1,27 @@ +/* +Module : common.cpp +Description : Shared functions and data. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX +*/ + +#include "types.h" + +hInt_t Zq::q; // should be in zq.cpp; here due to GHC #12152 + +hDim_t ipow(hDim_t base, hShort_t exp) +{ + hDim_t result = 1; + while (exp) { + if (exp & 1) { + result *= base; + } + exp >>= 1; + base *= base; + } + return result; +} diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/common.h b/Crypto/Lol/Cyclotomic/Tensor/CPP/common.h new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/common.h @@ -0,0 +1,20 @@ +/* +Module : common.h +Description : Shared functions and data. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX +*/ + +#ifndef COMMON_H_ +#define COMMON_H_ + +#include "types.h" + +// calculates base ** exp +hDim_t ipow(hDim_t base, hShort_t exp); + +#endif /* COMMON_H_ */ diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/crt.cpp b/Crypto/Lol/Cyclotomic/Tensor/CPP/crt.cpp new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/crt.cpp @@ -0,0 +1,598 @@ +/* +Module : crt.cpp +Description : Chinese remainder transform. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX +*/ + +#include "types.h" +#include "tensor.h" +#include "common.h" + +// If this macro is modified, make sure to update all functions below with +// cascading if/else statments so that temp space is allocated when necessary +// (i.e., for all primes >= DFTP_GENERIC_SIZE) +#define DFTP_GENERIC_SIZE 11 + +hDim_t bitrev (PrimeExponent pe, hDim_t j) { + hShort_t e; + hDim_t p = pe.prime; + hDim_t tempj = j; + hDim_t acc = 0; + + for(e = pe.exponent-1; e >= 0; e--) { + div_t qr = div(tempj,p); + acc += qr.rem * ipow(p,e); + tempj = qr.quot; + } + return acc; +} + +template void crtTwiddle (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, + PrimeExponent pe, ring* ru) +{ + hDim_t p = pe.prime; + hShort_t e = pe.exponent; + + pe.exponent -= 1; // used for an argument to bitrev + + if(p == 2) { + hDim_t mprime = 1<<(e-1); + hDim_t blockDim = rts*mprime; // size of block in block diagonal tensor matrix + + for(hDim_t i0 = 1; i0 < mprime; i0++) { // loops over i/(p-1) for i = 0..(m'-1), we can skip i0 = 0 + hDim_t temp2 = i0*rts; + ring twid = ru[bitrev(pe, i0)*tupSize]; + + for(hDim_t blockIdx = 0; blockIdx < lts; blockIdx++) { + hDim_t temp3 = blockIdx*blockDim + temp2; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + hDim_t idx = (temp3 + modOffset)*tupSize; + y[idx] *= twid; + } + } + } + } + else { // This loop is faster, probably due to the division in the loop above. + // cilk also slows it down + hDim_t mprime = ipow(p,e-1); + hDim_t blockDim = rts*(p-1)*mprime; // size of block in block diagonal tensor matrix + + for(hDim_t i0 = 1; i0 < mprime; i0++) { // loops over i/(p-1) for i = 0..(m'-1), we can skip i0 = 0 + hDim_t temp1 = i0*(p-1); + for(hDim_t i1 = 0; i1 < (p-1); i1++) { // loops over i%(p-1) for i = 0..(m'-1) + hDim_t temp2 = (temp1+i1)*rts; + ring twid = ru[bitrev(pe, i0)*(i1+1)*tupSize]; + + for(hDim_t blockIdx = 0; blockIdx < lts; blockIdx++) { + hDim_t temp3 = blockIdx*blockDim + temp2; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + hDim_t idx = (temp3 + modOffset)*tupSize; + y[idx] *= twid; + } + } + } + } + } +} + +// dim is power of p +template void dftTwiddle (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, + PrimeExponent pe, hDim_t dim, hDim_t rustride, ring* ru) +{ + hDim_t idx; + hDim_t p = pe.prime; + + pe.exponent -= 1; // used for an argument to bitrev + + if(p == 2) { + hDim_t mprime = dim>>1; // divides evenly + hDim_t temp1 = rts*dim; // for use in computing [modified] tensorOffset + for(hDim_t i0 = 1; i0 < mprime; i0++) { // loops over i/p for i = 0..(dim-1), but we skip i0=0 + hDim_t temp3 = rts*(i0*p+1); + ring twid = ru[bitrev(pe,i0)*rustride*tupSize]; + + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1 + temp3; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + idx = (temp2 + modOffset)*tupSize; + y[idx] *= twid; + } + } + } + } + else { + hDim_t mprime = dim/p; // divides evenly + hDim_t temp1 = rts*dim; // for use in computing [modified] tensorOffset + for(hDim_t i0 = 1; i0 < mprime; i0++) { // loops over i/p for i = 0..(dim-1), but we skip i0=0 + for(hDim_t i1 = 1; i1 < p; i1++) { // loops over i%p for i = 0..(dim-1), but we skip i1=0 + hDim_t temp3 = rts*(i0*p+i1); + ring twid = ru[bitrev(pe,i0)*i1*rustride*tupSize]; + + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1 + temp3; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + idx = (temp2 + modOffset)*tupSize; + y[idx] *= twid; + } + } + } + } + } +} + +//implied length of ru is rustride*p +//implied length of tempSpace is p, if p is not a special case +// temp is allowed to be NULL if p < DFTP_GENERIC_SIZE +template void dftp (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, + hDim_t p, hDim_t rustride, ring* ru, ring* tempSpace) +{ + hDim_t tensorOffset; + + if(p == 2) { + hDim_t temp1 = rts<<1; + + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + ring u = y[tensorOffset*tupSize]; + ring t = y[(tensorOffset+rts)*tupSize]; + y[tensorOffset*tupSize] = u + t; + y[(tensorOffset+rts)*tupSize] = u - t; + } + } + } + else if(p == 3) { + ring ru1 = ru[rustride*tupSize]; + ring ru2 = ru[(rustride<<1)*tupSize]; + hDim_t temp1 = rts*3; + + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + ring y1, y2, y3; + y1 = y[tensorOffset*tupSize]; + y2 = y[(tensorOffset+rts)*tupSize]; + y3 = y[(tensorOffset+(rts<<1))*tupSize]; + //q is <32 bits, so we can do 3 additions without overflow + y[tensorOffset*tupSize] += (y2 + y3); + y[(tensorOffset+rts)*tupSize] = y1 + (ru1*y2) + (ru2*y3); + y[(tensorOffset+(rts<<1))*tupSize] = y1 + (ru2*y2) + (ru1*y3); + } + } + } + else if(p == 5) { + hDim_t temp1 = rts*5; + ring ru1 = ru[rustride*tupSize]; + ring ru2 = ru[(rustride<<1)*tupSize]; + ring ru3 = ru[(rustride*3)*tupSize]; + ring ru4 = ru[(rustride<<2)*tupSize]; + + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + ring y1, y2, y3, y4, y5; + y1 = y[tensorOffset*tupSize]; + y2 = y[(tensorOffset+rts)*tupSize]; + y3 = y[(tensorOffset+(rts<<1))*tupSize]; + y4 = y[(tensorOffset+3*rts)*tupSize]; + y5 = y[(tensorOffset+(rts<<2))*tupSize]; + y[tensorOffset*tupSize] += y2 + y3 + y4 + y5; + y[(tensorOffset+rts)*tupSize] = y1 + (ru1*y2) + (ru2*y3) + (ru3*y4) + (ru4*y5); + y[(tensorOffset+(rts<<1))*tupSize] = y1 + (ru2*y2) + (ru4*y3) + (ru1*y4) + (ru3*y5); + y[(tensorOffset+rts*3)*tupSize] = y1 + (ru3*y2) + (ru1*y3) + (ru4*y4) + (ru2*y5); + y[(tensorOffset+(rts<<2))*tupSize] = y1 + (ru4*y2) + (ru3*y3) + (ru2*y4) + (ru1*y5); + } + } + } + else if(p == 7) { + hDim_t temp1 = rts*7; + ring ru1 = ru[rustride*tupSize]; + ring ru2 = ru[(rustride<<1)*tupSize]; + ring ru3 = ru[(rustride*3)*tupSize]; + ring ru4 = ru[(rustride<<2)*tupSize]; + ring ru5 = ru[(rustride*5)*tupSize]; + ring ru6 = ru[(rustride*6)*tupSize]; + + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + ring y1, y2, y3, y4, y5, y6, y7; + y1 = y[tensorOffset*tupSize]; + y2 = y[(tensorOffset+rts)*tupSize]; + y3 = y[(tensorOffset+(rts<<1))*tupSize]; + y4 = y[(tensorOffset+3*rts)*tupSize]; + y5 = y[(tensorOffset+(rts<<2))*tupSize]; + y6 = y[(tensorOffset+rts*5)*tupSize]; + y7 = y[(tensorOffset+rts*6)*tupSize]; + y[tensorOffset*tupSize] += y2 + y3 + y4 + y5 + y6 + y7; + y[(tensorOffset+rts)*tupSize] = y1 + (ru1*y2) + (ru2*y3) + (ru3*y4) + (ru4*y5) + (ru5*y6) + (ru6*y7); + y[(tensorOffset+(rts<<1))*tupSize] = y1 + (ru2*y2) + (ru4*y3) + (ru6*y4) + (ru1*y5) + (ru3*y6) + (ru5*y7); + y[(tensorOffset+rts*3)*tupSize] = y1 + (ru3*y2) + (ru6*y3) + (ru2*y4) + (ru5*y5) + (ru1*y6) + (ru4*y7); + y[(tensorOffset+(rts<<2))*tupSize] = y1 + (ru4*y2) + (ru1*y3) + (ru5*y4) + (ru2*y5) + (ru6*y6) + (ru3*y7); + y[(tensorOffset+rts*5)*tupSize] = y1 + (ru5*y2) + (ru3*y3) + (ru1*y4) + (ru6*y5) + (ru4*y6) + (ru2*y7); + y[(tensorOffset+rts*6)*tupSize] = y1 + (ru6*y2) + (ru5*y3) + (ru4*y4) + (ru3*y5) + (ru2*y6) + (ru1*y7); + } + } + } + else { + hDim_t temp1 = rts*p; + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + for(hDim_t row = 0; row < p; row++) { + tempSpace[row] = 0; + //p is small (<< 30 bits), so we can do p additions of mod-q values without overflow + for(hDim_t col = 0; col < p; col++) { + tempSpace[row] += (y[(tensorOffset+col*rts)*tupSize]*ru[((col*row) % p)*rustride*tupSize]); + } + } + + for(hDim_t row = 0; row < p; row++) { + y[(tensorOffset+rts*row)*tupSize] = tempSpace[row]; + } + } + } + } +} + +template void crtp (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, + hDim_t p, hDim_t rustride, ring* ru) +{ + hDim_t tensorOffset; + if(p == 2) { + return; + } + else if(p == 3) { + hDim_t temp1 = rts*2; + ring ru1 = ru[rustride*tupSize]; + ring ru2 = ru[(rustride<<1)*tupSize]; + + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + ring y1, y2; + y1 = y[tensorOffset*tupSize]; + y2 = y[(tensorOffset+rts)*tupSize]; + y[tensorOffset*tupSize] += (ru1*y2); + y[(tensorOffset+rts)*tupSize] = y1 + (ru2*y2); + } + } + } + else if(p == 5) { + hDim_t temp1 = rts*4; + ring ru1 = ru[rustride*tupSize]; + ring ru2 = ru[(rustride<<1)*tupSize]; + ring ru3 = ru[(rustride*3)*tupSize]; + ring ru4 = ru[(rustride<<2)*tupSize]; + + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + ring y1, y2, y3, y4; + y1 = y[tensorOffset*tupSize]; + y2 = y[(tensorOffset+rts)*tupSize]; + y3 = y[(tensorOffset+(rts<<1))*tupSize]; + y4 = y[(tensorOffset+3*rts)*tupSize]; + + y[tensorOffset*tupSize] += ((ru1*y2) + (ru2*y3) + (ru3*y4)); + y[(tensorOffset+rts)*tupSize] = y1 + (ru2*y2) + (ru4*y3) + (ru1*y4); + y[(tensorOffset+(rts<<1))*tupSize] = y1 + (ru3*y2) + (ru1*y3) + (ru4*y4); + y[(tensorOffset+rts*3)*tupSize] = y1 + (ru4*y2) + (ru3*y3) + (ru2*y4); + } + } + } + else if(p == 7) { + hDim_t temp1 = rts*6; + ring ru1 = ru[rustride*tupSize]; + ring ru2 = ru[(rustride<<1)*tupSize]; + ring ru3 = ru[(rustride*3)*tupSize]; + ring ru4 = ru[(rustride<<2)*tupSize]; + ring ru5 = ru[(rustride*5)*tupSize]; + ring ru6 = ru[(rustride*6)*tupSize]; + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + ring y1, y2, y3, y4, y5, y6; + y1 = y[tensorOffset*tupSize]; + y2 = y[(tensorOffset+rts)*tupSize]; + y3 = y[(tensorOffset+(rts<<1))*tupSize]; + y4 = y[(tensorOffset+3*rts)*tupSize]; + y5 = y[(tensorOffset+(rts<<2))*tupSize]; + y6 = y[(tensorOffset+rts*5)*tupSize]; + y[tensorOffset*tupSize] += ((ru1*y2) + (ru2*y3) + (ru3*y4) + (ru4*y5) + (ru5*y6)); + y[(tensorOffset+rts)*tupSize] = y1 + (ru2*y2) + (ru4*y3) + (ru6*y4) + (ru1*y5) + (ru3*y6); + y[(tensorOffset+(rts<<1))*tupSize] = y1 + (ru3*y2) + (ru6*y3) + (ru2*y4) + (ru5*y5) + (ru1*y6); + y[(tensorOffset+rts*3)*tupSize] = y1 + (ru4*y2) + (ru1*y3) + (ru5*y4) + (ru2*y5) + (ru6*y6); + y[(tensorOffset+(rts<<2))*tupSize] = y1 + (ru5*y2) + (ru3*y3) + (ru1*y4) + (ru6*y5) + (ru4*y6); + y[(tensorOffset+rts*5)*tupSize] = y1 + (ru6*y2) + (ru5*y3) + (ru4*y4) + (ru3*y5) + (ru2*y6); + } + } + } + else { + ring* tempSpace = (ring*)malloc((p-1)*sizeof(ring)); + hDim_t temp1 = rts*(p-1); + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + + for(hDim_t row = 1; row < p; row++) { + tempSpace[row-1] = 0; + for(hDim_t col = 0; col < p-1; col++) { + tempSpace[row-1] += (y[(tensorOffset+col*rts)*tupSize]*ru[((col*row) % p)*rustride*tupSize]); + } + } + + for(hDim_t row = 0; row < p-1; row++) { + y[(tensorOffset+rts*row)*tupSize] = tempSpace[row]; + } + } + } + free(tempSpace); + } +} + +//takes inverse rus +template void crtpinv (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, + hDim_t p, hDim_t rustride, ring* ruinv) +{ + hDim_t tensorOffset; + if(p == 2) { + return; + } + else if(p == 3) { + hDim_t temp1 = rts*2; + ring ru1 = ruinv[rustride*tupSize]; + ring ru2 = ruinv[(rustride<<1)*tupSize]; + + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + ring y1, y2, shift; + y1 = y[tensorOffset*tupSize]; + y2 = y[(tensorOffset+rts)*tupSize]; + + shift = (ru2*y1) + (ru1*y2); + + y[tensorOffset*tupSize] += y2 - shift; + y[(tensorOffset+rts)*tupSize] = (ru1*y1) + (ru2*y2) - shift; + } + } + } + else if(p == 5) { + hDim_t temp1 = rts*4; + ring ru1 = ruinv[rustride*tupSize]; + ring ru2 = ruinv[(rustride<<1)*tupSize]; + ring ru3 = ruinv[(rustride*3)*tupSize]; + ring ru4 = ruinv[(rustride<<2)*tupSize]; + + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + ring y1, y2, y3, y4, shift; + y1 = y[tensorOffset*tupSize]; + y2 = y[(tensorOffset+rts)*tupSize]; + y3 = y[(tensorOffset+(rts<<1))*tupSize]; + y4 = y[(tensorOffset+3*rts)*tupSize]; + + shift = (ru4*y1) + (ru3*y2) + (ru2*y3) + (ru1*y4); + + y[tensorOffset*tupSize] += y2 + y3 + y4 - shift; + y[(tensorOffset+rts)*tupSize] = (ru1*y1) + (ru2*y2) + (ru3*y3) + (ru4*y4) - shift; + y[(tensorOffset+(rts<<1))*tupSize] = (ru2*y1) + (ru4*y2) + (ru1*y3) + (ru3*y4) - shift; + y[(tensorOffset+rts*3)*tupSize] = (ru3*y1) + (ru1*y2) + (ru4*y3) + (ru2*y4) - shift; + } + } + } + else if(p == 7) { + hDim_t temp1 = rts*6; + ring ru1 = ruinv[rustride*tupSize]; + ring ru2 = ruinv[(rustride<<1)*tupSize]; + ring ru3 = ruinv[(rustride*3)*tupSize]; + ring ru4 = ruinv[(rustride<<2)*tupSize]; + ring ru5 = ruinv[(rustride*5)*tupSize]; + ring ru6 = ruinv[(rustride*6)*tupSize]; + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + ring y1, y2, y3, y4, y5, y6, shift; + y1 = y[tensorOffset*tupSize]; + y2 = y[(tensorOffset+rts)*tupSize]; + y3 = y[(tensorOffset+(rts<<1))*tupSize]; + y4 = y[(tensorOffset+3*rts)*tupSize]; + y5 = y[(tensorOffset+(rts<<2))*tupSize]; + y6 = y[(tensorOffset+rts*5)*tupSize]; + + shift = (ru6*y1) + (ru5*y2) + (ru4*y3) + (ru3*y4) + (ru2*y5) + (ru1*y6); + + y[tensorOffset*tupSize] += y2 + y3 + y4 + y5 + y6 - shift; + y[(tensorOffset+rts)*tupSize] = (ru1*y1) + (ru2*y2) + (ru3*y3) + (ru4*y4) + (ru5*y5) + (ru6*y6) - shift; + y[(tensorOffset+(rts<<1))*tupSize] = (ru2*y1) + (ru4*y2) + (ru6*y3) + (ru1*y4) + (ru3*y5) + (ru5*y6) - shift; + y[(tensorOffset+rts*3)*tupSize] = (ru3*y1) + (ru6*y2) + (ru2*y3) + (ru5*y4) + (ru1*y5) + (ru4*y6) - shift; + y[(tensorOffset+(rts<<2))*tupSize] = (ru4*y1) + (ru1*y2) + (ru5*y3) + (ru2*y4) + (ru6*y5) + (ru3*y6) - shift; + y[(tensorOffset+rts*5)*tupSize] = (ru5*y1) + (ru3*y2) + (ru1*y3) + (ru6*y4) + (ru4*y5) + (ru2*y6) - shift; + } + } + } + else { + ring* tempSpace = (ring*)malloc((p-1)*sizeof(ring)); + hDim_t temp1 = rts*(p-1); + for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(hDim_t modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + ring shift; + shift = 0; + for(hDim_t row = 0; row < p-1; row++) { + shift += (y[(tensorOffset+row*rts)*tupSize]*ruinv[(p-row-1)*rustride*tupSize]); + tempSpace[row] = 0; + for(hDim_t col = 0; col < p-1; col++) { + tempSpace[row] += (y[(tensorOffset+col*rts)*tupSize]*ruinv[((row*(col+1)) % p)*rustride*tupSize]); + } + } + + for(hDim_t row = 0; row < p-1; row++) { + y[(tensorOffset+rts*row)*tupSize] = tempSpace[row] - shift; + } + } + } + free(tempSpace); + } +} + +template void ppDFT (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, + PrimeExponent pe, hDim_t rustride, ring* ru, ring* temp) +{ + hDim_t p = pe.prime; + hShort_t e = pe.exponent; + + if(e == 0) { + return; + } + + hDim_t primeRuStride = rustride*ipow(p,e-1); + + hShort_t i; + + hDim_t ltsScale = ipow(p,e-1); + hDim_t rtsScale = 1; + hDim_t twidRuStride = rustride; + for(i = 0; i < e; i++) { + hDim_t rtsDim = rts*rtsScale; + dftp (y, tupSize, lts*ltsScale, rtsDim, p, primeRuStride, ru, temp); + dftTwiddle (y, tupSize, lts, rtsDim, pe, ltsScale*p, twidRuStride, ru); + + ltsScale /= p; + rtsScale *= p; + twidRuStride *= p; + pe.exponent -= 1; + } +} + +template void ppDFTInv (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, + PrimeExponent pe, hDim_t rustride, ring* ru, ring* temp) +{ + hDim_t p = pe.prime; + hShort_t e = pe.exponent; + + if(e == 0) { + return; + } + hDim_t primeRuStride = rustride*ipow(p,e-1); + + hShort_t i; + + hDim_t ltsScale = 1; + hDim_t rtsScale = ipow(p,e-1); + hDim_t twidRuStride = primeRuStride; + pe.exponent = 1; + for(i = 0; i < e; i++) { + hDim_t rtsDim = rts*rtsScale; + hDim_t ltsScaleP = ltsScale*p; + dftTwiddle (y, tupSize, lts, rtsDim, pe, ltsScaleP, twidRuStride, ru); + dftp (y, tupSize, lts*ltsScale, rtsDim, p, primeRuStride, ru, temp); + + ltsScale = ltsScaleP; + rtsScale /= p; + twidRuStride /= p; + pe.exponent += 1; + } +} + +template void ppcrt (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, + PrimeExponent pe, ring* ru) +{ + hDim_t p = pe.prime; + hDim_t e = pe.exponent; + hDim_t mprime = ipow(p,e-1); + ring* temp = 0; + if(p >= DFTP_GENERIC_SIZE) { + temp = (ring*)malloc(p*sizeof(ring)); + } + + crtp (y, tupSize, lts*mprime, rts, p, mprime, ru); + crtTwiddle (y, tupSize, lts, rts, pe, ru); + pe.exponent -= 1; + ppDFT (y, tupSize, lts, rts*(p-1), pe, p, ru, temp); + pe.exponent += 1; + + if(p >= DFTP_GENERIC_SIZE) { + free(temp); + } +} + +template void ppcrtinv (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, + PrimeExponent pe, ring* ru) +{ + hDim_t p = pe.prime; + hDim_t e = pe.exponent; + hDim_t mprime = ipow(p,e-1); + ring* temp = 0; + if(p >= DFTP_GENERIC_SIZE) { + temp = (ring*)malloc(p*sizeof(ring)); + } + + pe.exponent -= 1; + ppDFTInv (y, tupSize, lts, rts*(p-1), pe, p, ru, temp); + pe.exponent += 1; + crtTwiddle (y, tupSize, lts, rts, pe, ru); + crtpinv (y, tupSize, lts*mprime, rts, p, mprime, ru); + + if(p >= DFTP_GENERIC_SIZE) { + free(temp); + } +} + +extern "C" void tensorCRTRq (hShort_t tupSize, Zq* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, Zq** ru, hInt_t* qs) +{ + tensorFuserCRT (y, tupSize, ppcrt, totm, peArr, sizeOfPE, ru, qs); + canonicalizeZq(y,tupSize,totm,qs); +} + +//takes inverse rus +extern "C" void tensorCRTInvRq (hShort_t tupSize, Zq* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, + Zq** ruinv, Zq* mhatInv, hInt_t* qs) +{ + tensorFuserCRT (y, tupSize, ppcrtinv, totm, peArr, sizeOfPE, ruinv, qs); + for (hShort_t i = 0; i < tupSize; i++) { + Zq::q = qs[i]; + for (hDim_t j = 0; j < totm; j++) { + //careful here! I'm not setting the global q, so I can't rely on Zq multiplication + y[j*tupSize+i] = y[j*tupSize+i]*mhatInv[i]; + } + } + canonicalizeZq(y,tupSize,totm,qs); +} + +extern "C" void tensorCRTC (hShort_t tupSize, Complex* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, Complex** ru) +{ + tensorFuserCRT (y, tupSize, ppcrt, totm, peArr, sizeOfPE, ru, (hInt_t*)0); +} + +//takes inverse rus +extern "C" void tensorCRTInvC (hShort_t tupSize, Complex* y, hDim_t totm, PrimeExponent* peArr, + hShort_t sizeOfPE, Complex** ruinv, Complex* mhatInv) +{ + tensorFuserCRT (y, tupSize, ppcrtinv, totm, peArr, sizeOfPE, ruinv, (hInt_t*)0); + for (hShort_t i = 0; i < tupSize; i++) { + for (hDim_t j = 0; j < totm; j++) { + y[j*tupSize+i] *= mhatInv[i]; + } + } +} diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/g.cpp b/Crypto/Lol/Cyclotomic/Tensor/CPP/g.cpp new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/g.cpp @@ -0,0 +1,273 @@ +/* +Module : g.cpp +Description : Multiplication and division by 'g' in different bases. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX +*/ + +#include "types.h" +#include "tensor.h" +#include "common.h" + +template void gPow (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) +{ + if (p == 2) {return;} + hDim_t tmp1 = rts*(p-1); + hDim_t tmp2 = tmp1 - rts; + hDim_t blockOffset, modOffset; + hDim_t i; + for (blockOffset = 0; blockOffset < lts; ++blockOffset) { + hDim_t tmp3 = blockOffset * tmp1; + for (modOffset = 0; modOffset < rts; ++modOffset) { + hDim_t tensorOffset = tmp3 + modOffset; + ring last = y[(tensorOffset + tmp2)*tupSize]; + for (i = p-2; i != 0; --i) { + hDim_t idx = tensorOffset + i * rts; + y[idx*tupSize] += (last - y[(idx-rts)*tupSize]); + } + y[tensorOffset*tupSize] += last; + } + } +} + +template void gDec (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) +{ + if (p == 2) {return;} + hDim_t tmp1 = rts*(p-1); + hDim_t blockOffset; + hDim_t modOffset; + hDim_t i; + + for (blockOffset = 0; blockOffset < lts; ++blockOffset) { + hDim_t tmp2 = blockOffset * tmp1; + for (modOffset = 0; modOffset < rts; ++modOffset) { + hDim_t tensorOffset = tmp2 + modOffset; + ring acc = y[tensorOffset*tupSize]; + for (i = p-2; i != 0; --i) { + hDim_t idx = tensorOffset + i * rts; + acc += y[idx*tupSize]; + y[idx*tupSize] -= y[(idx-rts)*tupSize]; + } + y[tensorOffset*tupSize] += acc; + } + } +} + +template void gInvPow (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) +{ + if (p == 2) {return;} + hDim_t tmp1 = rts * (p-1); + hDim_t blockOffset, modOffset; + hDim_t i; + + for (blockOffset = 0; blockOffset < lts; ++blockOffset) { + hDim_t tmp2 = blockOffset * tmp1; + for (modOffset = 0; modOffset < rts; ++modOffset) { + hDim_t tensorOffset = tmp2 + modOffset; + ring lelts; + lelts = 0; + for (i = 0; i < p-1; ++i) { + lelts += y[(tensorOffset + i*rts)*tupSize]; + } + ring relts; + relts = 0; + for (i = p-2; i >= 0; --i) { + hDim_t idx = tensorOffset + i*rts; + ring z = y[idx*tupSize]; + ring lmul, rmul; + lmul = p-1-i; + rmul = i+1; + y[idx*tupSize] = lmul * lelts - rmul * relts; + lelts -= z; + relts += z; + } + } + } +} + +template void gInvDec (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) +{ + if (p == 2) {return;} + hDim_t blockOffset; + hDim_t modOffset; + hDim_t i; + hDim_t tmp1 = rts*(p-1); + + for (blockOffset = 0; blockOffset < lts; ++blockOffset) { + hDim_t tmp2 = blockOffset*tmp1; + for (modOffset = 0; modOffset < rts; ++modOffset) { + hDim_t tensorOffset = tmp2 + modOffset; + ring lastOut; + lastOut = 0; + for (i=1; i < p; ++i) { + ring ri; + ri = i; + lastOut += (ri * y[(tensorOffset + (i-1)*rts)*tupSize]); + } + ring rp; + rp = p; + ring acc = lastOut; + for (i = p-2; i > 0; --i) { + hDim_t idx = tensorOffset + i*rts; + ring tmp = acc; + acc -= y[idx*tupSize]*rp; + y[idx*tupSize] = tmp; + } + y[tensorOffset*tupSize] = acc; + } + } +} + +extern "C" void tensorGPowR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, gPow, totm, peArr, sizeOfPE, (hInt_t*)0); +} + +extern "C" void tensorGPowRq (hShort_t tupSize, Zq* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs) +{ + tensorFuserPrime (y, tupSize, gPow, totm, peArr, sizeOfPE, qs); + canonicalizeZq(y,tupSize,totm,qs); +} + +extern "C" void tensorGPowC (hShort_t tupSize, Complex* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, gPow, totm, peArr, sizeOfPE, (hInt_t*)0); +} + +extern "C" void tensorGDecR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, gDec, totm, peArr, sizeOfPE, (hInt_t*)0); +} + +extern "C" void tensorGDecRq (hShort_t tupSize, Zq* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs) +{ + tensorFuserPrime (y, tupSize, gDec, totm, peArr, sizeOfPE, qs); + canonicalizeZq(y,tupSize,totm,qs); +} + +extern "C" void tensorGDecC (hShort_t tupSize, Complex* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, gDec, totm, peArr, sizeOfPE, (hInt_t*)0); +} + +hInt_t oddRad(PrimeExponent* peArr, hShort_t sizeOfPE) { + hInt_t oddrad; + oddrad = 1; + for(int i = 0; i < sizeOfPE; i++) { + hShort_t p = peArr[i].prime; + if (p != 2) { + oddrad *= peArr[i].prime; + } + } + return oddrad; +} + +extern "C" hShort_t tensorGInvPowR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, gInvPow, totm, peArr, sizeOfPE, (hInt_t*)0); + + hInt_t oddrad = oddRad(peArr, sizeOfPE); + + for(int i = 0; i < tupSize*totm; i++) { + if (y[i] % oddrad) { + y[i] /= oddrad; + } + else { + return 0; + } + } + return 1; +} + +extern "C" hShort_t tensorGInvPowRq (hShort_t tupSize, Zq* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs) +{ + tensorFuserPrime (y, tupSize, gInvPow, totm, peArr, sizeOfPE, qs); + + hInt_t oddrad = oddRad(peArr, sizeOfPE); + + for(int i = 0; i < tupSize; i++) { + Zq::q = qs[i]; // global update + hInt_t ori = reciprocal(Zq::q, oddrad); + Zq oddradInv; + oddradInv = ori; + if (ori == 0) { + return 0; // error condition + } + for(hDim_t j = 0; j < totm; j++) { + y[j*tupSize+i] *= oddradInv; + } + } + + canonicalizeZq(y,tupSize,totm,qs); + return 1; +} + +extern "C" hShort_t tensorGInvPowC (hShort_t tupSize, Complex* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, gInvPow, totm, peArr, sizeOfPE, (hInt_t*)0); + + hInt_t oddrad = oddRad(peArr, sizeOfPE); + Complex oddradInv; + oddradInv = 1 / oddrad; + for(int i = 0; i < tupSize*totm; i++) { + y[i] *= oddradInv; + } + return 1; +} + +extern "C" hShort_t tensorGInvDecR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, gInvDec, totm, peArr, sizeOfPE, (hInt_t*)0); + + hInt_t oddrad = oddRad(peArr, sizeOfPE); + + for(int i = 0; i < tupSize*totm; i++) { + if (y[i] % oddrad) { + y[i] /= oddrad; + } + else { + return 0; + } + } + return 1; +} + +extern "C" hShort_t tensorGInvDecRq (hShort_t tupSize, Zq* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs) +{ + tensorFuserPrime (y, tupSize, gInvDec, totm, peArr, sizeOfPE, qs); + + hInt_t oddrad = oddRad(peArr, sizeOfPE); + + for(int i = 0; i < tupSize; i++) { + Zq::q = qs[i]; // global update + hInt_t ori = reciprocal(Zq::q, oddrad); + Zq oddradInv; + oddradInv = ori; + if (ori == 0) { + return 0; // error condition + } + for(hDim_t j = 0; j < totm; j++) { + y[j*tupSize+i] *= oddradInv; + } + } + + canonicalizeZq(y,tupSize,totm,qs); + return 1; +} + +extern "C" hShort_t tensorGInvDecC (hShort_t tupSize, Complex* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, gInvDec, totm, peArr, sizeOfPE, (hInt_t*)0); + + hInt_t oddrad = oddRad(peArr, sizeOfPE); + Complex oddradInv; + oddradInv = 1 / oddrad; + for(int i = 0; i < tupSize*totm; i++) { + y[i] *= oddradInv; + } + return 1; +} diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/l.cpp b/Crypto/Lol/Cyclotomic/Tensor/CPP/l.cpp new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/l.cpp @@ -0,0 +1,98 @@ +/* +Module : l.cpp +Description : Powerful <-> Decoding basis conversion. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX +*/ + +#include "types.h" +#include "tensor.h" + +template void lp (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) +{ + hDim_t blockOffset; + hDim_t modOffset; + int i; + + if(p == 2) {return;} + + hDim_t tmp1 = rts*(p-1); + for (blockOffset = 0; blockOffset < lts; ++blockOffset) { + hDim_t tmp2 = blockOffset*tmp1; + for (modOffset = 0; modOffset < rts; ++modOffset) { + hDim_t idx = tmp2 + modOffset + rts; + for (i = 1; i < p-1; ++i) { + y[idx*tupSize] += y[(idx-rts)*tupSize]; + idx += rts; + } + } + } +} + +template void lpInv (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) +{ + hDim_t blockOffset; + hDim_t modOffset; + int i; + + if(p == 2) {return;} + + hDim_t tmp1 = rts*(p-1); + for (blockOffset = 0; blockOffset < lts; ++blockOffset) { + hDim_t tmp2 = blockOffset*tmp1; + for (modOffset = 0; modOffset < rts; ++ modOffset) { + hDim_t tensorOffset = tmp2 + modOffset; + hDim_t idx = tensorOffset + (p-2) * rts; + for (i = p-2; i != 0; --i) { + y[idx*tupSize] -= y[(idx-rts)*tupSize] ; + idx -= rts; + } + } + } +} + +extern "C" void tensorLRq (hShort_t tupSize, Zq* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs) +{ + tensorFuserPrime (y, tupSize, lp, totm, peArr, sizeOfPE, qs); + canonicalizeZq(y,tupSize,totm,qs); +} + +extern "C" void tensorLR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, lp, totm, peArr, sizeOfPE, (hInt_t*)0); +} + +extern "C" void tensorLDouble (hShort_t tupSize, double* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, lp, totm, peArr, sizeOfPE, (hInt_t*)0); +} + +extern "C" void tensorLC (hShort_t tupSize, Complex* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, lp, totm, peArr, sizeOfPE, (hInt_t*)0); +} + +extern "C" void tensorLInvRq (hShort_t tupSize, Zq* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs) +{ + tensorFuserPrime (y, tupSize, lpInv, totm, peArr, sizeOfPE, qs); + canonicalizeZq(y,tupSize,totm,qs); +} + +extern "C" void tensorLInvR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, lpInv, totm, peArr, sizeOfPE, (hInt_t*)0); +} + +extern "C" void tensorLInvDouble (hShort_t tupSize, double* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, lpInv, totm, peArr, sizeOfPE, (hInt_t*)0); +} + +extern "C" void tensorLInvC (hShort_t tupSize, Complex* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + tensorFuserPrime (y, tupSize, lpInv, totm, peArr, sizeOfPE, (hInt_t*)0); +} diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/mul.cpp b/Crypto/Lol/Cyclotomic/Tensor/CPP/mul.cpp new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/mul.cpp @@ -0,0 +1,35 @@ +/* +Module : mul.cpp +Description : zipWith (*). +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX +*/ + +#include "types.h" + +template void zipWithStar (ring* a, ring* b, hShort_t tupSize, hDim_t totm, hInt_t* qs) +{ + for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) { + if(qs) { + Zq::q = qs[tupIdx]; + } + for(int i = 0; i < totm; i++) { + a[i*tupSize+tupIdx] *= b[i*tupSize+tupIdx]; + } + } +} + +//a = zipWith (*) a b +extern "C" void mulRq (hShort_t tupSize, Zq* a, Zq* b, hDim_t totm, hInt_t* qs) +{ + zipWithStar(a, b, tupSize, totm, qs); +} + +extern "C" void mulC (hShort_t tupSize, Complex* a, Complex* b, hDim_t totm) +{ + zipWithStar(a, b, tupSize, totm, (hInt_t*)0); +} diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/norm.cpp b/Crypto/Lol/Cyclotomic/Tensor/CPP/norm.cpp new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/norm.cpp @@ -0,0 +1,80 @@ +/* +Module : norm.cpp +Description : Compute g*norm(x)^2. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX +*/ + +#include "types.h" +#include "tensor.h" + +template void pNormSq (ring* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) +{ + hDim_t blockOffset; + hDim_t modOffset; + hDim_t i; + + if(p==2) {return;} + + hDim_t tmp1 = rts*(p-1); + for (blockOffset = 0; blockOffset < lts; ++blockOffset) { + hDim_t tmp2 = blockOffset*tmp1; + for (modOffset = 0; modOffset < rts; ++modOffset) { + hDim_t tensorOffset = tmp2 + modOffset; + ring sum = 0; + for (i = 0; i < p-1; ++i) { + sum += y[(tensorOffset + i*rts)*tupSize]; + } + for (i = 0; i < p-1; ++i) { + y[(tensorOffset + i*rts)*tupSize] += sum; + } + } + } +} + +extern "C" void tensorNormSqR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + hInt_t* tempSpace = (hInt_t*)malloc(totm*tupSize*sizeof(hInt_t)); + for(hDim_t i = 0; i < totm*tupSize; i++) { + tempSpace[i]=y[i]; + } + + tensorFuserPrime(y, tupSize, pNormSq, totm, peArr, sizeOfPE, (hInt_t*)0); + + //do dot product and return in index 0 + for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) { + hInt_t dotprod = 0; + for(hDim_t i = 0; i < totm; i++) { + dotprod += (tempSpace[i*tupSize+tupIdx]*y[i*tupSize+tupIdx]); + } + + y[tupIdx] = dotprod; + } + + free(tempSpace); +} + +extern "C" void tensorNormSqD (hShort_t tupSize, double* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) +{ + double* tempSpace = (double*)malloc(totm*tupSize*sizeof(double)); + for(hDim_t i = 0; i < totm*tupSize; i++) { + tempSpace[i]=y[i]; + } + tensorFuserPrime(y, tupSize, pNormSq, totm, peArr, sizeOfPE, (hInt_t*)0); + + //do dot product and return in index 0 + for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) { + double dotprod = 0; + for(hDim_t i = 0; i < totm; i++) { + dotprod += (tempSpace[i*tupSize+tupIdx]*y[i*tupSize+tupIdx]); + } + + y[tupIdx] = dotprod; + } + + free(tempSpace); +} diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/random.cpp b/Crypto/Lol/Cyclotomic/Tensor/CPP/random.cpp new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/random.cpp @@ -0,0 +1,64 @@ +/* +Module : random.cpp +Description : Gaussian sampling on a lattice. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX +*/ + +#include "types.h" +#include "tensor.h" +#include "common.h" +#include + +// I had been negating the ru-idx, but this was causing a *negative* mod, resulting in a hard-to-find bug +// current behavior (taking rus, rather than ruInv) matches RT +void primeD (double *y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hDim_t rustride, Complex* ru) +{ + if(p == 2) { + return; + } + hDim_t blockOffset, modOffset, tensorOffset; + double *tempSpace = (double*)malloc((p-1)*sizeof(double)); + hDim_t temp1 = rts*(p-1); + for(blockOffset = 0; blockOffset < lts; blockOffset++) { + hDim_t temp2 = blockOffset*temp1; + for(modOffset = 0; modOffset < rts; modOffset++) { + tensorOffset = temp2 + modOffset; + hDim_t row, col; + + for(row = 0; row < p-1; row++) { + double acc = 0; + for(col = 1; col <= (p>>1); col++) { + acc += 2 * ru[((row*col) % p)*rustride*tupSize].real * y[(tensorOffset+rts*(col-1))*tupSize]; + } + for(col = (p>>1)+1; col <= p-1; col++) { + acc += 2 * ru[((row*col) % p)*rustride*tupSize].imag * y[(tensorOffset+rts*(col-1))*tupSize]; + } + tempSpace[row] = acc/sqrt(2); + } + + for(row = 0; row < p-1; row++) { + y[(tensorOffset+rts*row)*tupSize] = tempSpace[row]; + } + } + } + free(tempSpace); +} + +void ppD (double *y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, Complex *ru) +{ + hDim_t p = pe.prime; + hDim_t e = pe.exponent; + hDim_t mprime = ipow(p,e-1); + primeD (y, tupSize, lts*mprime, rts, p, mprime, ru); +} + +//the contents of y will be destroyed, but should be initialized in Haskell-land to independent Guassians over the reals +extern "C" void tensorGaussianDec (hShort_t tupSize, double *y, hDim_t totm, PrimeExponent *peArr, hShort_t sizeOfPE, Complex** ru) +{ + tensorFuserCRT (y, tupSize, ppD, totm, peArr, sizeOfPE, ru, (hInt_t*)0); +} diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/tensor.h b/Crypto/Lol/Cyclotomic/Tensor/CPP/tensor.h new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/tensor.h @@ -0,0 +1,67 @@ +/* +Module : tensor.h +Description : Templates for the tensor DSL. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX +*/ + +#ifndef TENSOR_CPP_ +#define TENSOR_CPP_ + +#include "types.h" +#include "common.h" +#ifdef __cplusplus +template +using primeFunc = void (*) (ring*, hShort_t, hDim_t, hDim_t, hDim_t); + +template +using primeCRTFunc = void (*) (ringy*, hShort_t, hDim_t, hDim_t, PrimeExponent, ringru*); + +//for square transforms +template void tensorFuserPrime (ring* y, hShort_t tupSize, primeFunc f, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs) +{ + hDim_t lts = totm; + hDim_t rts = 1; + hShort_t i; + + for (i = 0; i < sizeOfPE; ++i) { + PrimeExponent pe = peArr[i]; + hDim_t ipow_pe = ipow(pe.prime, (pe.exponent-1)); + hDim_t dim = (pe.prime-1) * ipow_pe; // the totient of pe + lts /= dim; + for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) { + if(qs) { + Zq::q = qs[tupIdx]; // global update + } + (*f) (y+tupIdx, tupSize, lts*ipow_pe, rts, pe.prime); + } + rts *= dim; + } +} + +template void tensorFuserCRT (ringy* y, hShort_t tupSize, primeCRTFunc f, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, ringru** ru, hInt_t* qs) +{ + hDim_t lts = totm; + hDim_t rts = 1; + hShort_t i; + + for (i = 0; i < sizeOfPE; ++i) { + PrimeExponent pe = peArr[i]; + hDim_t ipow_pe = ipow(pe.prime, (pe.exponent-1)); + hDim_t dim = (pe.prime-1) * ipow_pe; // the totient of pe + lts /= dim; + for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) { + if(qs) { + Zq::q = qs[tupIdx]; // global update + } + (*f) (y+tupIdx, tupSize, lts, rts, pe, ru[i]+tupIdx); + } + rts *= dim; + } +} +#endif /* __cplusplus */ +#endif /* TENSOR_CPP_ */ diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/types.h b/Crypto/Lol/Cyclotomic/Tensor/CPP/types.h new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/types.h @@ -0,0 +1,170 @@ +/* +Module : types.h +Description : C++ types and function declarations. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX +*/ + +#ifndef TENSORTYPES_H_ +#define TENSORTYPES_H_ + +#include +#include +#include + +typedef int64_t hInt_t ; +typedef int32_t hDim_t ; +typedef int16_t hShort_t ; +typedef int8_t hByte_t ; + +typedef struct +{ + hShort_t prime; + hShort_t exponent; +} PrimeExponent; + + +hInt_t reciprocal (hInt_t a, hInt_t b); + +#define ASSERT(EXP) { \ + if (!(EXP)) { \ + fprintf (stderr, "Assertion in file '%s' line %d : " #EXP " is false\n", __FILE__, __LINE__); \ + exit(-1); \ + } \ +} + +//http://stackoverflow.com/questions/37572628 +#ifdef __cplusplus +//http://stackoverflow.com/a/4421719 +class Zq +{ +public: + hInt_t x; + + static hInt_t q; // declared here, defined in generalfuncs.cpp + + Zq& operator=(const hInt_t& c) + { + this->x = c % q; + return *this; + } + Zq& operator+=(const Zq& b) + { + this->x += b.x; + this->x %= q; + return *this; + } + Zq& operator-=(const Zq& b) + { + this->x -= b.x; + this->x %= q; + return *this; + } + Zq& operator*=(const Zq& b) + { + this->x *= b.x; + this->x %= q; + return *this; + } + Zq& operator/=(const Zq& b) + { + Zq binv; + binv = reciprocal(q,b.x); + ASSERT (binv.x); // binv == 0 indicates that x is not invertible mod q + *this *= binv; + return *this; + } +}; +inline Zq operator+(Zq a, const Zq& b) +{ + a += b; + return a; +} +inline Zq operator-(Zq a, const Zq& b) +{ + a -= b; + return a; +} +inline Zq operator*(Zq a, const Zq& b) +{ + a *= b; + return a; +} +inline Zq operator/(Zq a, const Zq& b) +{ + a /= b; + return a; +} + +void canonicalizeZq (Zq* y, hShort_t tupSize, hDim_t totm, hInt_t* qs); + +class Complex +{ +public: + double real; + double imag; + + Complex& operator=(const hInt_t& c) + { + this->real = c; + this->imag = 0; + return *this; + } + Complex& operator+=(const Complex& b) + { + this->real = this->real+b.real; + this->imag = this->imag+b.imag; + return *this; + } + Complex& operator-=(const Complex& b) + { + this->real = this->real-b.real; + this->imag = this->imag-b.imag; + return *this; + } + Complex& operator*=(const Complex& b) + { + double a = this->real; + this->real = (a*b.real)-(this->imag*b.imag); + this->imag = (a*b.imag)+(this->imag*b.real); + return *this; + } + Complex& operator/=(const Complex& b) + { + Complex bconj; + bconj.real = b.real; + bconj.imag = -b.imag; + *this *= bconj; + double den = (b.real*b.real+b.imag*b.imag); + this->real /= den; + this->imag /= den; + return *this; + } +}; +inline Complex operator+(Complex a, const Complex& b) +{ + a += b; + return a; +} +inline Complex operator-(Complex a, const Complex& b) +{ + a -= b; + return a; +} +inline Complex operator*(Complex a, const Complex& b) +{ + a *= b; + return a; +} +inline Complex operator/(Complex a, const Complex& b) +{ + a /= b; + return a; +} + +#endif /* __cplusplus */ +#endif /* TENSORTYPES_H_ */ diff --git a/Crypto/Lol/Cyclotomic/Tensor/CPP/zq.cpp b/Crypto/Lol/Cyclotomic/Tensor/CPP/zq.cpp new file mode 100644 --- /dev/null +++ b/Crypto/Lol/Cyclotomic/Tensor/CPP/zq.cpp @@ -0,0 +1,50 @@ +/* +Module : zq.cpp +Description : Implementation of Z_q-specific functions. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX +*/ + +#include "types.h" +#include "common.h" + +// a is the field size. we are looking for reciprocal of b +hInt_t reciprocal (hInt_t a, hInt_t b) +{ + hInt_t fieldSize = a; + + hInt_t y = 1; + hInt_t lasty = 0; + while (b != 0) { + hInt_t quotient = a / b; + hInt_t tmp = a % b; + a = b; + b = tmp; + tmp = y; + y = lasty - quotient*y; + lasty = tmp; + } + // if a!=1, then b is not invertible mod a + if(a!=1) { + return 0; + } + + // this actually returns EITHER the reciprocal OR reciprocal + fieldSize + hInt_t res = lasty + fieldSize; + return res; +} + +void canonicalizeZq (Zq* y, hShort_t tupSize, hDim_t totm, hInt_t* qs) { + for(int tupIdx = 0; tupIdx + Copyright (C) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -0,0 +1,5 @@ +This package contains a fast C++ implementation of the 'Tensor' interface for +Lol, via Haskell's FFI. + +You can test this package by running `stack test lol-cpp`, and benchmark by +running `stack bench lol-cpp`. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/benchmarks/BenchCPPMain.hs b/benchmarks/BenchCPPMain.hs new file mode 100644 --- /dev/null +++ b/benchmarks/BenchCPPMain.hs @@ -0,0 +1,83 @@ +{-| +Module : BenchCPPMain +Description : Main driver for CPP benchmarks. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX + +Main driver for CPP benchmarks. +-} + +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE NamedFieldPuns #-} +{-# LANGUAGE TypeOperators #-} + +module BenchCPPMain where + +import Crypto.Lol.Benchmarks +import Crypto.Lol.Benchmarks.Standard +import Crypto.Lol.Cyclotomic.Tensor.CPP +import Crypto.Lol.Factored +import qualified Crypto.Lol.Utils.PrettyPrint.Diagnostic as D +import qualified Crypto.Lol.Utils.PrettyPrint.Table as T +import Crypto.Random.DRBG + +import Data.Proxy + +-- choose which layers of Lol to benchmark +ls :: [String] +ls = [ + "STensor", + "Tensor", + "SUCyc", + "UCyc", + "Cyc" + ] + +-- choose which operations to benchmark +bs :: [String] +bs = [ + "unzipPow", + "unzipDec", + "unzipCRT", + "zipWith (*)", + "crt", + "crtInv", + "l", + "lInv", + "*g Pow", + "*g Dec", + "*g CRT", + "divg Pow", + "divg Dec", + "divg CRT", + "lift", + "error", + "twacePow", + "twaceDec", + "twaceCRT", + "embedPow", + "embedDec", + "embedCRT" + ] + +main :: IO () +main = diagnosticMain + +tableMain :: IO () +tableMain = do + let opts = (T.defaultOpts $ Just "UCyc"){T.benches=bs} + g1 <- defaultBenches (Proxy::Proxy CT) + mapM_ (T.prettyBenches opts) g1 + +diagnosticMain :: IO () +diagnosticMain = do + let opts = D.defaultOpts{D.levels=ls, D.benches=bs} + b1 <- benchGroup "Single Index" + [oneIdxBenches (Proxy::Proxy '(F64*F9*F25, Zq 14401)) (Proxy::Proxy CT) (Proxy::Proxy HashDRBG)] + b2 <- benchGroup "Twace-Embed" + [twoIdxBenches (Proxy::Proxy '(F64*F9*F25, F64*F9*F25, Zq 14401)) (Proxy::Proxy CT)] + mapM_ (D.prettyBenches opts) [b1,b2] diff --git a/lol-cpp.cabal b/lol-cpp.cabal new file mode 100644 --- /dev/null +++ b/lol-cpp.cabal @@ -0,0 +1,135 @@ +name: lol-cpp +-- The package version. See the Haskell package versioning policy (PVP) +-- for standards guiding when and how versions should be incremented. +-- http://www.haskell.org/haskellwiki/Package_versioning_policy +-- PVP summary: +-+------- breaking API changes +-- | | +----- non-breaking API additions +-- | | | +--- code changes with no API change +version: 0.0.0.1 +synopsis: A fast C++ backend for . +homepage: https://github.com/cpeikert/Lol +Bug-Reports: https://github.com/cpeikert/Lol/issues +license: GPL-2 +license-file: LICENSE +author: Eric Crockett , Chris Peikert +maintainer: Eric Crockett +copyright: Eric Crockett, Chris Peikert +category: Crypto +stability: experimental +build-type: Simple +extra-source-files: README, CHANGES.md, + Crypto/Lol/Cyclotomic/Tensor/CPP/*.h, + Crypto/Lol/Cyclotomic/Tensor/CPP/*.cpp +cabal-version: >= 1.10 +description: + Λ ∘ λ (Lol) is a general-purpose library for ring-based lattice cryptography. + This package provides a C++ implementation of Lol's Tensor interface. +source-repository head + type: git + location: https://github.com/cpeikert/Lol + +-- For information on compiling C with cabal: http://blog.ezyang.com/2010/06/setting-up-cabal-the-ffi-and-c2hs/ + +Flag llvm + Description: Compile via LLVM. This produces much better object code, + but you need to have the LLVM compiler installed. + -- If you enable this and get errors like "Error: can't resolve `.rodata' {.rodata section}" + -- then GHC doesn't like your version of LLVM! + Default: False + +Flag opt + Description: Turn on library optimizations + Default: True + +library + default-language: Haskell2010 + ghc-options: -fwarn-dodgy-imports + cc-options: -std=c++11 + Include-dirs: Crypto/Lol/Cyclotomic/Tensor/CPP + -- Due to #12152, the file containing the definition of `Zq::q` must be linked first, + -- otherwise dynamic linking (`cabal repl` or `stack ghci`) results in the error: + -- "Loading temp shared object failed: /tmp/ghc54651_0/libghc_1.so: undefined symbol _ZN2Zq1qE" + -- For `cabal repl`, we can simply reorder the list so that the file that should be linked + -- first comes first in the list. However `stack ghci` always links alphabetically, + -- so we really just have to define `Zq::q` in the first file alphabetically. + C-sources: Crypto/Lol/Cyclotomic/Tensor/CPP/common.cpp, + Crypto/Lol/Cyclotomic/Tensor/CPP/crt.cpp, + Crypto/Lol/Cyclotomic/Tensor/CPP/g.cpp, + Crypto/Lol/Cyclotomic/Tensor/CPP/l.cpp, + Crypto/Lol/Cyclotomic/Tensor/CPP/mul.cpp, + Crypto/Lol/Cyclotomic/Tensor/CPP/norm.cpp, + Crypto/Lol/Cyclotomic/Tensor/CPP/random.cpp, + Crypto/Lol/Cyclotomic/Tensor/CPP/zq.cpp + + if flag(llvm) + ghc-options: -fllvm -optlo-O3 + + -- ghc optimizations + if flag(opt) + -- makes lift much faster! + ghc-options: -funfolding-use-threshold1000 + exposed-modules: + Crypto.Lol.Cyclotomic.Tensor.CPP + + other-modules: + Crypto.Lol.Cyclotomic.Tensor.CPP.Backend + Crypto.Lol.Cyclotomic.Tensor.CPP.Extension + Crypto.Lol.Cyclotomic.Tensor.CPP.Instances + + build-depends: + arithmoi >= 0.4.1.3, + base >= 4.9 && < 5, + bytestring, + constraints, + containers >= 0.5.6.2, + crypto-api, + data-default >= 0.3.0, + deepseq >= 1.4.1.1, + lol >= 0.6.0.0, + monadcryptorandom, + MonadRandom >= 0.2, + mtl >= 2.2.1, + numeric-prelude >= 0.4.2, + protocol-buffers, + protocol-buffers-descriptor, + random >= 1.1, + reflection >= 1.5.1, + repa>=3.4, + singletons >= 1.1.2.1, + th-desugar >= 1.5.4, + tagged-transformer >= 0.7, + template-haskell >= 2.2.0.0, + transformers >= 0.4.2.0, + vector>=0.11, + vector-th-unbox >= 0.2.1.0 + + other-extensions: TemplateHaskell + +Benchmark bench-lol-cpp + type: exitcode-stdio-1.0 + default-language: Haskell2010 + main-is: BenchCPPMain.hs + ghc-options: -main-is BenchCPPMain + hs-source-dirs: benchmarks + + ghc-options: -O2 -funfolding-creation-threshold=15000 -funfolding-use-threshold=1000 + + build-depends: + base >= 4.9 && < 5, + DRBG, + lol >= 0.6.0.0, + lol-benches, + lol-cpp + +test-suite test-lol-cpp + type: exitcode-stdio-1.0 + default-language: Haskell2010 + main-is: TestCPPMain.hs + ghc-options: -main-is TestCPPMain + hs-source-dirs: tests + ghc-options: -threaded -O2 + + build-depends: + base >= 4.9 && < 5, + lol-cpp, + lol-tests diff --git a/tests/TestCPPMain.hs b/tests/TestCPPMain.hs new file mode 100644 --- /dev/null +++ b/tests/TestCPPMain.hs @@ -0,0 +1,21 @@ +{-| +Module : TestCPPMain +Description : Main driver for CPP tests. +Copyright : (c) Eric Crockett, 2011-2017 + Chris Peikert, 2011-2017 +License : GPL-2 +Maintainer : ecrockett0@email.com +Stability : experimental +Portability : POSIX + +Main driver for CPP tests. +-} + +module TestCPPMain where + +import Crypto.Lol.Cyclotomic.Tensor.CPP +import Crypto.Lol.Tests.Standard +import Data.Proxy + +main :: IO () +main = defaultTestMain (Proxy::Proxy CT)