lol 0.1.0.0 → 0.2.0.0

raw patch · 101 files changed

+10904/−10105 lines, 101 filesdep +DRBGdep +binarydep +bytestringdep −storable-tupledep −timedep −type-naturaldep ~MonadRandomdep ~basedep ~constraintsPVP ok

version bump matches the API change (PVP)

Dependencies added: DRBG, binary, bytestring, criterion, crypto-api, template-haskell

Dependencies removed: storable-tuple, time, type-natural

Dependency ranges changed: MonadRandom, base, constraints, deepseq, mtl, random, repa, singletons, transformers, vector

API changes (from Hackage documentation)

Files

• CHANGES.md +9/−0
• Crypto/Lol.hs +21/−0
• Crypto/Lol/CRTrans.hs +138/−0
• Crypto/Lol/Cyclotomic/Cyc.hs +583/−0
• Crypto/Lol/Cyclotomic/Linear.hs +90/−0
• Crypto/Lol/Cyclotomic/RescaleCyc.hs +23/−0
• Crypto/Lol/Cyclotomic/Tensor.hs +447/−0
• Crypto/Lol/Cyclotomic/Tensor/CTensor.hs +474/−0
• Crypto/Lol/Cyclotomic/Tensor/CTensor/Backend.hs +304/−0
• Crypto/Lol/Cyclotomic/Tensor/CTensor/Extension.hs +121/−0
• Crypto/Lol/Cyclotomic/Tensor/CTensor/basic.c +175/−0
• Crypto/Lol/Cyclotomic/Tensor/CTensor/crt.c +1386/−0
• Crypto/Lol/Cyclotomic/Tensor/CTensor/g.c +828/−0
• Crypto/Lol/Cyclotomic/Tensor/CTensor/generalfuncs.c +305/−0
• Crypto/Lol/Cyclotomic/Tensor/CTensor/l.c +434/−0
• Crypto/Lol/Cyclotomic/Tensor/CTensor/norm.c +86/−0
• Crypto/Lol/Cyclotomic/Tensor/CTensor/random.c +74/−0
• Crypto/Lol/Cyclotomic/Tensor/CTensor/tensorTypes.h +231/−0
• Crypto/Lol/Cyclotomic/Tensor/RepaTensor.hs +258/−0
• Crypto/Lol/Cyclotomic/Tensor/RepaTensor/CRT.hs +215/−0
• Crypto/Lol/Cyclotomic/Tensor/RepaTensor/Dec.hs +74/−0
• Crypto/Lol/Cyclotomic/Tensor/RepaTensor/Extension.hs +193/−0
• Crypto/Lol/Cyclotomic/Tensor/RepaTensor/GL.hs +134/−0
• Crypto/Lol/Cyclotomic/Tensor/RepaTensor/RTCommon.hs +277/−0
• Crypto/Lol/Cyclotomic/UCyc.hs +591/−0
• Crypto/Lol/Factored.hs +27/−0
• Crypto/Lol/FactoredDefs.hs +429/−0
• Crypto/Lol/Gadget.hs +119/−0
• Crypto/Lol/GaussRandom.hs +90/−0
• Crypto/Lol/LatticePrelude.hs +276/−0
• Crypto/Lol/PosBin.hs +18/−0
• Crypto/Lol/PosBinDefs.hs +147/−0
• Crypto/Lol/Reflects.hs +57/−0
• Crypto/Lol/Types/Complex.hs +89/−0
• Crypto/Lol/Types/FiniteField.hs +174/−0
• Crypto/Lol/Types/IZipVector.hs +61/−0
• Crypto/Lol/Types/IrreducibleChar2.hs +67/−0
• Crypto/Lol/Types/Numeric.hs +210/−0
• Crypto/Lol/Types/Random.hs +44/−0
• Crypto/Lol/Types/ZPP.hs +23/−0
• Crypto/Lol/Types/ZmStar.hs +69/−0
• Crypto/Lol/Types/ZqBasic.hs +345/−0
• README +0/−4
• benchmarks/CycBenches.hs +115/−0
• benchmarks/Main.hs +12/−0
• benchmarks/ZqBenches.hs +39/−0
• lol.cabal +104/−60
• src/Crypto/Lol.hs +0/−21
• src/Crypto/Lol/Applications/SymmSHE.hs +0/−498
• src/Crypto/Lol/CRTrans.hs +0/−181
• src/Crypto/Lol/Cyclotomic/Cyc.hs +0/−184
• src/Crypto/Lol/Cyclotomic/Linear.hs +0/−89
• src/Crypto/Lol/Cyclotomic/Tensor.hs +0/−412
• src/Crypto/Lol/Cyclotomic/Tensor/CTensor.hs +0/−731
• src/Crypto/Lol/Cyclotomic/Tensor/CTensor/Extension.hs +0/−122
• src/Crypto/Lol/Cyclotomic/Tensor/CTensor/basic.c +0/−162
• src/Crypto/Lol/Cyclotomic/Tensor/CTensor/crt.c +0/−1305
• src/Crypto/Lol/Cyclotomic/Tensor/CTensor/g.c +0/−685
• src/Crypto/Lol/Cyclotomic/Tensor/CTensor/generalfuncs.c +0/−293
• src/Crypto/Lol/Cyclotomic/Tensor/CTensor/l.c +0/−410
• src/Crypto/Lol/Cyclotomic/Tensor/CTensor/random.c +0/−72
• src/Crypto/Lol/Cyclotomic/Tensor/CTensor/tensorTypes.h +0/−224
• src/Crypto/Lol/Cyclotomic/Tensor/RepaTensor.hs +0/−199
• src/Crypto/Lol/Cyclotomic/Tensor/RepaTensor/CRT.hs +0/−185
• src/Crypto/Lol/Cyclotomic/Tensor/RepaTensor/Extension.hs +0/−193
• src/Crypto/Lol/Cyclotomic/Tensor/RepaTensor/GL.hs +0/−112
• src/Crypto/Lol/Cyclotomic/Tensor/RepaTensor/Gauss.hs +0/−51
• src/Crypto/Lol/Cyclotomic/Tensor/RepaTensor/RTCommon.hs +0/−233
• src/Crypto/Lol/Cyclotomic/UCyc.hs +0/−700
• src/Crypto/Lol/Cyclotomic/Utility.hs +0/−22
• src/Crypto/Lol/Factored.hs +0/−481
• src/Crypto/Lol/Gadget.hs +0/−93
• src/Crypto/Lol/GaussRandom.hs +0/−90
• src/Crypto/Lol/LatticePrelude.hs +0/−248
• src/Crypto/Lol/Reflects.hs +0/−39
• src/Crypto/Lol/Types/Complex.hs +0/−89
• src/Crypto/Lol/Types/FiniteField.hs +0/−122
• src/Crypto/Lol/Types/IZipVector.hs +0/−58
• src/Crypto/Lol/Types/IrreducibleChar2.hs +0/−60
• src/Crypto/Lol/Types/Numeric.hs +0/−213
• src/Crypto/Lol/Types/PrimeField.hs +0/−36
• src/Crypto/Lol/Types/ZPP.hs +0/−23
• src/Crypto/Lol/Types/ZmStar.hs +0/−69
• src/Crypto/Lol/Types/ZqBasic.hs +0/−256
• test-suite/CycTests.hs +0/−105
• test-suite/Main.hs +0/−19
• test-suite/SHETests.hs +0/−442
• test-suite/TensorTests.hs +0/−368
• test-suite/TestTypes.hs +0/−90
• test-suite/ZqTests.hs +0/−56
• tests/CycTests.hs +92/−0
• tests/Main.hs +14/−0
• tests/TensorTests.hs +343/−0
• tests/ZqTests.hs +56/−0
• utils/Apply.hs +17/−0
• utils/Benchmarks.hs +46/−0
• utils/Gen.hs +27/−0
• utils/Harness/Cyc.hs +102/−0
• utils/TestTypes.hs +39/−0
• utils/Tests.hs +49/−0
• utils/Utils.hs +133/−0

@@ -1,6 +1,15 @@ Changelog for lol project ================================ +0.2.0.0+-----+ * Added benchmarks+ * Better performance for C backend with RNS base ring+ * UCyc exposes bases as type for safety+ * Other safety improvements throughout+ * Easier index representation with TemplateHaskell+ * Split SymmSHE into new package lol-apps.+ 0.1.0.0 -----  * Fixed bug in Box-Muller sampling routine.

@@ -0,0 +1,21 @@++-- | Re-exports primary interfaces.++module Crypto.Lol +( module Crypto.Lol.Cyclotomic.Cyc+, module Crypto.Lol.Gadget+, module Crypto.Lol.LatticePrelude++, module Crypto.Lol.Types.ZqBasic+, module Crypto.Lol.Cyclotomic.Tensor.CTensor+, module Crypto.Lol.Cyclotomic.Tensor.RepaTensor+, module Crypto.Lol.Types.IrreducibleChar2) where++import Crypto.Lol.Cyclotomic.Cyc+import Crypto.Lol.Gadget+import Crypto.Lol.LatticePrelude++import Crypto.Lol.Types.ZqBasic+import Crypto.Lol.Cyclotomic.Tensor.CTensor+import Crypto.Lol.Cyclotomic.Tensor.RepaTensor+import Crypto.Lol.Types.IrreducibleChar2

@@ -0,0 +1,138 @@+{-# LANGUAGE FlexibleContexts, FlexibleInstances, RebindableSyntax,
+             ScopedTypeVariables, TypeFamilies #-}
+
+-- | Classes and helper methods for the Chinese remainder transform
+-- and ring extensions.
+
+module Crypto.Lol.CRTrans
+( CRTrans(..), CRTEmbed(..)
+, CRTInfo
+, crtInfoFact, crtInfoPPow, crtInfoPrime
+, gEmbPPow, gEmbPrime
+) where
+
+import Crypto.Lol.LatticePrelude
+
+import Control.Arrow
+import Data.Singletons
+import Data.Singletons.Prelude
+
+-- | Information that characterizes the (invertible) Chinese remainder
+-- transformation over a ring @r@, namely:
+--
+--     (1) a function that returns the @i@th power of some /principal/
+--     @m@th root of unity (for any integer @i@)
+--
+--     (2) the multiplicative inverse of @\\hat{m}@ in @r@.
+
+type CRTInfo r = (Int -> r, r)
+
+-- | A ring that (possibly) supports invertible Chinese remainder
+-- transformations of various indices.
+
+-- | The values of 'crtInfo' for different indices @m@ should be
+-- consistent, in the sense that if @omega@, @omega'@ are respectively
+-- the values returned for @m@, @m'@ where @m'@ divides @m@, then it
+-- should be the case that @omega^(m/m')=omega'@.
+
+class Ring r => CRTrans r where
+
+  -- | 'CRTInfo' for a given index @m@. The method itself may be
+  -- slow, but the function it returns should be fast, e.g., via
+  -- internal memoization.  The default implementation returns
+  -- 'Nothing'.
+  crtInfo :: Int -> Maybe (CRTInfo r)
+  crtInfo = const Nothing
+
+-- | A ring with a ring embedding into some ring @CRTExt r@ that has
+-- an invertible CRT transformation for /every/ positive index @m@.
+class (Ring r, Ring (CRTExt r)) => CRTEmbed r where
+  type CRTExt r
+
+  -- | Embeds from @r@ to @CRTExt r@
+  toExt :: r -> CRTExt r
+  -- | Projects from @CRTExt r@ to @r@
+  fromExt :: CRTExt r -> r
+
+-- CRTrans instance for product rings
+instance (CRTrans a, CRTrans b) => CRTrans (a,b) where
+  crtInfo i = do
+    (apow, aiInv) <- crtInfo i
+    (bpow, biInv) <- crtInfo i
+    return (apow &&& bpow, (aiInv, biInv))
+
+-- CRTEmbed instance for product rings
+instance (CRTEmbed a, CRTEmbed b) => CRTEmbed (a,b) where
+  type CRTExt (a,b) = (CRTExt a, CRTExt b)
+  toExt = toExt *** toExt
+  fromExt = fromExt *** fromExt
+
+omegaPowC :: (Transcendental a) => Int -> Int -> Complex a
+omegaPowC m i = cis (2*pi*fromIntegral i / fromIntegral m)
+
+-- | 'crtInfo' wrapper for 'Fact' types.
+crtInfoFact :: (Fact m, CRTrans r) => TaggedT m Maybe (CRTInfo r)
+crtInfoFact = (tagT . crtInfo) =<< pureT valueFact
+
+-- | 'crtInfo' wrapper for 'PPow' types.
+crtInfoPPow :: (PPow pp, CRTrans r) => TaggedT pp Maybe (CRTInfo r)
+crtInfoPPow = (tagT . crtInfo) =<< pureT valuePPow
+
+-- | 'crtInfo' wrapper for 'Prime' types.
+crtInfoPrime :: (Prim p, CRTrans r) => TaggedT p Maybe (CRTInfo r)
+crtInfoPrime = (tagT . crtInfo) =<< pureT valuePrime
+
+-- | A function that returns the 'i'th embedding of @g_{p^e} = g_p@ for
+-- @i@ in @Z*_{p^e}@.
+gEmbPPow :: forall pp r . (PPow pp, CRTrans r) => TaggedT pp Maybe (Int -> r)
+gEmbPPow = tagT $ case (sing :: SPrimePower pp) of
+  (SPP (STuple2 sp _)) -> withWitnessT gEmbPrime sp
+
+-- | A function that returns the @i@th embedding of @g_p@ for @i@ in @Z*_p@,
+-- i.e., @1-omega_p^i@.
+gEmbPrime :: (Prim p, CRTrans r) => TaggedT p Maybe (Int -> r)
+gEmbPrime = do
+  (f, _) <- crtInfoPrime
+  return $ \i -> one - f i      -- not checking that i /= 0 (mod p)
+
+-- the complex numbers have roots of unity of any order
+instance (Transcendental a) => CRTrans (Complex a) where
+  crtInfo m = Just (omegaPowC m, recip $ fromIntegral $ valueHat m)
+
+-- trivial CRTEmbed instance for complex numbers
+instance (Transcendental a) => CRTEmbed (Complex a) where
+  type CRTExt (Complex a) = Complex a
+  toExt = id
+  fromExt = id
+
+-- Default CRTrans instances for real and integer types, which do
+-- not have roots of unity (except in trivial cases). These are needed
+-- to use FastCyc with these integer types.
+instance CRTrans Double
+instance CRTrans Int
+instance CRTrans Int64
+instance CRTrans Integer
+-- can also do for Int8, Int16, Int32 etc.
+
+-- CRTEmbed instances for real and integer types, embedding into
+-- Complex.  These are needed to use FastCyc with these integer types.
+instance CRTEmbed Double where
+  type CRTExt Double = Complex Double
+  toExt = fromReal . realToField
+  fromExt = realToField . real
+
+instance CRTEmbed Int where
+  type CRTExt Int = Complex Double
+  toExt = fromIntegral
+  fromExt = fst . roundComplex
+
+instance CRTEmbed Int64 where
+  type CRTExt Int64 = Complex Double
+  toExt = fromIntegral
+  fromExt = fst . roundComplex
+
+instance CRTEmbed Integer where
+  -- CJP: sufficient precision?  Not in general.
+  type CRTExt Integer = Complex Double
+  toExt = fromIntegral
+  fromExt = fst . roundComplex

@@ -0,0 +1,583 @@+{-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts,+             FlexibleInstances, GADTs, MultiParamTypeClasses,+             NoImplicitPrelude, PolyKinds, RankNTypes, ScopedTypeVariables,+             TypeFamilies, TypeOperators, UndecidableInstances #-}++-- | An implementation of cyclotomic rings that hides the+-- internal representations of ring elements (e.g., the choice of+-- basis), and also offers more efficient storage and operations on+-- subring elements (including elements from the base ring itself).+--+-- For an implementation that allows (and requires) the programmer to+-- control the underlying representation, see+-- 'Crypto.Lol.Cyclotomic.UCyc.UCyc'.+--+-- __WARNING:__ as with all fixed-point arithmetic, the functions+-- associated with 'Cyc' may result in overflow (and thereby+-- incorrect answers and potential security flaws) if the input+-- arguments are too close to the bounds imposed by the base type.+-- The acceptable range of inputs for each function is determined by+-- the internal linear transforms and other operations it performs.++module Crypto.Lol.Cyclotomic.Cyc+(+-- * Data type and constraints+  Cyc, CElt, U.NFElt+-- * Constructors/deconstructors+, cycPow, cycDec, cycCRT, scalarCyc+, uncycPow, uncycDec, uncycCRT, unzipCyc, unzipCElt+-- * Core cyclotomic operations+, mulG, divG, gSqNorm, liftCyc, liftPow, liftDec+-- * Representation advice+, advisePow, adviseDec, adviseCRT+-- * Error sampling+, tGaussian, errorRounded, errorCoset+-- * Sub/extension rings+, embed, twace, coeffsCyc, powBasis, crtSet+-- * Rescaling cyclotomic elements+, R.RescaleCyc(..), R.Basis(..)+) where++import qualified Algebra.Additive     as Additive (C)+import qualified Algebra.Module       as Module (C)+import qualified Algebra.Ring         as Ring (C)+import qualified Algebra.ZeroTestable as ZeroTestable (C)++import Crypto.Lol.Cyclotomic.UCyc hiding (coeffsDec, coeffsPow, crtSet,+                                   divG, errorCoset, errorRounded, gSqNorm,+                                   mulG, powBasis, tGaussian, unzipCyc)++import           Crypto.Lol.CRTrans+import qualified Crypto.Lol.Cyclotomic.RescaleCyc as R+import           Crypto.Lol.Cyclotomic.Tensor     (CRTElt, TElt, Tensor)+import qualified Crypto.Lol.Cyclotomic.UCyc       as U+import           Crypto.Lol.Gadget+import           Crypto.Lol.LatticePrelude        as LP+import           Crypto.Lol.Types.FiniteField+import           Crypto.Lol.Types.ZPP++import Control.Applicative    hiding ((*>))+import Control.Arrow+import Control.DeepSeq+import Control.Monad.Identity   -- needed for coerce+import Control.Monad.Random+import Data.Coerce+import Data.Traversable++import Test.QuickCheck++-- | Represents a cyclotomic ring such as @Z[zeta]@,+-- @Zq[zeta]@, and @Q(zeta)@ in an explicit representation: @t@ is the+-- 'Tensor' type for storing coefficient tensors; @m@ is the+-- cyclotomic index; @r@ is the base ring of the coefficients (e.g.,+-- @Z@, @Zq@).+data Cyc t m r where+  Pow :: !(UCyc t m P r) -> Cyc t m r+  Dec :: !(UCyc t m D r) -> Cyc t m r+  CRT :: !(UCyc t m C r) -> Cyc t m r+  -- super-optimized storage of scalars+  Scalar :: !r -> Cyc t m r+  -- optimized storage of subring elements+  Sub :: (l `Divides` m) => !(Cyc t l r) -> Cyc t m r+  -- CJP: someday try to merge the above two++-- | Constraints needed for many operations involving 'Cyc' data.+type CElt t r = UCElt t r++---------- Constructors / deconstructors ----------++-- | Wrap a 'UCyc' as a 'Cyc'.+cycPow :: UCyc t m P r -> Cyc t m r+cycPow = Pow+{-# INLINABLE cycPow #-}++-- | Wrap a 'UCyc' as a 'Cyc'.+cycDec :: UCyc t m D r -> Cyc t m r+cycDec = Dec+{-# INLINABLE cycDec #-}++-- | Wrap a 'UCyc' as a 'Cyc'.+cycCRT :: UCyc t m C r -> Cyc t m r+cycCRT = CRT+{-# INLINABLE cycCRT #-}++-- | Embed a scalar from the base ring as a cyclotomic element.+scalarCyc :: r -> Cyc t m r+scalarCyc = Scalar+{-# INLINABLE scalarCyc #-}++-- | Unwrap a 'Cyc' as a 'UCyc' in powerful-basis representation.+uncycPow :: (Fact m, CElt t r) => Cyc t m r -> UCyc t m P r+{-# INLINABLE uncycPow #-}+uncycPow c = let (Pow u) = toPow' c in u++-- | Unwrap a 'Cyc' as a 'UCyc' in decoding-basis representation.+uncycDec :: (Fact m, CElt t r) => Cyc t m r -> UCyc t m D r+{-# INLINABLE uncycDec #-}+uncycDec c = let (Dec u) = toDec' c in u++-- | Unwrap a 'Cyc' as a 'UCyc' in CRT-basis representation.+uncycCRT :: (Fact m, CElt t r) => Cyc t m r -> UCyc t m C r+{-# INLINABLE uncycCRT #-}+uncycCRT c = let (CRT u) = toCRT' c in u++---------- Algebraic instances ----------++instance (Fact m, CElt t r) => ZeroTestable.C (Cyc t m r) where+  isZero (Pow u) = isZero u+  isZero (Dec u) = isZero u+  isZero (CRT u) = isZero u+  isZero (Scalar c) = isZero c+  isZero (Sub c) = isZero c+  {-# INLINABLE isZero #-}++instance (Eq r, Fact m, CElt t r) => Eq (Cyc t m r) where+  -- same representations+  (Scalar c1) == (Scalar c2) = c1 == c2+  (Pow u1) == (Pow u2) = u1 == u2+  (Dec u1) == (Dec u2) = u1 == u2+  (CRT u1) == (CRT u2) = u1 == u2+  -- compare Subs in compositum+  (Sub (c1 :: Cyc t l1 r)) == (Sub (c2 :: Cyc t l2 r)) =+    (embed' c1 :: Cyc t (FLCM l1 l2) r) == embed' c2+    \\ lcmDivides (Proxy::Proxy l1) (Proxy::Proxy l2)++  -- some other relatively efficient comparisons+  (Scalar c1) == (Pow u2) = scalarPow c1 == u2+  -- (Scalar c1) == (Dec u2) = scalarDec c1 == u2+  (Pow u1) == (Scalar c2) = u1 == scalarPow c2+  -- (Dec u1) == (Scalar c2) = u1 == scalarDec c2++  -- default: compare in powerful basis+  c1 == c2 = toPow' c1 == toPow' c2++  {-# INLINABLE (==) #-}++instance (Fact m, CElt t r) => Additive.C (Cyc t m r) where+  {-# INLINABLE zero #-}+  zero = Scalar zero++  {-# INLINABLE (+) #-}+  -- optimized addition of zero+  (Scalar c1) + c2 | isZero c1 = c2+  c1 + (Scalar c2) | isZero c2 = c1++  -- SAME CONSTRUCTORS+  (Scalar c1) + (Scalar c2) = Scalar (c1+c2)+  (Pow u1) + (Pow u2) = Pow $ u1 + u2+  (Dec u1) + (Dec u2) = Dec $ u1 + u2+  (CRT u1) + (CRT u2) = CRT $ u1 + u2+  -- Sub plus Sub: work in compositum+  (Sub (c1 :: Cyc t m1 r)) + (Sub (c2 :: Cyc t m2 r)) =+    (Sub $ (embed' c1 :: Cyc t (FLCM m1 m2) r) + embed' c2)+    \\ lcm2Divides (Proxy::Proxy m1) (Proxy::Proxy m2) (Proxy::Proxy m)++  -- SCALAR PLUS SOMETHING ELSE++  (Scalar c) + (Pow u) = Pow $ scalarPow c + u+  (Scalar c) + (Dec u) = Pow $ scalarPow c + toPow u -- workaround scalarDec+  (Scalar c) + (CRT u) = CRT $ scalarCRT c + u+  (Scalar c1) + (Sub c2) = Sub $ Scalar c1 + c2 -- must re-wrap Scalar!++  (Pow u) + (Scalar c) = Pow $ u + scalarPow c+  (Dec u) + (Scalar c) = Pow $ toPow u + scalarPow c -- workaround scalarDec+  (CRT u) + (Scalar c) = CRT $ u + scalarCRT c+  (Sub c1) + (Scalar c2) = Sub $ c1 + Scalar c2++  -- SUB PLUS NON-SUB, NON-SCALAR: work in full ring+  (Sub c1) + c2 = embed' c1 + c2+  c1 + (Sub c2) = c1 + embed' c2++  -- mixed Dec and Pow: use linear-time conversions+  (Dec u1) + (Pow u2) = Pow $ toPow u1 + u2+  (Pow u1) + (Dec u2) = Pow $ u1 + toPow u2++  -- one CRT: convert other to CRT+  (CRT u1) + (Pow u2) = CRT $ u1 + toCRT u2+  (CRT u1) + (Dec u2) = CRT $ u1 + toCRT u2+  (Pow u1) + (CRT u2) = CRT $ toCRT u1 + u2+  (Dec u1) + (CRT u2) = CRT $ toCRT u1 + u2++  {-# INLINABLE negate #-}+  negate (Pow u) = Pow $ negate u+  negate (Dec u) = Dec $ negate u+  negate (CRT u) = CRT $ negate u+  negate (Scalar c) = Scalar (negate c)+  negate (Sub c) = Sub $ negate c++instance (Fact m, CElt t r) => Ring.C (Cyc t m r) where+  {-# INLINABLE one #-}+  one = Scalar one++  {-# INLINABLE fromInteger #-}+  fromInteger = Scalar . fromInteger++  {-# INLINABLE (*) #-}++  -- optimized mul-by-zero+  v1@(Scalar c1) * _ | isZero c1 = v1+  _ * v2@(Scalar c2) | isZero c2 = v2++  -- both CRT+  (CRT u1) * (CRT u2) = CRT $ u1*u2++  -- at least one Scalar+  (Scalar c1) * (Scalar c2) = Scalar $ c1*c2+  (Scalar c) * (Pow u) = Pow $ c *> u+  (Scalar c) * (Dec u) = Dec $ c *> u+  (Scalar c) * (CRT u) = CRT $ c *> u+  (Scalar c1) * (Sub c2) = Sub $ Scalar c1 * c2++  (Pow u) * (Scalar c) = Pow $ c *> u+  (Dec u) * (Scalar c) = Dec $ c *> u+  (CRT u) * (Scalar c) = CRT $ c *> u+  (Sub c1) * (Scalar c2) = Sub $ c1 * Scalar c2++  -- TWO SUBS: work in a CRT rep for compositum+  (Sub (c1 :: Cyc t m1 r)) * (Sub (c2 :: Cyc t m2 r)) =+    -- re-wrap c1, c2 as Subs of the composition, and force them to CRT+    (Sub $ (toCRT' $ Sub c1 :: Cyc t (FLCM m1 m2) r) * toCRT' (Sub c2))+    \\ lcm2Divides (Proxy::Proxy m1) (Proxy::Proxy m2) (Proxy::Proxy m)++  -- ELSE: work in appropriate CRT rep+  c1 * c2 = toCRT' c1 * toCRT' c2++instance (GFCtx fp d, Fact m, CElt t fp) => Module.C (GF fp d) (Cyc t m fp) where+  -- CJP: optimize for Scalar if we can: r *> (Scalar c) is the tensor+  -- that has the coeffs of (r*c), followed by zeros.  (This assumes+  -- that the powerful basis has 1 as its first element, and that+  -- we're using pow to define the module mult.)++  -- can use any r-basis to define module mult, but must be+  -- consistent.+  r *> (Pow v) = Pow $ r LP.*> v+  r *> x = r *> toPow' x++---------- Core cyclotomic operations ----------++advisePow, adviseDec, adviseCRT :: (Fact m, CElt t r) => Cyc t m r -> Cyc t m r+{-# INLINABLE advisePow #-}+{-# INLINABLE adviseDec #-}+{-# INLINABLE adviseCRT #-}++-- | Same as 'adviseCRT', but for the powerful-basis representation.+advisePow = toPow'++-- | Same as 'adviseCRT', but for the powerful-basis representation.+adviseDec = toDec'++-- | Yield an equivalent element that /may/ be in a CRT+-- representation.  This can serve as an optimization hint. E.g.,+-- call 'adviseCRT' prior to multiplying the same value by many+-- other values.+adviseCRT = toCRT'++-- | Multiply by the special element @g@ of the @m@th cyclotomic.+mulG :: (Fact m, CElt t r) => Cyc t m r -> Cyc t m r+{-# INLINABLE mulG #-}+mulG (Pow u) = Pow $ U.mulG u+mulG (Dec u) = Dec $ U.mulG u+mulG (CRT u) = CRT $ U.mulG u+mulG (Scalar r) = CRT $ U.mulG $ scalarCRT r+mulG (Sub c) = mulG $ embed' c   -- must go to full ring++-- | Divide by @g@, returning 'Nothing' if not evenly divisible.+-- WARNING: this implementation is not a constant-time algorithm, so+-- information about the argument may be leaked through a timing+-- channel.+divG :: (Fact m, CElt t r) => Cyc t m r -> Maybe (Cyc t m r)+{-# INLINABLE divG #-}+divG (Pow u) = Pow <$> U.divG u+divG (Dec u) = Dec <$> U.divG u+divG (CRT u) = CRT <$> U.divG u+divG (Scalar r) = CRT <$> U.divG (scalarCRT r)+divG (Sub c) = divG $ embed' c  -- must go to full ring++-- | Sample from the "tweaked" Gaussian error distribution @t*D@ in+-- the decoding basis, where @D@ has scaled variance @v@.  (Note: This+-- implementation uses Double precision to generate the Gaussian+-- sample, which may not be sufficient for rigorous proof-based+-- security.)+tGaussian :: (Fact m, OrdFloat q, Random q, Tensor t, TElt t q,+              ToRational v, MonadRandom rnd)+             => v -> rnd (Cyc t m q)+tGaussian = (Dec <$>) . U.tGaussian+{-# INLINABLE tGaussian #-}++-- | Yield the scaled squared norm of @g_m \cdot e@ under+-- the canonical embedding, namely,+-- @\hat{m}^{ -1 } \cdot || \sigma(g_m \cdot e) ||^2@ .+gSqNorm :: forall t m r . (Fact m, CElt t r) => Cyc t m r -> r+{-# INLINABLE gSqNorm #-}+gSqNorm (Pow u) = U.gSqNorm $ toDec u+gSqNorm (Dec u) = U.gSqNorm u+gSqNorm (CRT u) = U.gSqNorm $ toDec u+-- CJP: this is suboptimal (gSqNorm of scalar 1 is known)+-- also workaround scalarDec+gSqNorm (Scalar u) = U.gSqNorm (toDec $ scalarPow u :: UCyc t m D r)+gSqNorm (Sub c) = U.gSqNorm (U.embedDec $ uncycDec c :: UCyc t m D r)++-- | Generate an LWE error term with given scaled variance,+-- deterministically rounded with respect to the decoding basis.+errorRounded :: (ToInteger z, Tensor t, Fact m, TElt t z,+                 ToRational v, MonadRandom rnd) => v -> rnd (Cyc t m z)+{-# INLINABLE errorRounded #-}+errorRounded = (Dec <$>) . U.errorRounded++-- | Generate an LWE error term with given scaled variance @* p^2@ over+-- the given coset, deterministically rounded with respect to the+-- decoding basis.+errorCoset ::+  (Mod zp, z ~ ModRep zp, Lift zp z, Fact m,+   CElt t zp, TElt t z, ToRational v, MonadRandom rnd)+  => v -> Cyc t m zp -> rnd (Cyc t m z)+errorCoset v = (Dec <$>) . U.errorCoset v . uncycDec+{-# INLINABLE errorCoset #-}++---------- Inter-ring operations ----------++-- | Embed (lazily) into an extension ring.+embed :: forall t m m' r . (m `Divides` m') => Cyc t m r -> Cyc t m' r+{-# INLINABLE embed #-}+embed (Scalar c) = Scalar c           -- keep as scalar+embed (Sub (c :: Cyc t l r)) = Sub c  -- keep as subring element+  \\ transDivides (Proxy::Proxy l) (Proxy::Proxy m) (Proxy::Proxy m')+embed c = Sub c++-- | Force to a non-'Sub' constructor (for internal use only).+embed' :: forall t r l m . (l `Divides` m, CElt t r) => Cyc t l r -> Cyc t m r+{-# INLINE embed' #-}+embed' (Pow u) = Pow $ embedPow u+embed' (Dec u) = Dec $ embedDec u+embed' (CRT u) = either Pow CRT $ embedCRT u+embed' (Scalar c) = Scalar c+embed' (Sub (c :: Cyc t k r)) = embed' c+  \\ transDivides (Proxy::Proxy k) (Proxy::Proxy l) (Proxy::Proxy m)++-- | The "tweaked trace" (twace) function+-- @Tw(x) = (mhat \/ m'hat) * Tr(g' \/ g * x)@,+-- which fixes @R@ pointwise (i.e., @twace . embed == id@).+twace :: forall t m m' r . (m `Divides` m', CElt t r)+         => Cyc t m' r -> Cyc t m r+{-# INLINABLE twace #-}+twace (Pow u) = Pow $ U.twacePow u+twace (Dec u) = Dec $ U.twaceDec u+twace (CRT u) = either Pow CRT $ twaceCRT u+twace (Scalar u) = Scalar u+twace (Sub (c :: Cyc t l r)) = Sub (twace c :: Cyc t (FGCD l m) r)+                               \\ gcdDivides (Proxy::Proxy l) (Proxy::Proxy m)++-- | Return the given element's coefficient vector with respect to+-- the (relative) powerful/decoding basis of the cyclotomic+-- extension @O_m' / O_m@.+coeffsCyc :: (m `Divides` m', CElt t r) => R.Basis -> Cyc t m' r -> [Cyc t m r]+{-# INLINABLE coeffsCyc #-}+coeffsCyc R.Pow c' = Pow <$> U.coeffsPow (uncycPow c')+coeffsCyc R.Dec c' = Dec <$> U.coeffsDec (uncycDec c')++-- | The relative powerful basis of @O_m' / O_m@.+powBasis :: (m `Divides` m', CElt t r) => Tagged m [Cyc t m' r]+powBasis = (Pow <$>) <$> U.powBasis+{-# INLINABLE powBasis #-}++-- | The relative mod-@r@ CRT set of the extension.+crtSet :: (m `Divides` m', ZPP r, CElt t r, TElt t (ZpOf r))+          => Tagged m [Cyc t m' r]+crtSet = (Pow <$>) <$> U.crtSet+{-# INLINABLE crtSet #-}+++---------- Lattice operations and instances ----------++instance (Reduce a b, Fact m, CElt t a, CElt t b)+    -- CJP: need these specific constraints to get Reduce instance for Sub case+         => Reduce (Cyc t m a) (Cyc t m b) where+  {-# INLINABLE reduce #-}+  reduce (Pow u) = Pow $ reduce u+  reduce (Dec u) = Dec $ reduce u+  reduce (CRT u) = Pow $ reduce $ toPow u+  reduce (Scalar c) = Scalar $ reduce c+  reduce (Sub (c :: Cyc t l a)) = Sub (reduce c :: Cyc t l b)++type instance LiftOf (Cyc t m r) = Cyc t m (LiftOf r)++-- | Lift using the specified basis.+liftCyc :: (Lift b a, Fact m, TElt t a, CElt t b)+           => R.Basis -> Cyc t m b -> Cyc t m a+{-# INLINABLE liftCyc #-}++liftCyc R.Pow = liftPow+liftCyc R.Dec = liftDec++liftPow, liftDec :: (Lift b a, Fact m, TElt t a, CElt t b)+                    => Cyc t m b -> Cyc t m a+{-# INLINABLE liftPow #-}+{-# INLINABLE liftDec #-}++-- | Lift using the powerful basis.+liftPow (Pow u) = Pow $ lift u+liftPow (Dec u) = Pow $ lift $ toPow u+liftPow (CRT u) = Pow $ lift $ toPow u+-- optimized for subrings; these are correct for powerful basis but+-- not for decoding+liftPow (Scalar c) = Scalar $ lift c+liftPow (Sub c) = Sub $ liftPow c++-- | Lift using the decoding basis.+liftDec (Pow u) = Dec $ lift $ toDec u+liftDec (Dec u) = Dec $ lift u+liftDec (CRT u) = Dec $ lift $ toDec u+liftDec (Scalar c) = Dec $ lift $ toDec $ scalarPow c -- workaround scalarDec+liftDec (Sub c) = liftDec $ embed' c++-- | Unzip for unrestricted types.+unzipCyc :: (Tensor t, Fact m) => Cyc t m (a,b) -> (Cyc t m a, Cyc t m b)+{-# INLINABLE unzipCyc #-}+unzipCyc (Pow u) = Pow *** Pow $ U.unzipCyc u+unzipCyc (Dec u) = Dec *** Dec $ U.unzipCyc u+unzipCyc (CRT u) = CRT *** CRT $ U.unzipCyc u+unzipCyc (Scalar c) = Scalar *** Scalar $ c+unzipCyc (Sub c) = Sub *** Sub $ unzipCyc c++-- | Type-restricted (and potentially more efficient) unzip.+unzipCElt :: (Tensor t, Fact m, CElt t (a,b), CElt t a, CElt t b)+             => Cyc t m (a,b) -> (Cyc t m a, Cyc t m b)+{-# INLINABLE unzipCElt #-}+unzipCElt (Pow u) = Pow *** Pow $ U.unzipUCElt u+unzipCElt (Dec u) = Dec *** Dec $ U.unzipUCElt u+unzipCElt (CRT u) = CRT *** CRT $ U.unzipUCElt u+unzipCElt (Scalar c) = Scalar *** Scalar $ c+unzipCElt (Sub c) = Sub *** Sub $ unzipCElt c++-- generic RescaleCyc instance++instance {-# OVERLAPS #-} (Rescale a b, CElt t a, TElt t b)+    => R.RescaleCyc (Cyc t) a b where++  -- Optimized for subring constructors, for powerful basis.+  -- Analogs for decoding basis are not quite correct, because (* -1)+  -- doesn't commute with 'rescale' due to tiebreakers!+  rescaleCyc R.Pow (Scalar c) = Scalar $ rescale c+  rescaleCyc R.Pow (Sub c) = Sub $ R.rescalePow c++  rescaleCyc R.Pow c = Pow $ fmapPow rescale $ uncycPow c+  rescaleCyc R.Dec c = Dec $ fmapDec rescale $ uncycDec c+  {-# INLINABLE rescaleCyc #-}++-- specialized instance for product rings: ~2x faster algorithm+instance (Mod a, Field b, Lift a (ModRep a), Reduce (LiftOf a) b,+         CElt t (a,b), CElt t a, CElt t b, CElt t (LiftOf a))+         => R.RescaleCyc (Cyc t) (a,b) b where++  -- optimized for subrings and powerful basis (see comments in other+  -- instance for why this doesn't work for decoding basis)+  rescaleCyc R.Pow (Scalar c) = Scalar $ rescale c+  rescaleCyc R.Pow (Sub c) = Sub $ R.rescalePow c++  rescaleCyc bas c = let aval = proxy modulus (Proxy::Proxy a)+                         (a,b) = unzipCElt c+                         z = liftCyc bas a+                     in Scalar (recip (reduce aval)) * (b - reduce z)+  {-# INLINABLE rescaleCyc #-}+++instance (Gadget gad zq, Fact m, CElt t zq) => Gadget gad (Cyc t m zq) where+  gadget = (scalarCyc <$>) <$> gadget+  -- specialization fo 'encode', done efficiently+  encode s = ((* adviseCRT s) <$>) <$> gadget+  {-# INLINABLE gadget #-}+  {-# INLINABLE encode #-}++-- promote Decompose, using the powerful basis+instance (Decompose gad zq, Fact m, CElt t zq, CElt t (DecompOf zq))+         => Decompose gad (Cyc t m zq) where++  type DecompOf (Cyc t m zq) = Cyc t m (DecompOf zq)++  -- faster implementations: decompose directly in subring, which is+  -- correct because we decompose in powerful basis+  decompose (Scalar c) = pasteT $ Scalar <$> peelT (decompose c)+  decompose (Sub c) = pasteT $ Sub <$> peelT (decompose c)++  -- traverse: Traversable (UCyc t m P) and Applicative (Tagged gad ZL)+  decompose (Pow u) = fromZL $ Pow <$> traverse (toZL . decompose) u+  decompose c = decompose $ toPow' c++  {-# INLINABLE decompose #-}++toZL :: Tagged s [a] -> TaggedT s ZipList a+toZL = coerce++fromZL :: TaggedT s ZipList a -> Tagged s [a]+fromZL = coerce++-- promote Correct, using the decoding basis+instance (Correct gad zq, Fact m, CElt t zq) => Correct gad (Cyc t m zq) where+  -- sequence: Monad [] and Traversable (UCyc t m D)+  -- sequenceA: Applicative (UCyc t m D) and Traversable (TaggedT gad [])+  correct bs = Dec *** (Dec <$>) $+               second sequence $ U.unzipCyc $ (correct . pasteT) <$>+               sequenceA (uncycDec <$> peelT bs)+  {-# INLINABLE correct #-}++---------- Change of representation (internal use only) ----------++toPow', toDec', toCRT' :: (Fact m, CElt t r) => Cyc t m r -> Cyc t m r+{-# INLINE toPow' #-}+{-# INLINE toDec' #-}+{-# INLINE toCRT' #-}++-- | Force to powerful-basis representation (for internal use only).+toPow' c@(Pow _) = c+toPow' (Dec u) = Pow $ toPow u+toPow' (CRT u) = Pow $ toPow u+toPow' (Scalar c) = Pow $ scalarPow c+toPow' (Sub c) = toPow' $ embed' c++-- | Force to decoding-basis representation (for internal use only).+toDec' (Pow u) = Dec $ toDec u+toDec' c@(Dec _) = c+toDec' (CRT u) = Dec $ toDec u+toDec' (Scalar c) = Dec $ toDec $ scalarPow c -- workaround scalarDec+toDec' (Sub c) = toDec' $ embed' c++-- | Force to CRT representation (for internal use only).+toCRT' (Pow u) = CRT $ toCRT u+toCRT' (Dec u) = CRT $ toCRT u+toCRT' c@(CRT _) = c+toCRT' (Scalar c) = CRT $ scalarCRT c+-- CJP: below is the fastest algorithm for when both source and target+-- have the same CRTr/CRTe choice.  It is not the fastest when the+-- choices are different (it will do an unnecessary CRT if input is+-- non-CRT), but this is an unusual case.  Note: both toCRT' are+-- necessary in generaly, because embed' may not preserve CRT+-- representation!+toCRT' (Sub c) = toCRT' $ embed' $ toCRT' c++---------- Utility instances ----------++instance (Tensor t, Fact m, NFData r, TElt t r,+          NFData (CRTExt r), TElt t (CRTExt r)) => NFData (Cyc t m r) where+  rnf (Pow u) = rnf u+  rnf (Dec u) = rnf u+  rnf (CRT u) = rnf u+  rnf (Scalar u) = rnf u+  rnf (Sub c) = rnf c++instance (Random r, Tensor t, Fact m, CRTElt t r) => Random (Cyc t m r) where+  random g = let (u,g') = random g+             in (either Pow CRT u, g')+  {-# INLINABLE random #-}++  randomR _ = error "randomR non-sensical for Cyc"++instance (Arbitrary (UCyc t m P r)) => Arbitrary (Cyc t m r) where+  arbitrary = Pow <$> arbitrary+  shrink = shrinkNothing++instance (Show r, Show (CRTExt r), Tensor t, Fact m, TElt t r, TElt t (CRTExt r)) => Show (Cyc t m r) where+  show (Scalar c) = "Cyc Scalar: " ++ show c+  show (Pow u) = "Cyc: " ++ show u+  show (Dec u) = "Cyc: " ++ show u+  show (CRT u) = "Cyc: " ++ show u+  show (Sub c) = "Cyc Sub: " ++ show c

@@ -0,0 +1,90 @@+{-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts,+             GeneralizedNewtypeDeriving, KindSignatures,+             MultiParamTypeClasses, NoImplicitPrelude, RoleAnnotations,+             ScopedTypeVariables, StandaloneDeriving, TypeFamilies,+             TypeOperators, UndecidableInstances #-}++-- | Functions from one cyclotomic ring to another that are linear+-- over a common subring.++module Crypto.Lol.Cyclotomic.Linear+( Linear, ExtendLinIdx+, linearDec, evalLin, extendLin+) where++import Crypto.Lol.Cyclotomic.Cyc+import Crypto.Lol.LatticePrelude++import Algebra.Additive as Additive (C)++import Control.Applicative+import Control.DeepSeq++-- | An @E@-linear function from @R@ to @S@.++-- CJP: also have constructor for relative Pow basis of R/E?  So far+-- not needed.+newtype Linear t z (e::Factored) (r::Factored) (s::Factored) = RD [Cyc t s z]++deriving instance (NFData (Cyc t s z)) => NFData (Linear t z e r s)++-- some params are phantom but matter for safety+type role Linear representational nominal representational representational nominal++-- | Construct an @E@-linear function given a list of its output values+-- (in @S@) on the relative decoding basis of @R/E@.  The number of+-- elements in the list must not exceed the size of the basis.+linearDec :: forall t z e r s .+             (e `Divides` r, e `Divides` s, CElt t z)+             => [Cyc t s z] -> Linear t z e r s+linearDec ys = let ps = proxy powBasis (Proxy::Proxy e) `asTypeOf` ys+               in if length ys <= length ps then RD (adviseCRT <$> ys)+               else error $ "linearDec: too many entries: "+                        ++ show (length ys) ++ " versus "+                        ++ show (length ps)++-- | Evaluates the given linear function on the input.+evalLin :: forall t z e r s .+           (e `Divides` r, e `Divides` s, CElt t z)+           => Linear t z e r s -> Cyc t r z -> Cyc t s z+evalLin (RD ys) r = sum (zipWith (*) ys $+                         embed <$> (coeffsCyc Dec r :: [Cyc t e z]))++instance Additive (Cyc t s z) => Additive.C (Linear t z e r s) where+  zero = RD []++  (RD as) + (RD bs) = RD $ sumall as bs+    where sumall [] ys = ys+          sumall xs [] = xs+          sumall (x:xs) (y:ys) = x+y : sumall xs ys++  negate (RD as) = RD $ negate <$> as++instance (Reduce z zq, Fact s, CElt t z, CElt t zq)+         => Reduce (Linear t z e r s) (Linear t zq e r s) where+  reduce (RD ys) = RD $ reduce <$> ys++type instance LiftOf (Linear t zp e r s) = Linear t (LiftOf zp) e r s++instance (CElt t zp, CElt t z, z ~ LiftOf zp, Lift zp z, Fact s)+         => Lift' (Linear t zp e r s) where+  lift (RD ys) = RD $ liftCyc Pow <$> ys++-- | A convenient constraint synonym for extending a linear function+-- to larger rings.+type ExtendLinIdx e r s e' r' s' =+  (Fact r, e ~ FGCD r e', r' ~ FLCM r e', -- these imply R'=R\otimes_E E'+   e' `Divides` s', s `Divides` s') -- lcm(s,e')|s' <=> (S+E') \subseteq S'++-- | Extend an @E@-linear function @R->S@ to an @E'@-linear function+-- @R\'->S\'@.  (Mathematically, such extension only requires+-- @lcm(r,e\') | r\'@ (not equality), but this generality would+-- significantly complicate the implementation, and for our purposes+-- there's no reason to use any larger @r'@.)+extendLin :: (ExtendLinIdx e r s e' r' s', CElt t z)+           => Linear t z e r s -> Linear t z e' r' s'+-- CJP: this simple implementation works because R/E and R'/E' have+-- identical decoding bases, because R' \cong R \otimes_E E'.  If we+-- relax the constraint on E then we'd have to change the+-- implementation to something more difficult.+extendLin (RD ys) = RD (embed <$> ys)

@@ -0,0 +1,23 @@+{-# LANGUAGE MultiParamTypeClasses #-}++module Crypto.Lol.Cyclotomic.RescaleCyc where++import Crypto.Lol.Factored++data Basis = Pow | Dec++-- | Represents cyclotomic rings that are rescalable over their base+-- rings.  (This is a class because it allows for more efficient+-- specialized implementations.)++class RescaleCyc c a b where+  -- | Rescale in the given basis.+  rescaleCyc :: Fact m => Basis -> c m a -> c m b++-- | Specialized convenience functions.+rescalePow, rescaleDec :: (RescaleCyc c a b, Fact m) => c m a -> c m b+{-# INLINE rescalePow #-}+{-# INLINE rescaleDec #-}+rescalePow = rescaleCyc Pow+rescaleDec = rescaleCyc Dec+

@@ -0,0 +1,447 @@+{-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts,+             NoImplicitPrelude, PolyKinds, RankNTypes, ScopedTypeVariables, +             TupleSections, TypeFamilies, TypeOperators, +             UndecidableInstances #-}++-- | Interface for cyclotomic tensors, and helper functions for tensor+-- indexing.++module Crypto.Lol.Cyclotomic.Tensor+( Tensor(..), CRTElt+-- * Top-level CRT functions+, hasCRTFuncs+, scalarCRT, mulGCRT, divGCRT, crt, crtInv, twaceCRT, embedCRT+-- * Tensor indexing+, Matrix, indexM, twCRTs+, zmsToIndexFact+, indexInfo+, extIndicesPowDec, extIndicesCRT, extIndicesCoeffs+, baseIndicesPow, baseIndicesDec, baseIndicesCRT+, digitRev+)+where++import Crypto.Lol.CRTrans+import Crypto.Lol.LatticePrelude as LP hiding (lift, (*>))+import Crypto.Lol.Types.FiniteField++import           Control.Applicative+import           Control.DeepSeq+import           Control.Monad.Random+import           Data.Constraint+import           Data.Singletons.Prelude hiding ((:-))+import           Data.Traversable+import           Data.Tuple           (swap)+import qualified Data.Vector          as V+import qualified Data.Vector.Unboxed  as U++-- | Synonym for constraints required for CRT-related functions.+type CRTElt t r = (ZeroTestable r, IntegralDomain r, CRTrans r, TElt t r)++-- | 'Tensor' encapsulates all the core linear transformations needed+-- for cyclotomic ring arithmetic.++-- | The type @t m r@ represents a cyclotomic coefficient tensor of+-- index @m@ over base ring @r@.  Most of the methods represent linear+-- transforms corresponding to operations in particular bases.+-- CRT-related methods are wrapped in 'Maybe' because they are+-- well-defined only when a CRT basis exists over the ring @r@ for+-- index @m@.++-- | The superclass constraints are for convenience, to ensure that we+-- can sample error tensors of 'Double's.++-- | __WARNING:__ as with all fixed-point arithmetic, the methods+-- in 'Tensor' may result in overflow (and thereby incorrect answers+-- and potential security flaws) if the input arguments are too close+-- to the bounds imposed by the base type.  The acceptable range of+-- inputs for each method is determined by the linear transform it+-- implements.++class (TElt t Double, TElt t (Complex Double))+      => Tensor (t :: Factored -> * -> *) where++  -- | Constraints needed by @t@ to hold type @r@.+  type TElt t r :: Constraint++  -- | Properties that hold for any index. Use with '\\'.+  entailIndexT :: Tagged (t m r)+                  (Fact m :- (Applicative (t m), Traversable (t m)))+  +  -- | Properties that hold for any (legal) fully-applied tensor. Use+  -- with '\\'.+  entailEqT :: Tagged (t m r)+               ((Eq r, Fact m, TElt t r) :- Eq (t m r))+  entailZTT :: Tagged (t m r)+               ((ZeroTestable r, Fact m, TElt t r) :- ZeroTestable (t m r))+  entailNFDataT :: Tagged (t m r)+                   ((NFData r, Fact m, TElt t r) :- NFData (t m r))+  entailRandomT :: Tagged (t m r)+                   ((Random r, Fact m, TElt t r) :- Random (t m r))+  entailShowT :: Tagged (t m r)+                 ((Show r, Fact m, TElt t r) :- (Show (t m r)))+  entailModuleT :: Tagged (GF fp d, t m fp)+                   ((GFCtx fp d, Fact m, TElt t fp) :- Module (GF fp d) (t m fp))++  -- | Convert a scalar to a tensor in the powerful basis.+  scalarPow :: (Additive r, Fact m, TElt t r) => r -> t m r++  {- CJP: suppressed to do annoyingly complicated algorithm++  -- | Convert a scalar to a tensor in the decoding basis.+  scalarDec :: (Additive r, Fact m, TElt t r) => r -> t m r+  -}++  -- | 'l' converts from decoding-basis representation to+  -- powerful-basis representation; 'lInv' is its inverse.+  l, lInv :: (Additive r, Fact m, TElt t r) => t m r -> t m r++  -- | Multiply by @g@ in the powerful/decoding basis+  mulGPow, mulGDec :: (Ring r, Fact m, TElt t r) => t m r -> t m r++  -- | Divide by @g@ in the powerful/decoding basis.  The 'Maybe'+  -- output indicates that the operation may fail, which happens+  -- exactly when the input is not divisible by @g@.+  divGPow, divGDec :: (ZeroTestable r, IntegralDomain r, Fact m, TElt t r)+                      => t m r -> Maybe (t m r)++  -- | A tuple of all the operations relating to the CRT basis, in a+  -- single 'Maybe' value for safety.  Clients should typically not+  -- use this method directly, but instead call the corresponding+  -- top-level functions: the elements of the tuple correpond to the+  -- functions 'scalarCRT', 'mulGCRT', 'divGCRT', 'crt', 'crtInv'.+  crtFuncs :: (Fact m, CRTElt t r) =>+              Maybe (    r -> t m r, -- scalarCRT+                     t m r -> t m r, -- mulGCRT+                     t m r -> t m r, -- divGCRT+                     t m r -> t m r, -- crt+                     t m r -> t m r) -- crtInv++  -- | Sample from the "tweaked" Gaussian error distribution @t*D@+  -- in the decoding basis, where @D@ has scaled variance @v@.+  tGaussianDec :: (OrdFloat q, Random q, TElt t q,+                   ToRational v, Fact m, MonadRandom rnd)+                  => v -> rnd (t m q)++  -- | Given the coefficient tensor of @e@ with respect to the+  -- decoding basis of @R@, yield the (scaled) squared norm of @g_m+  -- \cdot e@ under the canonical embedding, namely,+  --  @\hat{m}^{ -1 } \cdot || \sigma(g_m \cdot e) ||^2@ .+  gSqNormDec :: (Ring r, Fact m, TElt t r) => t m r -> r++  -- | The @twace@ linear transformation, which is the same in both the+  -- powerful and decoding bases.+  twacePowDec :: (Ring r, m `Divides` m', TElt t r) => t m' r -> t m r++  -- | The @embed@ linear transformations, for the powerful and+  -- decoding bases.+  embedPow, embedDec :: (Additive r, m `Divides` m', TElt t r)+                        => t m r -> t m' r++  -- | A tuple of all the extension-related operations involving the+  -- CRT bases, for safety.  Clients should typically not use this+  -- method directly, but instead call the corresponding top-level+  -- functions: the elements of the tuple correpond to the functions+  -- 'twaceCRT', 'embedCRT'.+  crtExtFuncs :: (m `Divides` m', CRTElt t r) =>+                 Maybe (t m' r -> t m  r, -- twaceCRT+                        t m  r -> t m' r) -- embedCRT++  -- | Map a tensor in the powerful\/decoding\/CRT basis, representing+  -- an @O_m'@ element, to a vector of tensors representing @O_m@+  -- elements in the same kind of basis.+  coeffs :: (Ring r, m `Divides` m', TElt t r) => t m' r -> [t m r]++  -- | The powerful extension basis w.r.t. the powerful basis.+  powBasisPow :: (Ring r, TElt t r, m `Divides` m') => Tagged m [t m' r]++  -- | A list of tensors representing the mod-@p@ CRT set of the+  -- extension.+  crtSetDec :: (m `Divides` m', PrimeField fp, Coprime (PToF (CharOf fp)) m',+                TElt t fp)+               => Tagged m [t m' fp]++  -- | Potentially optimized version of 'fmap' when the input and+  -- output element types satisfy 'TElt'.+  fmapT :: (Fact m, TElt t a, TElt t b) => (a -> b) -> t m a -> t m b+  -- | Potentially optimized monadic 'fmap'.+  fmapTM :: (Monad mon, Fact m, TElt t a, TElt t b)+             => (a -> mon b) -> t m a -> mon (t m b)++  -- | Potentially optimized 'zipWith'.+  zipWithT :: (Fact m, TElt t a, TElt t b, TElt t c)+              => (a -> b -> c) -> t m a -> t m b -> t m c++  -- | Unzip for types that satisfy 'TElt'.+  unzipTElt :: (Fact m, TElt t (a,b), TElt t a, TElt t b) +               => t m (a,b) -> (t m a, t m b)+  -- | Unzip for arbitrary types.+  unzipT :: (Fact m) => t m (a,b) -> (t m a, t m b)++-- | Convenience value indicating whether 'crtFuncs' exists.+hasCRTFuncs :: forall t m r . (Tensor t, Fact m, CRTElt t r)+               => TaggedT (t m r) Maybe ()+{-# INLINABLE hasCRTFuncs #-}+hasCRTFuncs = tagT $ do+  (_ :: r -> t m r,_,_,_,_) <- crtFuncs+  return ()++-- | Yield a tensor for a scalar in the CRT basis.  (This function is+-- simply an appropriate entry from 'crtFuncs'.)+scalarCRT :: (Tensor t, Fact m, CRTElt t r) => Maybe (r -> t m r)+{-# INLINABLE scalarCRT #-}+scalarCRT = (\(f,_,_,_,_) -> f) <$> crtFuncs+++mulGCRT, divGCRT, crt, crtInv ::+  (Tensor t, Fact m, CRTElt t r) => Maybe (t m r -> t m r)+{-# INLINABLE mulGCRT #-}+{-# INLINABLE divGCRT #-}+{-# INLINABLE crt #-}+{-# INLINABLE crtInv #-}++-- | Multiply by @g@ in the CRT basis. (This function is simply an+-- appropriate entry from 'crtFuncs'.)+mulGCRT = (\(_,f,_,_,_) -> f) <$> crtFuncs+-- | Divide by @g@ in the CRT basis.  (This function is simply an+-- appropriate entry from 'crtFuncs'.)+divGCRT = (\(_,_,f,_,_) -> f) <$> crtFuncs+-- | The CRT transform.  (This function is simply an appropriate entry+-- from 'crtFuncs'.)+crt = (\(_,_,_,f,_) -> f) <$> crtFuncs+-- | The inverse CRT transform.  (This function is simply an+-- appropriate entry from 'crtFuncs'.)+crtInv = (\(_,_,_,_,f) -> f) <$> crtFuncs++-- | The "tweaked trace" function for tensors in the CRT basis:+-- For cyclotomic indices m | m',+-- @Tw(x) = (mhat\/m\'hat) * Tr(g\'\/g * x)@.+-- (This function is simply an appropriate entry from 'crtExtFuncs'.)+twaceCRT :: forall t r m m' . (Tensor t, m `Divides` m', CRTElt t r)+            => Maybe (t m' r -> t m r)+twaceCRT = proxyT hasCRTFuncs (Proxy::Proxy (t m' r)) *>+           proxyT hasCRTFuncs (Proxy::Proxy (t m  r)) *>+           (fst <$> crtExtFuncs)++-- | Embed a tensor with index @m@ in the CRT basis to a tensor with+-- index @m'@ in the CRT basis.+-- (This function is simply an appropriate entry from 'crtExtFuncs'.)+embedCRT :: forall t r m m' . (Tensor t, m `Divides` m', CRTElt t r)+            => Maybe (t m r -> t m' r)+embedCRT = proxyT hasCRTFuncs (Proxy::Proxy (t m' r)) *>+           proxyT hasCRTFuncs (Proxy::Proxy (t m  r)) *>+           (snd <$> crtExtFuncs)++fMatrix :: forall m r mon . (Fact m, Monad mon, Ring r)+           => (forall pp . (PPow pp) => TaggedT pp mon (MatrixC r))+           -> TaggedT m mon (Matrix r)+fMatrix mat = tagT $ go $ sUnF (sing :: SFactored m)+  where go :: Sing (pplist :: [PrimePower]) -> mon (Matrix r)+        go spps = case spps of+          SNil -> return MNil+          (SCons spp rest) -> do+            rest' <- go rest+            mat' <- withWitnessT mat spp+            return $ MKron rest' mat'++-- deeply embedded DSL for Kronecker products of matrices++data MatrixC r = +  MC (Int -> Int -> r)           -- yields element i,j+  Int Int                        -- dims++-- | A Kronecker product of zero of more matrices over @r@.+data Matrix r = MNil | MKron (Matrix r) (MatrixC r)++-- | Extract the @(i,j)@ element of a 'Matrix'.+indexM :: Ring r => Matrix r -> Int -> Int -> r+indexM MNil 0 0 = LP.one+indexM (MKron m (MC mc r c)) i j =+  let (iq,ir) = i `divMod` r+      (jq,jr) = j `divMod` c+      in indexM m iq jq * mc ir jr++-- | The "tweaked" CRT^* matrix: @CRT^* . diag(sigma(g_m))@.+twCRTs :: (Fact m, CRTrans r) => TaggedT m Maybe (Matrix r)+twCRTs = fMatrix twCRTsPPow++-- | The "tweaked" CRT^* matrix (for prime powers): @CRT^* * diag(sigma(g_p))@.+twCRTsPPow :: (PPow pp, CRTrans r) => TaggedT pp Maybe (MatrixC r)+twCRTsPPow = do+  phi    <- pureT totientPPow+  iToZms <- pureT indexToZmsPPow+  jToPow <- pureT indexToPowPPow+  (wPow, _) <- crtInfoPPow+  gEmb <- gEmbPPow++  return $ MC (\j i -> let i' = iToZms i+                       in wPow (jToPow j * negate i') * gEmb i') phi phi++-- Reindexing functions++-- | Base-p digit reversal; input and output are in @[p^e]@.+digitRev :: PP -> Int -> Int+digitRev (_,0) 0 = 0+-- CJP: use accumulator to avoid multiple exponentiations?+digitRev (p,e) j +  | e >= 1 = let (q,r) = j `divMod` p+             in r * (p^(e-1)) + digitRev (p,e-1) q++indexToPowPPow, indexToZmsPPow :: PPow pp => Tagged pp (Int -> Int)+indexToPowPPow = indexToPow <$> ppPPow+indexToZmsPPow = indexToZms <$> ppPPow++-- | Convert a @Z_m^*@ index to a linear tensor index in @[m]@.+zmsToIndexFact :: Fact m => Tagged m (Int -> Int)+zmsToIndexFact = zmsToIndex <$> ppsFact++-- | For a prime power @p^e@, map a tensor index to the corresponding+-- power j in @[phi(p^e)]@, as in the powerful basis.+indexToPow :: PP -> Int -> Int+-- CJP: use accumulator to avoid multiple exponentiations?+indexToPow (p,e) j = let (jq,jr) = j `divMod` (p-1)+                     in p^(e-1)*jr + digitRev (p,e-1) jq++-- | For a prime power @p^e@, map a tensor index to the corresponding+-- element i in @Z_{p^e}^*@.+indexToZms :: PP -> Int -> Int+indexToZms (p,_) i = let (i1,i0) = i `divMod` (p-1)+                       in p*i1 + i0 + 1 ++-- | Convert a Z_m^* index to a linear tensor index.+zmsToIndex :: [PP] -> Int -> Int+zmsToIndex [] _ = 0+zmsToIndex (pp:rest) i = zmsToIndexPP pp (i `mod` valuePP pp)+                         + (totientPP pp) * zmsToIndex rest i++-- | Inverse of 'indexToZms'.+zmsToIndexPP :: PP -> Int -> Int+zmsToIndexPP (p,_) i = let (i1,i0) = i `divMod` p+                       in (p-1)*i1 + i0 - 1++-- Index correspondences for ring extensions++-- | Correspondences between the linear indexes into a basis of O_m',+-- and pair indices into (extension basis) \otimes (basis of O_m).+-- The work the same for Pow,Dec,CRT bases because all these bases+-- have that factorization.  The first argument is the list of+-- @(phi(m),phi(m'))@ pairs for the (merged) prime powers of @m@,@m'@.+toIndexPair :: [(Int,Int)] -> Int -> (Int,Int)+fromIndexPair :: [(Int,Int)] -> (Int,Int) -> Int++toIndexPair [] 0 = (0,0)+toIndexPair ((phi,phi'):rest) i' =+  let (i'q,i'r) = i' `divMod` phi'+      (i'rq,i'rr) = i'r `divMod` phi+      (i'q1,i'q0) = toIndexPair rest i'q+  in (i'rq + i'q1*(phi' `div` phi), i'rr + i'q0*phi)++fromIndexPair [] (0,0) = 0+fromIndexPair ((phi,phi'):rest) (i1,i0) =+  let (i0q,i0r) = i0 `divMod` phi+      (i1q,i1r) = i1 `divMod` (phi' `div` phi)+      i = fromIndexPair rest (i1q,i0q)+  in (i0r + i1r*phi) + i*phi'++-- | A collection of useful information for working with tensor+-- extensions.  The first component is a list of triples @(p,e,e')@+-- where @e@, @e'@ are respectively the exponents of prime @p@ in @m@,+-- @m'@.  The next two components are @phi(m)@ and @phi(m')@.  The+-- final component is a pair @(phi(p^e), phi(p^e'))@ for each triple+-- in the first component.+indexInfo :: forall m m' . (m `Divides` m')+             => Tagged '(m, m') ([(Int,Int,Int)], Int, Int, [(Int,Int)])+indexInfo = let pps = proxy ppsFact (Proxy::Proxy m)+                pps' = proxy ppsFact (Proxy::Proxy m')+                mpps = mergePPs pps pps'+                phi = totientPPs pps+                phi' = totientPPs pps'+                tots = totients mpps+            in tag (mpps, phi, phi', tots)++-- | A vector of @phi(m)@ entries, where the @i@th entry is the index+-- into the powerful\/decoding basis of @O_m'@ of the+-- @i@th entry of the powerful\/decoding basis of @O_m@.+extIndicesPowDec :: (m `Divides` m') => Tagged '(m, m') (U.Vector Int)+extIndicesPowDec = do+  (_, phi, _, tots) <- indexInfo+  return $ U.generate phi (fromIndexPair tots . (0,))++-- | A vector of @phi(m)@ blocks of @phi(m')\/phi(m)@ consecutive+-- entries. Each block contains all those indices into the CRT basis+-- of @O_m'@ that "lie above" the corresponding index into the CRT+-- basis of @O_m@.+extIndicesCRT :: forall m m' . (m `Divides` m')+                 => Tagged '(m, m') (U.Vector Int)+extIndicesCRT = do+  (_, phi, phi', tots) <- indexInfo+  return $ U.generate phi'+           (fromIndexPair tots . swap . (`divMod` (phi' `div` phi)))++baseWrapper :: forall m m' a . (m `Divides` m', U.Unbox a)+               => ([(Int,Int,Int)] -> Int -> a)+               -> Tagged '(m, m') (U.Vector a)+baseWrapper f = do+  (mpps, _, phi', _) <- indexInfo+  return $ U.generate phi' (f mpps)++-- | A lookup table for 'toIndexPair' applied to indices @[phi(m')]@.+baseIndicesPow :: forall m m' . (m `Divides` m')+                  => Tagged '(m, m') (U.Vector (Int,Int))+-- | A lookup table for 'baseIndexDec' applied to indices @[phi(m')]@.+baseIndicesDec :: forall m m' . (m `Divides` m')+                  => Tagged '(m, m') (U.Vector (Maybe (Int,Bool)))++-- | Same as 'baseIndicesPow', but only includes the second component+-- of each pair.+baseIndicesCRT :: forall m m' . (m `Divides` m')+                  => Tagged '(m, m') (U.Vector Int)++baseIndicesPow = baseWrapper (toIndexPair . totients)++-- this one is more complicated; requires the prime powers+baseIndicesDec = baseWrapper baseIndexDec++baseIndicesCRT =+  baseWrapper (\pps -> snd . toIndexPair (totients pps))+++-- | The @i0@th entry of the @i1@th vector is 'fromIndexPair' @(i1,i0)@.+extIndicesCoeffs :: forall m m' . (m `Divides` m')+                    => Tagged '(m, m') (V.Vector (U.Vector Int))+extIndicesCoeffs = do+  (_, phi, phi', tots) <- indexInfo+  return $ V.generate (phi' `div` phi)+           (\i1 -> U.generate phi (\i0 -> fromIndexPair tots (i1,i0)))++-- | Convenient reindexing functions++-- | Maps an index of the extension ring array to its corresponding+-- index in the base ring array (if it exists), with sign, under the+-- decoding basis.+baseIndexDec :: [(Int,Int,Int)] -> Int -> Maybe (Int, Bool)+baseIndexDec [] 0 = Just (0,False)+baseIndexDec ((p,e,e'):rest) i'+   = let (i'q, i'r) = i' `divMod` totientPP (p,e')+         phi = totientPP (p,e)+         curr+           | p>2 && e==0 && e' > 0 = case i'r of+               0 -> Just (0,False)+               1 -> Just (0,True)+               _ -> Nothing+           | otherwise = if i'r < phi then Just (i'r,False) else Nothing+     in do+       (i,b) <- curr+       (j,b') <- baseIndexDec rest i'q+       return (i + phi*j, b /= b')++-- the first list of pps must "divide" the other.  result is a list of+-- all (prime, min e, max e).+mergePPs :: [PP] -> [PP] -> [(Int,Int,Int)]+mergePPs [] pps = LP.map (\(p,e) -> (p,0,e)) pps+mergePPs allpps@((p,e):pps) ((p',e'):pps')+  | p == p' && e <= e' = (p,  e, e') : mergePPs pps pps'+  | p > p'  = (p', 0, e') : mergePPs allpps pps'++totients :: [(Int, Int, Int)] -> [(Int,Int)]+totients = LP.map (\(p,e,e') -> (totientPP (p,e), totientPP (p,e')))

@@ -0,0 +1,474 @@+{-# LANGUAGE ConstraintKinds, DataKinds, GADTs,+             FlexibleContexts, FlexibleInstances, TypeOperators, PolyKinds,+             GeneralizedNewtypeDeriving, InstanceSigs, RoleAnnotations,+             MultiParamTypeClasses, NoImplicitPrelude, StandaloneDeriving,+             ScopedTypeVariables, TupleSections, TypeFamilies, RankNTypes,+             TypeSynonymInstances, UndecidableInstances,+             RebindableSyntax #-}++-- | Wrapper for a C implementation of the 'Tensor' interface.++module Crypto.Lol.Cyclotomic.Tensor.CTensor+( 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 (liftM)+import Control.Monad.Identity (Identity(..), runIdentity)+import Control.Monad.Random+import Control.Monad.Trans (lift)++import Data.Coerce+import Data.Constraint  hiding ((***))+import Data.Foldable as F+import Data.Int+import Data.Maybe+import Data.Traversable as T+import Data.Vector.Generic           as V (zip, unzip, fromList, toList)+import Data.Vector.Storable          as SV (Vector, (!), replicate, replicateM, thaw, convert, foldl',+                                            unsafeSlice, mapM, fromList, toList,+                                            generate, foldl1',+                                            unsafeWith, zipWith, map, length, unsafeFreeze, thaw)+import Data.Vector.Storable.Internal (getPtr)+import Data.Vector.Storable.Mutable  as SM hiding (replicate)++import Foreign.Marshal.Utils (with)+import Foreign.Ptr+import Foreign.Storable        (Storable (..))+import Test.QuickCheck         hiding (generate)++import Crypto.Lol.CRTrans+import Crypto.Lol.Cyclotomic.Tensor+import Crypto.Lol.Cyclotomic.Tensor.CTensor.Backend+import Crypto.Lol.Cyclotomic.Tensor.CTensor.Extension+import Crypto.Lol.GaussRandom+import Crypto.Lol.LatticePrelude as LP hiding (replicate, unzip, zip, lift)+import Crypto.Lol.Reflects+import Crypto.Lol.Types.FiniteField+import Crypto.Lol.Types.IZipVector+import Crypto.Lol.Types.ZqBasic++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++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++deriving instance Show r => Show (CT m r)++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 r) => (CT' l r -> CT' m r) -> (CT l r -> CT m r)+wrap f (CT v) = CT $ f v+wrap f (ZV v) = CT $ f $ zvToCT' v++wrapM :: (Storable r, Monad mon) => (CT' l r -> mon (CT' m r))+         -> (CT l r -> mon (CT m r))+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 C function for zipWith (*)++instance (Additive r, Storable r, Fact m, Dispatch r)+  => Additive.C (CT m r) where+  (CT a@(CT' _)) + (CT b@(CT' _)) = CT $ (untag $ cZipDispatch dadd) 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, Fact m)+         => 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 $ untag $ basicDispatch dl+  lInv = wrap $ untag $ basicDispatch dlinv++  mulGPow = wrap mulGPow'+  mulGDec = wrap $ untag $ basicDispatch dmulgdec++  divGPow = wrapM divGPow'+  -- we divide by p in the C code (for divGDec only(?)), do NOT call checkDiv!+  divGDec = wrapM $ Just . untag (basicDispatch dginvdec)++  crtFuncs = (,,,,) <$>+    Just (CT . repl) <*>+    (wrap <$> untag (cZipDispatch dmul) <$> untagT gCoeffsCRT) <*>+    (wrap <$> untag (cZipDispatch dmul) <$> untagT gInvCoeffsCRT) <*>+    (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 (CT v) = CT $ coerce (SV.map f) v+  fmapT f v@(ZV _) = fmapT f $ toCT v++  fmapTM f (CT (CT' v)) = (CT . CT') <$> SV.mapM f v+  fmapTM f v@(ZV _) = fmapTM f $ toCT v++  zipWithT f (CT (CT' v1)) (CT (CT' v2)) = CT $ CT' $ SV.zipWith f v1 v2+  zipWithT f v1 v2 = zipWithT f (toCT v1) (toCT v2)++  unzipTElt (CT (CT' v)) = (CT . CT') *** (CT . CT') $ unzip v+  unzipTElt v = unzipTElt $ toCT v++  unzipT v@(CT _) = unzipT $ toZV v+  unzipT (ZV v) = ZV *** ZV $ unzipIZV v++  {-# 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 #-}+  {-# INLINABLE crtExtFuncs #-}+  {-# INLINABLE coeffs #-}+  {-# INLINABLE powBasisPow #-}+  {-# INLINABLE crtSetDec #-}+  {-# INLINABLE fmapT #-}+  {-# INLINABLE fmapTM #-}+  {-# INLINABLE zipWithT #-}+  {-# INLINABLE unzipTElt #-}+  {-# INLINABLE 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 :: (Fact m, Fact m') +  => 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 :: +  (Fact m, Fact m')+  => Tagged '(m,m') [Vector r] -> Tagged m [CT' m' r]+coerceBasis = coerce++mulGPow' :: (TElt CT r, Fact m, Additive r) => CT' m r -> CT' m r+mulGPow' = untag $ basicDispatch dmulgpow++divGPow' :: forall m r . (TElt CT r, Fact m, IntegralDomain r, ZeroTestable r) => CT' m r -> Maybe (CT' m r)+divGPow' = untag $ checkDiv $ basicDispatch dginvpow++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 :: forall m r .+     (Storable r, Fact m, Additive r)+      => (Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ())+         -> Tagged m (CT' m r -> CT' m r)+basicDispatch f = return $ unsafePerformIO . withBasicArgs f++gSqNormDec' :: forall m r .+     (Storable r, Fact m, Additive r, Dispatch r)+      => Tagged m (CT' m r -> r)+gSqNormDec' = return $ (!0) . unCT . unsafePerformIO . withBasicArgs dnorm++ctCRT :: forall m r . (Storable r, CRTrans r, Dispatch r,+          Fact m)+         => TaggedT m Maybe (CT' m r -> CT' m r)+ctCRT = do -- in TaggedT m Maybe+  ru' <- ru+  return $ \x -> unsafePerformIO $ +    withPtrArray ru' (flip withBasicArgs x . dcrt)++-- CTensor CRT^(-1) functions take inverse rus+ctCRTInv :: (Storable r, CRTrans r, Dispatch r,+          Fact m)+         => TaggedT m Maybe (CT' m r -> CT' m r)+ctCRTInv = do -- in Maybe+  mhatInv <- snd <$> crtInfoFact+  ruinv' <- ruInv+  return $ \x -> unsafePerformIO $ +    withPtrArray ruinv' (\ruptr -> with mhatInv (flip withBasicArgs x . dcrtinv ruptr))++checkDiv :: forall m r . +  (IntegralDomain r, Storable r, ZeroTestable r, +   Fact m)+    => Tagged m (CT' m r -> CT' m r) -> Tagged m (CT' m r -> Maybe (CT' m r))+checkDiv f = do+  f' <- f+  oddRad' <- fromIntegral <$> oddRadicalFact+  return $ \x -> +    let (CT' y) = f' x+    in CT' <$> SV.mapM (`divIfDivis` oddRad') y++divIfDivis :: (IntegralDomain r, ZeroTestable r) => r -> r -> Maybe r+divIfDivis num den = let (q,r) = num `divMod` den+                     in if isZero r then Just q else Nothing++cZipDispatch :: (Storable r, Fact m, Additive r)+  => (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)+  ruinv' <- mapTaggedT (return . fromMaybe (error "complexGaussianRoots")) ruInv+  totm <- pureT totientFact+  m <- pureT valueFact+  rad <- pureT radicalFact+  yin <- lift $ realGaussians (var * fromIntegral (m `div` rad)) totm+  return $ unsafePerformIO $ +    withPtrArray ruinv' (\ruptr -> withBasicArgs (dgaussdec ruptr) (CT' yin))++instance (Arbitrary r, Fact m, Storable r) => Arbitrary (CT' m r) where+  arbitrary = replM arbitrary+  shrink = shrinkNothing++instance (Storable r, Arbitrary (CT' m r)) => Arbitrary (CT m r) where+  arbitrary = CT <$> arbitrary++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++--{-# INLINE scalarPow' #-}+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 :: forall r m . +   (CRTrans r, Fact m, Storable r)+   => TaggedT m Maybe [Vector r]+--{-# INLINE ru #-}+ru = do+  mval <- pureT valueFact+  wPow <- fst <$> crtInfoFact+  LP.map+    (\(p,e) -> do+        let pp = p^e+            pow = mval `div` pp+        generate pp (wPow . (*pow))) <$>+      pureT ppsFact++--{-# INLINE ruInv #-}+ruInv = do+  mval <- pureT valueFact+  wPow <- fst <$> crtInfoFact+  LP.map+    (\(p,e) -> do+        let pp = p^e+            pow = mval `div` pp+        generate pp (\i -> wPow $ -i*pow)) <$>+      pureT ppsFact++gCoeffsCRT, gInvCoeffsCRT :: (TElt CT r, CRTrans r, Fact m, ZeroTestable r, IntegralDomain r)+  => TaggedT m Maybe (CT' m r)+gCoeffsCRT = ctCRT <*> return (mulGPow' $ scalarPow' LP.one)+-- It's necessary to call 'fromJust' here: otherwise +-- sequencing functions in 'crtFuncs' relies on 'divGPow' having an+-- implementation in C, which is not true for all types which have a C+-- implementation of, e.g. 'crt'. In particular, 'Complex Double' has C support+-- for 'crt', but not for 'divGPow'.+-- This really breaks the contract of Tensor, so it's probably a bad idea.+--   Someone can get the "crt" and can even pull the function "divGCRT" from Tensor,+--   but it will fail when they try to apply it.+-- As an implementation note if I ever do fix this: the division by rad(m) can be+-- tricky for Double/Complex Doubles, so be careful! This is why we have a custom+-- Complex wrapper around NP.Complex.+gInvCoeffsCRT = ($ fromJust $ divGPow' $ scalarPow' LP.one) <$> ctCRT++-- we can't put this in Extension with the rest of the twace/embed fucntions because it needs access to +-- the C backend+twaceCRT' :: forall m m' r .+             (TElt CT r, CRTrans r, m `Divides` m', ZeroTestable r, IntegralDomain r)+             => TaggedT '(m, m') Maybe (Vector r -> Vector r)+twaceCRT' = tagT $ do -- Maybe monad+  (CT' g') <- proxyT gCoeffsCRT (Proxy::Proxy m')+  (CT' gInv) <- proxyT gInvCoeffsCRT (Proxy::Proxy m)+  embed <- proxyT embedCRT' (Proxy::Proxy '(m,m'))+  indices <- pure $ proxy extIndicesCRT (Proxy::Proxy '(m,m'))+  (_, m'hatinv) <- proxyT crtInfoFact (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

@@ -0,0 +1,304 @@+{-# LANGUAGE ConstraintKinds, FlexibleContexts, FlexibleInstances, MultiParamTypeClasses, +             PolyKinds, ScopedTypeVariables, TypeFamilies, UndecidableInstances #-}++module Crypto.Lol.Cyclotomic.Tensor.CTensor.Backend+(Dispatch+,dcrt+,dcrtinv+,dgaussdec+,dl+,dlinv+,dnorm+,dmulgpow+,dmulgdec+,dginvpow+,dginvdec+,dadd+,dmul+,marshalFactors+,CPP+,withArray+,withPtrArray+) where++import Control.Applicative+import Control.Monad++import Crypto.Lol.LatticePrelude as LP (Complex, Proxy(..), proxy, (++), map, mapM_, PP, Tagged, tag)+import Crypto.Lol.Reflects+import Crypto.Lol.Types.ZqBasic++import Data.Int+import Data.Vector.Storable          as SV (Vector, (!), replicate, replicateM, thaw, convert, foldl',+                                            unsafeToForeignPtr0, unsafeSlice, mapM, fromList,+                                            generate, foldl1',+                                            unsafeWith, zipWith, map, length, unsafeFreeze, thaw)+import Data.Vector.Storable.Internal (getPtr)++import           Foreign.ForeignPtr (touchForeignPtr)+import           Foreign.Marshal.Array (withArray)+import           Foreign.Marshal.Utils (with)+import           Foreign.Ptr (plusPtr, castPtr, Ptr)+import           Foreign.Storable        (Storable (..))+import qualified Foreign.Storable.Record as Store++marshalFactors :: [PP] -> Vector CPP+marshalFactors = SV.fromList . LP.map (\(p,e) -> CPP (fromIntegral p) (fromIntegral e))++-- http://stackoverflow.com/questions/6517387/vector-vector-foo-ptr-ptr-foo-io-a-io-a+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++data CPP = CPP {p' :: !Int32, e' :: !Int16}+-- stolen from http://hackage.haskell.org/packages/archive/numeric-prelude/0.4.0.3/doc/html/src/Number-Complex.html#T+-- the NumericPrelude Storable instance for complex numbers+instance Storable CPP where+   sizeOf    = Store.sizeOf store+   alignment = Store.alignment store+   peek      = Store.peek store+   poke      = Store.poke store++store :: Store.Dictionary CPP+store = Store.run $+   liftA2 CPP+      (Store.element p')+      (Store.element e')++instance Show CPP where+    show (CPP p e) = "(" LP.++ show p LP.++ "," LP.++ show e LP.++ ")"++instance (Storable a, Storable b,+          CTypeOf a ~ CTypeOf b) +  -- enforces right associativity and that each type of +  -- the tuple has the same C repr, so using an array repr is safe+  => 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++type family CTypeOf x where+  CTypeOf (a,b) = CTypeOf a+  CTypeOf (ZqBasic (q :: k) Int64) = ZqB64D+  CTypeOf Double = DoubleD+  CTypeOf Int64 = Int64D+  CTypeOf (Complex Double) = ComplexD++-- 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 (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))++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+  dcrt      :: Ptr (Ptr r) ->           Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ()+  dcrtinv   :: Ptr (Ptr r) -> Ptr r ->  Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ()+  dgaussdec :: Ptr (Ptr (Complex r)) -> Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ()++  dl        :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ()+  dlinv     :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ()+  dnorm     :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ()+  dmulgpow  :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ()+  dmulgdec  :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ()+  dginvpow  :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ()+  dginvdec  :: Ptr r -> Int64 -> Ptr CPP -> Int16 -> IO ()+  +  dadd :: Ptr r -> Ptr r -> Int64 -> IO ()+  dmul :: Ptr r -> Ptr r -> Int64 -> IO ()++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)+  dadd aout bout totm = +    let qs = proxy getModuli (Proxy::Proxy r)+        numPairs = proxy numComponents (Proxy::Proxy r)+    in with qs $ \qsptr ->+        addRq numPairs (castPtr aout) (castPtr bout) totm (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"++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 = error "cannot call CT mulGDec on type Complex Double"+  dginvpow pout = +    tensorGInvPowC (proxy numComponents (Proxy::Proxy r)) (castPtr pout)+  dginvdec = error "cannot call CT divGDec on type Complex Double"+  dadd aout bout = +    addC (proxy numComponents (Proxy::Proxy r)) (castPtr aout) (castPtr bout)+  dmul aout bout = +    mulC (proxy numComponents (Proxy::Proxy r)) (castPtr aout) (castPtr bout)+  dgaussdec = error "cannot call CT gaussianDec on type Comple Double"++instance (Tuple r, CTypeOf r ~ DoubleD) => Dispatch' DoubleD r where+  dcrt = error "cannot call CT Crt on type Double"+  dcrtinv = error "cannot call CT CrtInv on type Double"+  dl pout = +    tensorLDouble (proxy numComponents (Proxy::Proxy r)) (castPtr pout)+  dlinv pout = +    tensorLInvDouble (proxy numComponents (Proxy::Proxy r)) (castPtr pout)+  dnorm = error "cannot call CT normSq on type Double"+  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"+  dadd aout bout = +    addD (proxy numComponents (Proxy::Proxy r)) (castPtr aout) (castPtr bout)+  dmul = error "cannot call CT (*) on type Double"+  dgaussdec ruptr pout totm pfac numFacts = +    tensorGaussianDec (proxy numComponents (Proxy::Proxy r)) (castPtr pout) totm pfac numFacts (castPtr ruptr)++instance (Tuple r, CTypeOf r ~ Int64D) => Dispatch' Int64D r where+  dcrt = error "cannot call CT Crt on type Int64"+  dcrtinv = error "cannot call CT CrtInv on type Int64"+  dl pout = +    tensorLR (proxy numComponents (Proxy::Proxy r)) (castPtr pout)+  dlinv pout = +    tensorLInvR (proxy numComponents (Proxy::Proxy r)) (castPtr pout)+  dnorm pout = +    tensorNormSqR (proxy numComponents (Proxy::Proxy r)) (castPtr pout)+  dmulgpow pout = +    tensorGPowR (proxy numComponents (Proxy::Proxy r)) (castPtr pout)+  dmulgdec pout = +    tensorGDecR (proxy numComponents (Proxy::Proxy r)) (castPtr pout)+  dginvpow pout = +    tensorGInvPowR (proxy numComponents (Proxy::Proxy r)) (castPtr pout)+  dginvdec pout = +    tensorGInvDecR (proxy numComponents (Proxy::Proxy r)) (castPtr pout)+  dadd aout bout = +    addR (proxy numComponents (Proxy::Proxy r)) (castPtr aout) (castPtr bout)+  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 "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 "tensorGInvPowR" tensorGInvPowR ::   Int16 -> Ptr Int64 -> Int64 -> Ptr CPP -> Int16          -> IO ()+foreign import ccall unsafe "tensorGInvPowRq" tensorGInvPowRq :: Int16 -> Ptr (ZqBasic q Int64) -> Int64 -> Ptr CPP -> Int16 -> Ptr Int64 -> IO ()+foreign import ccall unsafe "tensorGInvPowC" tensorGInvPowC ::   Int16 -> Ptr (Complex Double) -> Int64 -> Ptr CPP -> Int16          -> IO ()+foreign import ccall unsafe "tensorGInvDecR" tensorGInvDecR ::   Int16 -> Ptr Int64 -> Int64 -> Ptr CPP -> Int16          -> IO ()+foreign import ccall unsafe "tensorGInvDecRq" tensorGInvDecRq :: Int16 -> Ptr (ZqBasic q Int64) -> Int64 -> Ptr CPP -> Int16 -> Ptr Int64 -> IO ()++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 ()++foreign import ccall unsafe "addRq" addRq :: Int16 -> Ptr (ZqBasic q Int64) -> Ptr (ZqBasic q Int64) -> Int64 -> Ptr Int64 -> IO ()+foreign import ccall unsafe "addR" addR :: Int16 -> Ptr Int64 -> Ptr Int64 -> Int64 -> IO ()+foreign import ccall unsafe "addC" addC :: Int16 -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int64 -> IO ()+foreign import ccall unsafe "addD" addD :: Int16 -> Ptr Double -> Ptr Double -> Int64 -> IO ()

@@ -0,0 +1,121 @@+{-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts, NoImplicitPrelude, +             PolyKinds, RebindableSyntax, ScopedTypeVariables,+             TypeFamilies, TypeOperators #-}++-- | CT-specific functions for embedding/twacing in various bases++module Crypto.Lol.Cyclotomic.Tensor.CTensor.Extension+( embedPow', embedDec', embedCRT'+, twacePowDec' -- , twaceCRT'+, coeffs', powBasisPow'+, crtSetDec'+, backpermute'+) where++import Crypto.Lol.CRTrans+import Crypto.Lol.LatticePrelude as LP hiding (null, lift)+import Crypto.Lol.Cyclotomic.Tensor as T+import Crypto.Lol.Types.FiniteField+import Crypto.Lol.Types.ZmStar+import Crypto.Lol.Reflects++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.Generic as G (generate, Vector, (!), length)+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Storable as SV+++-- | /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 <a,b,c,d> <0,3,2,3,1,0> = <a,d,c,d,b,a>+backpermute' :: (Vector v a)+             => U.Vector Int -- ^ @is@ index vector (of length @n@)+             -> v a   -- ^ @xs@ value vector+             -> v a+--{-# INLINE backpermute' #-}+backpermute' is v = generate (G.length is) (\i -> v ! (is ! i))++embedPow', embedDec' :: (Additive r, Vector v r, m `Divides` m')+                     => Tagged '(m, m') (v r -> v r)+-- | 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 ! 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 m m' v r . (CRTrans r, Vector v r, m `Divides` m')+          => TaggedT '(m, m') Maybe (v r -> v r)+embedCRT' = +  (lift (proxyT crtInfoFact (Proxy::Proxy m') :: Maybe (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' :: (Vector v r, m `Divides` m')+        => Tagged '(m, m') (v r -> [v 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 v . (Vector v r, m `Divides` m')+             => Tagged '(m, m') (v r -> v r)+twacePowDec' = backpermute' <$> extIndicesPowDec+++-- EAC: twaceCRT is defined in CTensor because it needs access to C-backend functions+++-- | 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' :: Matrix (GF fp d)+            = fromMaybe (error "internal error: crtSetDec': twCRTs") $ proxyT twCRTs m'p+          zmsToIdx = proxy T.zmsToIndexFact m'p+          elt j i = indexM 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'+                                      (sum $ LP.map (elt j) is))) cosets

@@ -0,0 +1,175 @@+#include "tensorTypes.h"+#ifdef CINTRIN+#include <immintrin.h>+#endif++#ifdef STATS+int mulCtr = 0;+struct timespec mulTime = {0,0};++int addCtr = 0;+struct timespec addTime = {0,0};+#endif++//a = zipWith (*) a b+void mulRq (hShort_t tupSize, hInt_t* a, hInt_t* b, hDim_t totm, hInt_t* qs) {+#ifdef STATS+    mulCtr++;+    struct timespec s1,t1;+    clock_gettime(CLOCK_REALTIME, &s1);+#endif+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {+        hInt_t q = qs[tupIdx];+        for(int i = 0; i < totm; i++) {+            a[i*tupSize+tupIdx] = (a[i*tupSize+tupIdx]*b[i*tupSize+tupIdx])%q;+        }+    }+#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    mulTime = tsAdd(mulTime, tsSubtract(t1,s1));+#endif+}+/*+void mulMq (hShort_t tupSize, hInt_t* a, const hInt_t* b, const hDim_t totm, const hByte_t logr, const hInt_t k, const hInt_t q) {+#ifdef STATS+    mulCtr++;+    struct timespec s1,t1;+    clock_gettime(CLOCK_REALTIME, &s1);+#endif+    hInt_t mask = (1<<logr)-1; // R-1++    for(int i = 0; i < totm; i++) {+        hInt_t x = a[i]*b[i];+        hInt_t s = k*(x & mask);+        hInt_t m = s & mask;+        a[i] = (x+m*q)>>logr;+    }++#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    mulTime = tsAdd(mulTime, tsSubtract(t1,s1));+#endif+}+*/+void mulC (hShort_t tupSize, complex_t* a, complex_t* b, hDim_t totm) {+#ifdef STATS+    mulCtr++;+    struct timespec s1,t1;+    clock_gettime(CLOCK_REALTIME, &s1);+#endif+    for(int i = 0; i < totm*tupSize; i++)+    {+        CMPLX_IMUL(a[i],b[i]);+    }+#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    mulTime = tsAdd(mulTime, tsSubtract(t1,s1));+#endif+}++//a = zipWith (+) a b+void addRq (hShort_t tupSize, hInt_t* a, const hInt_t* b, const hDim_t totm, const hInt_t* qs) {+#ifdef STATS+    addCtr++;+    struct timespec s1,t1;+    clock_gettime(CLOCK_REALTIME, &s1);+#endif+/*+need to check this code more carefully when tupSize > 1+#ifdef CINTRIN++    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {+        hInt_t q = qs[tupIdx];+        __m128i qss = _mm_set1_epi64x(q);+        for(int i = 0; i < totm; i+=2) {+            __m128i xs = _mm_load_si128((const __m128i*)(a+(i*tupSize+tupIdx)));+            __m128i ys = _mm_load_si128((const __m128i*)(b+(i*tupSize+tupIdx)));+            __m128i zs = _mm_add_epi64(xs,ys);+            zs = _mm_rem_epi64(zs,qss);+            _mm_store_si128((__m128i*)(a+(i*tupSize+tupIdx)),zs);+        }+    }+#else+*/+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {+        hInt_t q = qs[tupIdx];+        for(int i = 0; i < totm; i++) {+            hInt_t temp = a[i*tupSize+tupIdx]+b[i*tupSize+tupIdx];+            if (temp >= q) a[i*tupSize+tupIdx]=temp-q;+            else a[i*tupSize+tupIdx] = temp;+        }+    }+//#endif+#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    addTime = tsAdd(addTime, tsSubtract(t1,s1));+#endif+}+/*+void addMq (hInt_t* a, const hInt_t* b, const hDim_t totm, const hByte_t logr, const hInt_t k, const hInt_t q) {+#ifdef STATS+    addCtr++;+    struct timespec s1,t1;+    clock_gettime(CLOCK_REALTIME, &s1);+#endif+    hInt_t twoq = q<<1;+    for(int i = 0; i < totm; i++) {+        hInt_t temp = (a[i]+b[i]);+        if (temp >= twoq) a[i]=temp-twoq;+        else a[i] = temp;+    }++#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    addTime = tsAdd(addTime, tsSubtract(t1,s1));+#endif+}+*/+//a = zipWith (+) a b+void addR (hShort_t tupSize, hInt_t* a, hInt_t* b, hDim_t totm) {+#ifdef STATS+    addCtr++;+    struct timespec s1,t1;+    clock_gettime(CLOCK_REALTIME, &s1);+#endif+    for(int i = 0; i < totm*tupSize; i++)    {+        a[i] += b[i];+    }+#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    addTime = tsAdd(addTime, tsSubtract(t1,s1));+#endif+}++void addC (hShort_t tupSize, complex_t* a, complex_t* b, hDim_t totm) {+#ifdef STATS+    addCtr++;+    struct timespec s1,t1;+    clock_gettime(CLOCK_REALTIME, &s1);+#endif+    for(int i = 0; i < totm*tupSize; i++)+    {+        CMPLX_IADD(a[i],b[i]);+    }+#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    addTime = tsAdd(addTime, tsSubtract(t1,s1));+#endif+}++void addD (hShort_t tupSize, double* a, double* b, hDim_t totm) {+#ifdef STATS+    addCtr++;+    struct timespec s1,t1;+    clock_gettime(CLOCK_REALTIME, &s1);+#endif+    for(int i = 0; i < totm*tupSize; i++)+    {+        a[i]+=b[i];+    }+#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    addTime = tsAdd(addTime, tsSubtract(t1,s1));+#endif+}+

@@ -0,0 +1,1386 @@+#include "tensorTypes.h"+#include <time.h>+#include <stdlib.h>++// there should be a special cases that do NOT require temp space to be allocated for all primes *smaller* than DFTP_GENERIC_SIZE+#define DFTP_GENERIC_SIZE 11++#ifdef STATS+int crtRqCtr = 0;+int crtInvRqCtr = 0;+int crtCCtr = 0;+int crtInvCCtr = 0;++struct timespec crttime1 = {0,0};+struct timespec crttime2 = {0,0};+struct timespec crttime3 = {0,0};+struct timespec crttime4 = {0,0};++struct timespec crtInvRqTime = {0,0};+struct timespec crtCTime = {0,0};+struct timespec crtInvCTime = {0,0};+#endif++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;+}++void crtTwiddleRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, +                   PrimeExponent pe, hInt_t* ru, hInt_t q)+{+    hDim_t p = pe.prime;+    hShort_t e = pe.exponent;+    +#ifdef DEBUG_MODE+    ASSERT(e != 0);+#endif+    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;+            hInt_t 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] = (y[idx]*twid) % q;+                }+            }+        }+    }+    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;+                hInt_t 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] = (y[idx]*twid) % q;+                    }+                }+            }+        }+    }+}++// dim is power of p+void dftptwidRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, +                 PrimeExponent pe, hDim_t dim, hDim_t rustride, hInt_t* ru, hInt_t q)+{+    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);+            hInt_t 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] = (y[idx]*twid) % q;+                }+            }+        }+    }+    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);+                hInt_t 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] = (y[idx]*twid) % q;+                    }+                }+            }+        }+    }+}++//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+void dftpRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, +             hDim_t p, hDim_t rustride, hInt_t* ru, hInt_t* tempSpace, hInt_t q)+{+    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;+                hInt_t u = y[tensorOffset*tupSize];+                hInt_t t = y[(tensorOffset+rts)*tupSize];+                y[tensorOffset*tupSize] = (u + t) % q;+                y[(tensorOffset+rts)*tupSize] = (u - t) % q;+            }+        }+    }+    else if(p == 3)+    {+        hInt_t ru1 = ru[rustride*tupSize];+        hInt_t 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;+                hInt_t 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]          = (y1 + y2 + y3) % q;+                y[(tensorOffset+rts)*tupSize]      = (y1 + ((ru1*y2) % q) + ((ru2*y3) % q)) % q;+                y[(tensorOffset+(rts<<1))*tupSize] = (y1 + ((ru2*y2) % q) + ((ru1*y3) % q)) % q;+            }   +        }++    }+    else if(p == 5)+    {+        hDim_t temp1 = rts*5;+        hInt_t ru1 = ru[rustride*tupSize];+        hInt_t ru2 = ru[(rustride<<1)*tupSize];+        hInt_t ru3 = ru[(rustride*3)*tupSize];+        hInt_t 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;+                hInt_t 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]          = (y1 + y2 + y3 + y4 + y5) % q;+                y[(tensorOffset+rts)*tupSize]      = (y1 + ((ru1*y2) % q) + ((ru2*y3) % q) + ((ru3*y4) % q) + ((ru4*y5) % q)) % q;+                y[(tensorOffset+(rts<<1))*tupSize] = (y1 + ((ru2*y2) % q) + ((ru4*y3) % q) + ((ru1*y4) % q) + ((ru3*y5) % q)) % q;+                y[(tensorOffset+rts*3)*tupSize]    = (y1 + ((ru3*y2) % q) + ((ru1*y3) % q) + ((ru4*y4) % q) + ((ru2*y5) % q)) % q;+                y[(tensorOffset+(rts<<2))*tupSize] = (y1 + ((ru4*y2) % q) + ((ru3*y3) % q) + ((ru2*y4) % q) + ((ru1*y5) % q)) % q;+            }+        }+    }+    else if(p == 7)+    {+        hDim_t temp1 = rts*7;+        hInt_t ru1 = ru[rustride*tupSize];+        hInt_t ru2 = ru[(rustride<<1)*tupSize];+        hInt_t ru3 = ru[(rustride*3)*tupSize];+        hInt_t ru4 = ru[(rustride<<2)*tupSize];+        hInt_t ru5 = ru[(rustride*5)*tupSize];+        hInt_t 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;+                hInt_t 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]          = (y1 +     y2 +     y3 +     y4 +     y5 +     y6 +     y7) % q;+                y[(tensorOffset+rts)*tupSize]      = (y1 + ((ru1*y2) % q) + ((ru2*y3) % q) + ((ru3*y4) % q) + ((ru4*y5) % q) + ((ru5*y6) % q) + ((ru6*y7) % q)) % q;+                y[(tensorOffset+(rts<<1))*tupSize] = (y1 + ((ru2*y2) % q) + ((ru4*y3) % q) + ((ru6*y4) % q) + ((ru1*y5) % q) + ((ru3*y6) % q) + ((ru5*y7) % q)) % q;+                y[(tensorOffset+rts*3)*tupSize]    = (y1 + ((ru3*y2) % q) + ((ru6*y3) % q) + ((ru2*y4) % q) + ((ru5*y5) % q) + ((ru1*y6) % q) + ((ru4*y7) % q)) % q;+                y[(tensorOffset+(rts<<2))*tupSize] = (y1 + ((ru4*y2) % q) + ((ru1*y3) % q) + ((ru5*y4) % q) + ((ru2*y5) % q) + ((ru6*y6) % q) + ((ru3*y7) % q)) % q;+                y[(tensorOffset+rts*5)*tupSize]    = (y1 + ((ru5*y2) % q) + ((ru3*y3) % q) + ((ru1*y4) % q) + ((ru6*y5) % q) + ((ru4*y6) % q) + ((ru2*y7) % q)) % q;+                y[(tensorOffset+rts*6)*tupSize]    = (y1 + ((ru6*y2) % q) + ((ru5*y3) % q) + ((ru4*y4) % q) + ((ru3*y5) % q) + ((ru2*y6) % q) + ((ru1*y7) % q)) % q;+            }   +        }+    }+    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++)+                {+                    hInt_t acc = 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++)+                    {+                        acc += ((y[(tensorOffset+col*rts)*tupSize]*ru[((col*row) % p)*rustride*tupSize])%q);+                    }+                    tempSpace[row] = acc % q;+                }+                +                for(hDim_t row = 0; row < p; row++)+                {+                    y[(tensorOffset+rts*row)*tupSize] = tempSpace[row];+                }+            }+        }+    }+}++void crtpRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, +             hDim_t p, hDim_t rustride, hInt_t* ru, hInt_t q)+{+    hDim_t tensorOffset;+    if(p == 2)+    {+        return;+    }+    else if(p == 3)+    {+        hDim_t temp1 = rts*2;+        hInt_t ru1 = ru[rustride*tupSize];+        hInt_t 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;+                hInt_t y1, y2;+                y1 = y[tensorOffset*tupSize];+                y2 = y[(tensorOffset+rts)*tupSize];+                y[tensorOffset*tupSize]     = (y1 + ((ru1*y2)%q)) % q;+                y[(tensorOffset+rts)*tupSize] = (y1 + ((ru2*y2)%q)) % q;+            }   +        }+    }+    else if(p == 5)+    {+        hDim_t temp1 = rts*4;+        hInt_t ru1 = ru[rustride*tupSize];+        hInt_t ru2 = ru[(rustride<<1)*tupSize];+        hInt_t ru3 = ru[(rustride*3)*tupSize];+        hInt_t 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;+                hInt_t 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]          = (y1 + ((ru1*y2) % q) + ((ru2*y3) % q) + ((ru3*y4) % q)) % q;+                y[(tensorOffset+rts)*tupSize]      = (y1 + ((ru2*y2) % q) + ((ru4*y3) % q) + ((ru1*y4) % q)) % q;+                y[(tensorOffset+(rts<<1))*tupSize] = (y1 + ((ru3*y2) % q) + ((ru1*y3) % q) + ((ru4*y4) % q)) % q;+                y[(tensorOffset+rts*3)*tupSize]    = (y1 + ((ru4*y2) % q) + ((ru3*y3) % q) + ((ru2*y4) % q)) % q;+            }   +        }+    }+    else if(p == 7)+    {+        hDim_t temp1 = rts*6;+        hInt_t ru1 = ru[rustride*tupSize];+        hInt_t ru2 = ru[(rustride<<1)*tupSize];+        hInt_t ru3 = ru[(rustride*3)*tupSize];+        hInt_t ru4 = ru[(rustride<<2)*tupSize];+        hInt_t ru5 = ru[(rustride*5)*tupSize];+        hInt_t 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;+                hInt_t 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]          = (y1 + ((ru1*y2) % q) + ((ru2*y3) % q) + ((ru3*y4) % q) + ((ru4*y5) % q) + ((ru5*y6) % q)) % q;+                y[(tensorOffset+rts)*tupSize]      = (y1 + ((ru2*y2) % q) + ((ru4*y3) % q) + ((ru6*y4) % q) + ((ru1*y5) % q) + ((ru3*y6) % q)) % q;+                y[(tensorOffset+(rts<<1))*tupSize] = (y1 + ((ru3*y2) % q) + ((ru6*y3) % q) + ((ru2*y4) % q) + ((ru5*y5) % q) + ((ru1*y6) % q)) % q;+                y[(tensorOffset+rts*3)*tupSize]    = (y1 + ((ru4*y2) % q) + ((ru1*y3) % q) + ((ru5*y4) % q) + ((ru2*y5) % q) + ((ru6*y6) % q)) % q;+                y[(tensorOffset+(rts<<2))*tupSize] = (y1 + ((ru5*y2) % q) + ((ru3*y3) % q) + ((ru1*y4) % q) + ((ru6*y5) % q) + ((ru4*y6) % q)) % q;+                y[(tensorOffset+rts*5)*tupSize]    = (y1 + ((ru6*y2) % q) + ((ru5*y3) % q) + ((ru4*y4) % q) + ((ru3*y5) % q) + ((ru2*y6) % q)) % q;+            }+        }+    }+    else+    {+        hInt_t* tempSpace = (hInt_t*)malloc((p-1)*sizeof(hInt_t));+        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++)+                {+                    hInt_t acc = 0;+                    for(hDim_t col = 0; col < p-1; col++)+                    {+                        acc += ((y[(tensorOffset+col*rts)*tupSize]*ru[((col*row) % p)*rustride*tupSize]) % q);+                    }+                    tempSpace[row-1] = acc % q;+                }+                +                for(hDim_t row = 0; row < p-1; row++)+                {+                    y[(tensorOffset+rts*row)*tupSize] = tempSpace[row];+                }+            }+        }+        free(tempSpace);+    }+}++//takes inverse rus+void crtpinvRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, +                hDim_t p, hDim_t rustride, hInt_t* ruinv, hInt_t q)+{+    if(p ==2)+    {+        // need this case so that we can divide overall by mhat^(-1)+        return;+    }+    else+    {+        hDim_t tensorOffset,i;+        hInt_t* tempSpace = (hInt_t*)malloc((p-1)*sizeof(hInt_t));+        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(i = 0; i < p-1; i++)+                {+                    hInt_t sum = 0;+                    int j;+                    for(j = 0; j < p-1; j++)+                    {+                        int ruIdx = ((j+1)*i) % p;+                        sum += ((y[(tensorOffset+j*rts)*tupSize] * ruinv[ruIdx*rustride*tupSize]) % q);+                    }+                    tempSpace[i] = sum % q;+                }++                hInt_t shift = 0;+                for(i = 0; i < p-1; i++)+                {+                    // we were given the inverse rus, so we need to negate the indices+                    shift += ((y[(tensorOffset+i*rts)*tupSize] * ruinv[rustride*(p-(i+1))*tupSize]) % q);+                }++                for(i = 0; i < p-1; i++)+                {+                    y[(tensorOffset+i*rts)*tupSize] = (tempSpace[i] - shift) % q; +                }+            }+        }+    }+}++void ppDFTRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, +              PrimeExponent pe, hDim_t rustride, hInt_t* ru, hInt_t q, hInt_t* 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;+        dftpRq (y, tupSize, lts*ltsScale, rtsDim, p, primeRuStride, ru, temp, q);+        dftptwidRq (y, tupSize, lts, rtsDim, pe, ltsScale*p, twidRuStride, ru, q);+        +        ltsScale /= p;+        rtsScale *= p;+        twidRuStride *= p;+        pe.exponent -= 1;+    }+}++void ppDFTInvRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, +                 PrimeExponent pe, hDim_t rustride, hInt_t* ru, hInt_t q, hInt_t* 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;+        dftptwidRq (y, tupSize, lts, rtsDim, pe, ltsScaleP, twidRuStride, ru, q);+        dftpRq (y, tupSize, lts*ltsScale, rtsDim, p, primeRuStride, ru, temp, q);+        +        ltsScale = ltsScaleP;+        rtsScale /= p;+        twidRuStride /= p;+        pe.exponent += 1;+    }+}++void ppcrtRq (void* y, hShort_t tupSize, hDim_t lts, hDim_t rts, +              PrimeExponent pe, void* ru, hInt_t* qs)+{+    hDim_t p = pe.prime;+    hDim_t e = pe.exponent;+#ifdef DEBUG_MODE+    ASSERT(e != 0);+#endif+    hDim_t mprime = ipow(p,e-1);+    +#ifdef DEBUG_MODE+    printf("lts is %" PRId32 "\trts is %" PRId32 "\n", lts, rts);+    printf("rus for p=%" PRId32 ", e=%" PRId16 "\t[", pe.prime, pe.exponent);+    hDim_t i;+    for(i = 0; i < ipow(p,e); i++) {+        printf("%" PRId64 ",", ((hInt_t*)ru)[i*tupSize+tupIdx]);+    }+    printf("]\n");+#endif++    hInt_t* temp = 0;+    if(p >= DFTP_GENERIC_SIZE)+    {+        temp = (hInt_t*)malloc(p*sizeof(hInt_t));+    }+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {+        hInt_t* z = ((hInt_t*)y)+tupIdx;+        hInt_t q = qs[tupIdx];+        hInt_t* ruOffset = ((hInt_t*)ru)+tupIdx;+        crtpRq (z, tupSize, lts*mprime, rts, p, mprime, ruOffset, q);+        crtTwiddleRq (z, tupSize, lts, rts, pe, ruOffset, q);+        pe.exponent -= 1;+        ppDFTRq (z,  tupSize, lts, rts*(p-1), pe, p, ruOffset, q, temp);+        pe.exponent += 1;+    }++    if(p >= DFTP_GENERIC_SIZE)+    {+        free(temp);+    }+}++void ppcrtinvRq (void* y, hShort_t tupSize, hDim_t lts, hDim_t rts, +                 PrimeExponent pe, void* ru, hInt_t* qs)+{+    hDim_t p = pe.prime;+    hDim_t e = pe.exponent;+#ifdef DEBUG_MODE+    ASSERT(e != 0);+#endif+    hDim_t mprime = ipow(p,e-1);+#ifdef DEBUG_MODE+    printf("lts is %" PRId32 "\trts is %" PRId32 "\n", lts, rts);+    printf("rus for p=%" PRId32 ", e=%" PRId16 "\t[", pe.prime, pe.exponent);+    hDim_t i;+    for(i = 0; i < ipow(p,e); i++) {+        printf("%" PRId64 ",", ((hInt_t*)ru)[i*tupSize+tupIdx]);+    }+    printf("]\n");+#endif++    hInt_t* temp = 0;+    if(p >= DFTP_GENERIC_SIZE)+    {+        temp = (hInt_t*)malloc(p*sizeof(hInt_t));+    }+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {+        hInt_t* z = ((hInt_t*)y)+tupIdx;+        hInt_t q = qs[tupIdx];+        hInt_t* ruOffset = ((hInt_t*)ru)+tupIdx;++        pe.exponent -= 1;+        ppDFTInvRq (z, tupSize, lts, rts*(p-1), pe, p, ruOffset, q, temp);+        pe.exponent += 1;+        crtTwiddleRq (z, tupSize, lts, rts, pe, ruOffset, q);+        crtpinvRq (z, tupSize, lts*mprime, rts, p, mprime, ruOffset, q);+    }++    if(p >= DFTP_GENERIC_SIZE)+    {+        free(temp);+    }+}++// EAC: Somebody who knows C/C++ should find a better way to handle pointers-to-pointers in a generic way+void tensorCRTRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t** ru, hInt_t* qs)+{+    hDim_t i;+#ifdef STATS+    struct timespec s1,s2,s3,s4,t1,t2,t3,t4;++    crtRqCtr++;++    clock_gettime(CLOCK_REALTIME, &s1);+    clock_gettime(CLOCK_MONOTONIC, &s2);+    clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &s3);+    clock_gettime(CLOCK_THREAD_CPUTIME_ID, &s4);+#endif+#ifdef DEBUG_MODE+    printf("\n\nEntered tensorCRTRq\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\ttupSize=%" PRId16 "\n[", totm, sizeOfPE,tupSize);+    for(i = 0; i < tupSize; i++) {+        printf("%" PRId64 ",", qs[i]);+    }+    printf("]\n[");+    for(i = 0; i < totm*tupSize; i++) {+        printf("%" PRId64 ",", y[i]);+    }+    printf("]\n[");+    for(i = 0; i < sizeOfPE; i++) {+        printf("(%" PRId32 ",%" PRId16 "),", peArr[i].prime, peArr[i].exponent);+    }+    printf("]\n");+#endif++    void** rus = (void**)malloc(sizeOfPE*sizeof(void*));+    for(i = 0; i < sizeOfPE; i++)+    {+        rus[i] = (void*) (ru[i]);+    }++  tensorFuserCRT (y, tupSize, ppcrtRq, totm, peArr, sizeOfPE, rus, qs);++    for(i = 0; i < tupSize; i++) {+        hInt_t q = qs[i];+      for(hDim_t j = 0; j < totm; j++)+      {+          if(y[j*tupSize+i]<0)+          {+              y[j*tupSize+i]+=q;+          }+#ifdef DEBUG_MODE+          if(y[j*tupSize+i]<0)+          {+              printf("TENSOR CRT^T INV\n");+          }+#endif+        }+  }++  free(rus);+#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    clock_gettime(CLOCK_MONOTONIC, &t2);+    clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &t3);+    clock_gettime(CLOCK_THREAD_CPUTIME_ID, &t4);++    crttime1 = tsAdd(crttime1, tsSubtract(t1,s1));+    crttime2 = tsAdd(crttime2, tsSubtract(t2,s2));+    crttime3 = tsAdd(crttime3, tsSubtract(t3,s3));+    crttime4 = tsAdd(crttime4, tsSubtract(t4,s4));+#endif+}++//takes inverse rus+void tensorCRTInvRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, +                    hInt_t** ruinv, hInt_t* mhatInv, hInt_t* qs)+{+  hDim_t i;+#ifdef STATS+    struct timespec s1,t1;+    crtInvRqCtr++;+    clock_gettime(CLOCK_REALTIME, &s1);+#endif+#ifdef DEBUG_MODE+    printf("\n\nEntered tensorCRTInvRq\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\ttupSize=%" PRId16 "\nqs[", totm, sizeOfPE, tupSize);+    for(i = 0; i < tupSize; i++) {+        printf("%" PRId64 ",", qs[i]);+    }+    printf("]\nmhatInv[");+    for(i = 0; i < tupSize; i++) {+        printf("%" PRId64 ",", mhatInv[i]);+    }+    printf("]\ny[");+    for(i = 0; i < totm*tupSize; i++) {+        printf("%" PRId64 ",", y[i]);+    }+    printf("]\npps[");+    for(i = 0; i < sizeOfPE; i++) {+        printf("(%" PRId32 ",%" PRId16 "),", peArr[i].prime, peArr[i].exponent);+    }+    printf("]\n");+#endif++  void** rus = (void**)malloc(sizeOfPE*sizeof(void*));+    for(i = 0; i < sizeOfPE; i++)+    {+        rus[i] = (void*) (ruinv[i]);+    }+  tensorFuserCRT (y, tupSize, ppcrtinvRq, totm, peArr, sizeOfPE, rus, qs);++    for (i = 0; i < tupSize; i++) {+        hInt_t q = qs[i];+      for (hDim_t j = 0; j < totm; j++)+      {+          y[j*tupSize+i] = (y[j*tupSize+i]*mhatInv[i])%q;+          if(y[j*tupSize+i] < 0)+          {+              y[j*tupSize+i] +=q;+          }+#ifdef DEBUG_MODE+          if(y[j*tupSize+i]<0)+          {+              printf("TENSOR CRT INV\n");+          }+#endif+        }+  }++  free(rus);+#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    crtInvRqTime = tsAdd(crtInvRqTime, tsSubtract(t1,s1));+#endif+}+++++++++++++++++++++++++++++++++++++++++void crtTwiddleC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, complex_t* ru)+{+    hDim_t idx;+    hDim_t p = pe.prime;+    hShort_t e = pe.exponent;+    +#ifdef DEBUG_MODE+    ASSERT(e != 0);+#endif++    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;+            complex_t 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++)+                {+                    idx = temp3 + modOffset;+                    CMPLX_IMUL(y[idx*tupSize],twid);+                }+            }+        }+    }+    else+    {+        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;+                complex_t 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++)+                    {+                        idx = temp3 + modOffset;+                        CMPLX_IMUL(y[idx*tupSize],twid);+                    }+                }+            }+        }+    }+}+    +// dim is power of p+void dftptwidC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, hDim_t dim, hDim_t rustride, complex_t* 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);+            complex_t 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;+                    CMPLX_IMUL(y[idx*tupSize],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);+                complex_t 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;+                        CMPLX_IMUL(y[idx*tupSize],twid);+                    }+                }+            }+        }+    }+}++//implied length of ru is rustride*p+//implied length of tempSpace is p, if p is not a special case+void dftpC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hDim_t rustride, complex_t* ru, complex_t* tempSpace)+{+    hDim_t blockOffset, modOffset, tensorOffset;+    +    if(p == 2)+    {+        hDim_t temp1 = rts<<1;+        for(blockOffset = 0; blockOffset < lts; blockOffset++)+        {+            hDim_t temp2 = blockOffset*temp1;+            for(modOffset = 0; modOffset < rts; modOffset++)+            {+                tensorOffset = temp2 + modOffset;+                complex_t u = y[tensorOffset*tupSize];+                complex_t t = y[(tensorOffset+rts)*tupSize];+                y[tensorOffset*tupSize] = CMPLX_ADD(u,t);+                y[(tensorOffset+rts)*tupSize] = CMPLX_SUB(u,t);+            }+        }+    }+    else if(p == 3)+    {+        hDim_t temp1 = rts*3;+        complex_t ru1 = ru[rustride*tupSize];+        complex_t ru2 = ru[(rustride<<1)*tupSize];++        for(blockOffset = 0; blockOffset < lts; blockOffset++)+        {+            hDim_t temp2 = blockOffset*temp1;+            for(modOffset = 0; modOffset < rts; modOffset++)+            {+                tensorOffset = temp2 + modOffset;+                complex_t y1, y2, y3;+                y1 = y[tensorOffset*tupSize];+                y2 = y[(tensorOffset+rts)*tupSize];+                y3 = y[(tensorOffset+(rts<<1))*tupSize];+                y[tensorOffset*tupSize]            = CMPLX_ADD3(y1,               y2,                y3);+                y[(tensorOffset+rts)*tupSize]      = CMPLX_ADD3(y1, CMPLX_MUL(ru1,y2), CMPLX_MUL(ru2,y3));+                y[(tensorOffset+(rts<<1))*tupSize] = CMPLX_ADD3(y1, CMPLX_MUL(ru2,y2), CMPLX_MUL(ru1,y3));+            }   +        }+    }+    else if(p == 5)+    {+        hDim_t temp1 = rts*5;+        complex_t ru1 = ru[rustride*tupSize];+        complex_t ru2 = ru[(rustride<<1)*tupSize];+        complex_t ru3 = ru[(rustride*3)*tupSize];+        complex_t ru4 = ru[(rustride<<2)*tupSize];++        for(blockOffset = 0; blockOffset < lts; blockOffset++)+        {+            hDim_t temp2 = blockOffset*temp1;+            for(modOffset = 0; modOffset < rts; modOffset++)+            {+                tensorOffset = temp2 + modOffset;+                complex_t 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]          = CMPLX_ADD5(y1,               y2,                y3,                y4,                y5);+                y[(tensorOffset+rts)*tupSize]      = CMPLX_ADD5(y1, CMPLX_MUL(ru1,y2), CMPLX_MUL(ru2,y3), CMPLX_MUL(ru3,y4), CMPLX_MUL(ru4,y5));+                y[(tensorOffset+(rts<<1))*tupSize] = CMPLX_ADD5(y1, CMPLX_MUL(ru2,y2), CMPLX_MUL(ru4,y3), CMPLX_MUL(ru1,y4), CMPLX_MUL(ru3,y5));+                y[(tensorOffset+rts*3)*tupSize]    = CMPLX_ADD5(y1, CMPLX_MUL(ru3,y2), CMPLX_MUL(ru1,y3), CMPLX_MUL(ru4,y4), CMPLX_MUL(ru2,y5));+                y[(tensorOffset+(rts<<2))*tupSize] = CMPLX_ADD5(y1, CMPLX_MUL(ru4,y2), CMPLX_MUL(ru3,y3), CMPLX_MUL(ru2,y4), CMPLX_MUL(ru1,y5));+            }   +        }+    }+    else if(p == 7)+    {+        hDim_t temp1 = rts*7;+        complex_t ru1 = ru[rustride*tupSize];+        complex_t ru2 = ru[(rustride<<1)*tupSize];+        complex_t ru3 = ru[(rustride*3)*tupSize];+        complex_t ru4 = ru[(rustride<<2)*tupSize];+        complex_t ru5 = ru[(rustride*5)*tupSize];+        complex_t ru6 = ru[(rustride*6)*tupSize];++        for(blockOffset = 0; blockOffset < lts; blockOffset++)+        {+            hDim_t temp2 = blockOffset*temp1;+            for(modOffset = 0; modOffset < rts; modOffset++)+            {+                tensorOffset = temp2 + modOffset;+                complex_t 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]          = CMPLX_ADD7(y1,               y2,                y3,                y4,                y5,                y6,                y7);+                y[(tensorOffset+rts)*tupSize]      = CMPLX_ADD7(y1, CMPLX_MUL(ru1,y2), CMPLX_MUL(ru2,y3), CMPLX_MUL(ru3,y4), CMPLX_MUL(ru4,y5), CMPLX_MUL(ru5,y6), CMPLX_MUL(ru6,y7));+                y[(tensorOffset+(rts<<1))*tupSize] = CMPLX_ADD7(y1, CMPLX_MUL(ru2,y2), CMPLX_MUL(ru4,y3), CMPLX_MUL(ru6,y4), CMPLX_MUL(ru1,y5), CMPLX_MUL(ru3,y6), CMPLX_MUL(ru5,y7));+                y[(tensorOffset+rts*3)*tupSize]    = CMPLX_ADD7(y1, CMPLX_MUL(ru3,y2), CMPLX_MUL(ru6,y3), CMPLX_MUL(ru2,y4), CMPLX_MUL(ru5,y5), CMPLX_MUL(ru1,y6), CMPLX_MUL(ru4,y7));+                y[(tensorOffset+(rts<<2))*tupSize] = CMPLX_ADD7(y1, CMPLX_MUL(ru4,y2), CMPLX_MUL(ru1,y3), CMPLX_MUL(ru5,y4), CMPLX_MUL(ru2,y5), CMPLX_MUL(ru6,y6), CMPLX_MUL(ru3,y7));+                y[(tensorOffset+rts*5)*tupSize]    = CMPLX_ADD7(y1, CMPLX_MUL(ru5,y2), CMPLX_MUL(ru3,y3), CMPLX_MUL(ru1,y4), CMPLX_MUL(ru6,y5), CMPLX_MUL(ru4,y6), CMPLX_MUL(ru2,y7));+                y[(tensorOffset+rts*6)*tupSize]    = CMPLX_ADD7(y1, CMPLX_MUL(ru6,y2), CMPLX_MUL(ru5,y3), CMPLX_MUL(ru4,y4), CMPLX_MUL(ru3,y5), CMPLX_MUL(ru2,y6), CMPLX_MUL(ru1,y7));+            }   +        }+    }+    else+    {+        hDim_t temp1 = rts*p;+        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; row++)+                {+                    complex_t acc = ((complex_t){0,0});+                    for(col = 0; col < p; col++)+                    {+                        CMPLX_IADD(acc, CMPLX_MUL(y[(tensorOffset+col*rts)*tupSize], ru[((col*row) % p)*rustride*tupSize]));+                    }+                    tempSpace[row] = acc;+                }+                +                for(row = 0; row < p; row++)+                {+                    y[(tensorOffset+rts*row)*tupSize] = tempSpace[row];   +                }+            }+        }+    }+}++void crtpC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hDim_t rustride, complex_t* ru)+{+    hDim_t blockOffset, modOffset, tensorOffset;+    +    if(p == 2)+    {+        return;+    }+    else if(p == 3)+    {+        hDim_t temp1 = rts*2;+        complex_t ru1 = ru[rustride*tupSize];+        complex_t ru2 = ru[(rustride<<1)*tupSize];++        for(blockOffset = 0; blockOffset < lts; blockOffset++)+        {+            hDim_t temp2 = blockOffset*temp1;+            for(modOffset = 0; modOffset < rts; modOffset++)+            {+                tensorOffset = temp2 + modOffset;+                complex_t y1, y2;+                y1 = y[tensorOffset*tupSize];+                y2 = y[(tensorOffset+rts)*tupSize];+                y[tensorOffset*tupSize]     = CMPLX_ADD(y1, CMPLX_MUL(ru1,y2));+                y[(tensorOffset+rts)*tupSize] = CMPLX_ADD(y1, CMPLX_MUL(ru2,y2));+            }   +        }+    }+    else if(p == 5)+    {+        hDim_t temp1 = rts*4;+        complex_t ru1 = ru[rustride*tupSize];+        complex_t ru2 = ru[(rustride<<1)*tupSize];+        complex_t ru3 = ru[(rustride*3)*tupSize];+        complex_t ru4 = ru[(rustride<<2)*tupSize];++        for(blockOffset = 0; blockOffset < lts; blockOffset++)+        {+            hDim_t temp2 = blockOffset*temp1;+            for(modOffset = 0; modOffset < rts; modOffset++)+            {+                tensorOffset = temp2 + modOffset;+                complex_t 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]          = CMPLX_ADD4(y1, CMPLX_MUL(ru1,y2), CMPLX_MUL(ru2,y3), CMPLX_MUL(ru3,y4));+                y[(tensorOffset+rts)*tupSize]      = CMPLX_ADD4(y1, CMPLX_MUL(ru2,y2), CMPLX_MUL(ru4,y3), CMPLX_MUL(ru1,y4));+                y[(tensorOffset+(rts<<1))*tupSize] = CMPLX_ADD4(y1, CMPLX_MUL(ru3,y2), CMPLX_MUL(ru1,y3), CMPLX_MUL(ru4,y4));+                y[(tensorOffset+rts*3)*tupSize]    = CMPLX_ADD4(y1, CMPLX_MUL(ru4,y2), CMPLX_MUL(ru3,y3), CMPLX_MUL(ru2,y4));+            }   +        }+    }+    else if(p == 7)+    {+        hDim_t temp1 = rts*6;+        complex_t ru1 = ru[rustride*tupSize];+        complex_t ru2 = ru[(rustride<<1)*tupSize];+        complex_t ru3 = ru[(rustride*3)*tupSize];+        complex_t ru4 = ru[(rustride<<2)*tupSize];+        complex_t ru5 = ru[(rustride*5)*tupSize];+        complex_t ru6 = ru[(rustride*6)*tupSize];++        for(blockOffset = 0; blockOffset < lts; blockOffset++)+        {+            hDim_t temp2 = blockOffset*temp1;+            for(modOffset = 0; modOffset < rts; modOffset++)+            {+                tensorOffset = temp2 + modOffset;+                complex_t 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]          = CMPLX_ADD6(y1, CMPLX_MUL(ru1,y2), CMPLX_MUL(ru2,y3), CMPLX_MUL(ru3,y4), CMPLX_MUL(ru4,y5), CMPLX_MUL(ru5,y6));+                y[(tensorOffset+rts)*tupSize]      = CMPLX_ADD6(y1, CMPLX_MUL(ru2,y2), CMPLX_MUL(ru4,y3), CMPLX_MUL(ru6,y4), CMPLX_MUL(ru1,y5), CMPLX_MUL(ru3,y6));+                y[(tensorOffset+(rts<<1))*tupSize] = CMPLX_ADD6(y1, CMPLX_MUL(ru3,y2), CMPLX_MUL(ru6,y3), CMPLX_MUL(ru2,y4), CMPLX_MUL(ru5,y5), CMPLX_MUL(ru1,y6));+                y[(tensorOffset+rts*3)*tupSize]    = CMPLX_ADD6(y1, CMPLX_MUL(ru4,y2), CMPLX_MUL(ru1,y3), CMPLX_MUL(ru5,y4), CMPLX_MUL(ru2,y5), CMPLX_MUL(ru6,y6));+                y[(tensorOffset+(rts<<2))*tupSize] = CMPLX_ADD6(y1, CMPLX_MUL(ru5,y2), CMPLX_MUL(ru3,y3), CMPLX_MUL(ru1,y4), CMPLX_MUL(ru6,y5), CMPLX_MUL(ru4,y6));+                y[(tensorOffset+rts*5)*tupSize]    = CMPLX_ADD6(y1, CMPLX_MUL(ru6,y2), CMPLX_MUL(ru5,y3), CMPLX_MUL(ru4,y4), CMPLX_MUL(ru3,y5), CMPLX_MUL(ru2,y6));+            }   +        }+    }+    else+    {+        complex_t* tempSpace = (complex_t*)malloc((p-1)*sizeof(complex_t));+        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 = 1; row < p; row++)+                {+                    complex_t acc = ((complex_t){0,0});+                    for(col = 0; col < p-1; col++)+                    {+                        CMPLX_IADD(acc, CMPLX_MUL(y[(tensorOffset+col*rts)*tupSize], ru[((col*row) % p)*rustride*tupSize]));+                    }+                    tempSpace[row-1] = acc;+                }+                +                for(row = 0; row < p-1; row++)+                {+                    y[(tensorOffset+rts*row)*tupSize] = tempSpace[row];   +                }+            }+        }+        free(tempSpace);+    }+}++//takes inverse rus+void crtpinvC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hDim_t rustride, complex_t* ruinv)+{+    if(p ==2)+    {+        // need this case so that we can divide overall by mhat^(-1)+        return;+    }+    else+    {+        hDim_t tensorOffset,i;+        complex_t* tempSpace = (complex_t*)malloc(p*sizeof(complex_t));+        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(i = 0; i < p-1; i++)+                {+                    complex_t sum = ((complex_t){0,0});+                    int j;+                    for(j = 0; j < p-1; j++)+                    {+                        int ruIdx = (((j+1)*i) % p)*rustride;+                        CMPLX_IADD(sum, CMPLX_MUL(y[(tensorOffset+j*rts)*tupSize],ruinv[ruIdx*tupSize]));+                    }+                    tempSpace[i] = sum;+                }++                complex_t shift = ((complex_t){0,0});+                for(i = 0; i < p-1; i++)+                {+                    // we were given the inverse rus, so we need to negate the indices+                    int ruIdx = p-(i+1);+                    CMPLX_IADD(shift, CMPLX_MUL(y[(tensorOffset+i*rts)*tupSize], ruinv[rustride*ruIdx*tupSize]));+                }++                for(i = 0; i < p-1; i++)+                {+                    y[(tensorOffset+i*rts)*tupSize] = CMPLX_SUB(tempSpace[i], shift); +                }+            }+        }+    }+}++void ppDFTC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, hDim_t rustride, complex_t* ru, complex_t* 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;+        dftpC (y, tupSize, lts*ltsScale, rtsDim, p, primeRuStride, ru, temp);+        dftptwidC (y, tupSize, lts, rtsDim, pe, ltsScale*p, twidRuStride, ru);+        +        ltsScale /= p;+        rtsScale *= p;+        twidRuStride *= p;+        pe.exponent -= 1;+    }+}++void ppDFTInvC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, hDim_t rustride, complex_t* ru, complex_t* 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;+        dftptwidC (y, tupSize, lts, rtsDim, pe, ltsScaleP, twidRuStride, ru);+        dftpC (y, tupSize, lts*ltsScale, rtsDim, p, primeRuStride, ru, temp);+        +        ltsScale = ltsScaleP;+        rtsScale /= p;+        twidRuStride /= p;+        pe.exponent += 1;+    }+}++void ppcrtC (void* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, void* ru, hInt_t* qs)+{+    hDim_t p = pe.prime;+    hDim_t e = pe.exponent;+#ifdef DEBUG_MODE+    ASSERT(e != 0);+#endif+    hDim_t mprime = ipow(p,e-1);++#ifdef DEBUG_MODE+    printf("rus for p=%" PRId32 ", e=%" PRId16 "\t[", pe.prime, pe.exponent);+    hDim_t i;+    for(i = 0; i < ipow(p,e); i++) {+        printf("(%f,%f),", ((complex_t*)ru)[i].real, ((complex_t*)ru)[i].imag);+    }+    printf("]\n");+#endif++    complex_t* temp = 0;+    if(p >= DFTP_GENERIC_SIZE)+    {+        temp = (complex_t*)malloc(p*sizeof(complex_t));+    }++    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {+        complex_t* z = ((complex_t*)y)+tupIdx;+        complex_t* ruOffset = ((complex_t*)ru)+tupIdx;+        crtpC (z, tupSize, lts*mprime, rts, p, mprime, ruOffset);+        crtTwiddleC (z, tupSize, lts, rts, pe, ruOffset);+        pe.exponent -= 1;+        ppDFTC (z, tupSize, lts, rts*(p-1), pe, p, ruOffset, temp);+        pe.exponent += 1;+    }++    if(p >= DFTP_GENERIC_SIZE)+    {+        free(temp);+    }+}++void ppcrtinvC (void* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, void* ru, hInt_t* qs)+{+    hDim_t p = pe.prime;+    hDim_t e = pe.exponent;+#ifdef DEBUG_MODE+    ASSERT(e != 0);+#endif+    hDim_t mprime = ipow(p,e-1);+    +    complex_t* temp = 0;+    if(p >= DFTP_GENERIC_SIZE)+    {+        temp = (complex_t*)malloc(p*sizeof(complex_t));+    }++    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {+        complex_t* z = ((complex_t*)y)+tupIdx;+        complex_t* ruOffset = ((complex_t*)ru)+tupIdx;+        pe.exponent -= 1;+        ppDFTInvC (z, tupSize, lts, rts*(p-1), pe, p, ruOffset, temp);+        pe.exponent += 1;+        crtTwiddleC (z, tupSize, lts, rts, pe, ruOffset);+        crtpinvC (z, tupSize, lts*mprime, rts, p, mprime, ruOffset);+    }++    if(p >= DFTP_GENERIC_SIZE)+    {+        free(temp);+    }+}++void tensorCRTC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, complex_t** ru)+{+#ifdef STATS+    struct timespec s1,t1;+    crtCCtr++;+    clock_gettime(CLOCK_REALTIME, &s1);+#endif+#ifdef DEBUG_MODE+    printf("\n\nEntered tensorCRTC\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\n[", totm, sizeOfPE);+    hDim_t j;+    for(j = 0; j < totm; j++) {+        printf("(%f,%f),", y[j].real, y[j].imag);+    }+    printf("]\n[");+    for(j = 0; j < sizeOfPE; j++) {+        printf("(%" PRId32 ",%" PRId16 "),", peArr[j].prime, peArr[j].exponent);+    }+    printf("]\n");+#endif+    void** rus = (void**)malloc(sizeOfPE*sizeof(void*));+    hShort_t i;+    for(i = 0; i < sizeOfPE; i++)+    {+        rus[i] = (void*) (ru[i]);+    }+  tensorFuserCRT (y, tupSize, ppcrtC, totm, peArr, sizeOfPE, rus, (hInt_t*)0);+  free(rus);+#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    crtCTime = tsAdd(crtCTime, tsSubtract(t1,s1));+#endif+}++//takes inverse rus+void tensorCRTInvC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, +                    hShort_t sizeOfPE, complex_t** ruinv, complex_t* mhatInv)+{+#ifdef STATS+    struct timespec s1,t1;+    crtInvCCtr++;+    clock_gettime(CLOCK_REALTIME, &s1);+#endif+  hDim_t i;+  +  void** rus = (void**)malloc(sizeOfPE*sizeof(void*));+    for(i = 0; i < sizeOfPE; i++)+    {+        rus[i] = (void*) (ruinv[i]);+    }+  +  tensorFuserCRT (y, tupSize, ppcrtinvC, totm, peArr, sizeOfPE, rus, (hInt_t*)0);++    for (i = 0; i < tupSize; i++) {+      for (hDim_t j = 0; j < totm; j++)+      {+          CMPLX_IMUL(y[j*tupSize+i], mhatInv[i]);+      }+    }+  +  free(rus);+#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    crtInvCTime = tsAdd(crtInvCTime, tsSubtract(t1,s1));+#endif+}

@@ -0,0 +1,828 @@+#include "tensorTypes.h"
+
+
+void gPowR (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p)
+{
+  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;
+      hInt_t 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;
+    }
+  }
+}
+
+void gPowRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hInt_t q)
+{
+  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;
+      hInt_t last = y[(tensorOffset + tmp2)*tupSize];
+      for (i = p-2; i != 0; --i)
+      {
+        hDim_t idx = tensorOffset + i * rts;
+        y[idx*tupSize] = (y[idx*tupSize] + last - y[(idx-rts)*tupSize]) % q;
+      }
+      y[tensorOffset*tupSize] = (y[tensorOffset*tupSize] + last) % q;
+    }
+  }
+}
+
+void gPowC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p)
+{
+  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;
+      complex_t last = y[(tensorOffset + tmp2)*tupSize];
+      for (i = p-2; i != 0; --i)
+      {
+        hDim_t idx = tensorOffset + i * rts;
+        CMPLX_IADD(y[idx*tupSize],last);
+        CMPLX_ISUB(y[idx*tupSize],y[(idx-rts)*tupSize]);
+      }
+      CMPLX_IADD(y[tensorOffset*tupSize],last);
+    }
+  }
+}
+
+void ppGPowR (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs)
+{
+#ifdef DEBUG_MODE
+  ASSERT (q==0);
+#endif
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+     
+    if (p != 2)
+    {
+      for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+        gPowR (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);
+      }
+    }
+}
+
+void ppGPowRq (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs)
+{
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if (p != 2)
+    {
+      for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+        gPowRq (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p, qs[tupIdx]);
+      }
+    }
+}
+
+void ppGPowC (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs)
+{
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if (p != 2)
+    {
+      for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+        gPowC (((complex_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);
+      }
+    }
+}
+
+void gDecR (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p)
+{
+  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;
+      hInt_t 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;
+    }
+  }
+}
+
+void gDecRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hInt_t q)
+{
+  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;
+      hInt_t acc = y[tensorOffset*tupSize];
+      for (i = p-2; i != 0; --i)
+      {
+        hDim_t idx = tensorOffset + i * rts;
+        // acc is at most p*q << 64 bits, so no need to mod
+        acc = acc + y[idx*tupSize];
+        y[idx*tupSize] = (y[idx*tupSize] - y[(idx-rts)*tupSize]) % q;
+      }
+      y[tensorOffset*tupSize] = (y[tensorOffset*tupSize] + acc) % q;
+    }
+  }
+}
+
+void ppGDecR (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs)
+{
+#ifdef DEBUG_MODE
+  ASSERT (q==0);
+#endif
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if (p != 2)
+    {
+      for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+        gDecR (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);
+      }
+    }
+}
+
+void ppGDecRq (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs)
+{
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if (p != 2)
+    {
+      for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+        gDecRq (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p, qs[tupIdx]);
+      }
+    }
+}
+
+
+void gInvPowR (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p)
+{
+  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;
+      hInt_t lelts = 0;
+      for (i = 0; i < p-1; ++i)
+      {
+        lelts += y[(tensorOffset + i*rts)*tupSize];
+      }
+      hInt_t relts = 0;
+      for (i = p-2; i >= 0; --i)
+      {
+        hDim_t idx = tensorOffset + i*rts;
+        hInt_t z = y[idx*tupSize];
+        y[idx*tupSize] = (p-1-i) * lelts - (i+1)*relts;
+        lelts -= z;
+        relts += z;
+      }
+    }
+  }
+}
+
+void gInvPowRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hInt_t q)
+{
+  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;
+      hInt_t lelts = 0;
+      //lelts is at most p*q, so we can mod once at the end
+      for (i = 0; i < p-1; ++i)
+      {
+        lelts = lelts + y[(tensorOffset + i*rts)*tupSize];
+      }
+      lelts = lelts % q;
+      //in the next loop, lelts <= p*q and relts <= p*q
+      //products are <= p*p*q, and diff is <= 2*p*p*q
+      //so we assume 2*p^2 << 31 bits
+      hInt_t relts = 0;
+      for (i = p-2; i >= 0; --i)
+      {
+        hDim_t idx = tensorOffset + i*rts;
+        hInt_t z = y[idx*tupSize];
+        y[idx*tupSize] = (((p-1-i) * lelts) - ((i+1)*relts)) % q;
+        lelts -= z;
+        relts += z;
+      }
+    }
+  }
+}
+
+void gInvPowC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p)
+{
+  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;
+      complex_t lelts = ((complex_t){0, 0});;
+      for (i = 0; i < p-1; ++i)
+      {
+        CMPLX_IADD(lelts,y[(tensorOffset + i*rts)*tupSize]);
+      }
+      complex_t relts = ((complex_t){0, 0});;
+      for (i = p-2; i >= 0; --i)
+      {
+        hDim_t idx = tensorOffset + i*rts;
+        complex_t z = y[idx*tupSize];
+
+        complex_t c1 = ((complex_t){p-1-i, 0});
+        complex_t c2 = ((complex_t){i+1, 0});
+        complex_t t1 = CMPLX_MUL(c1, lelts);
+        complex_t t2 = CMPLX_MUL(c2, relts);
+        y[idx*tupSize] = CMPLX_SUB(t1,t2);
+        CMPLX_ISUB(lelts,z);
+        CMPLX_IADD(relts,z);
+      }
+    }
+  }
+}
+
+void ppGInvPowR (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs)
+{
+#ifdef DEBUG_MODE
+  ASSERT (q==0);
+#endif
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if (p != 2)
+    {
+      for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+        gInvPowR (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);
+      }
+    }
+}
+
+void ppGInvPowRq (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs)
+{
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if (p != 2)
+    {
+      for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+        gInvPowRq (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p, qs[tupIdx]);
+      }
+    }
+}
+
+void ppGInvPowC (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs)
+{
+#ifdef DEBUG_MODE
+  ASSERT (q==0);
+#endif
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if (p != 2)
+    {
+      for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+        gInvPowC (((complex_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);
+      }
+    }
+}
+
+//do not call for p=2!
+void gCRTRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hInt_t* gcoeffs, hInt_t q)
+{
+    hDim_t gindex;
+    hDim_t blockOffset, modOffset, idx;
+    hDim_t temp1 = rts*(p-1);
+    
+    for(gindex = 0; gindex < p-1; gindex++)
+    {
+        hInt_t coeff = gcoeffs[gindex*tupSize];
+        hDim_t temp3 = gindex*rts;
+        for(blockOffset = 0; blockOffset < lts; blockOffset++)
+        {
+            hDim_t temp2 = blockOffset*temp1 + temp3;
+            for(modOffset = 0; modOffset < rts; modOffset++)
+            {
+                idx = temp2 + modOffset;
+                y[idx*tupSize] = (y[idx*tupSize]*coeff)%q;
+            }
+        }
+    }
+}
+
+//do not call for p=2!
+void gCRTC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, complex_t* gcoeffs)
+{
+    hDim_t gindex;
+    hDim_t blockOffset, modOffset, idx;
+    hDim_t temp1 = rts*(p-1);
+    
+    for(gindex = 0; gindex < p-1; gindex++)
+    {
+        complex_t coeff = gcoeffs[gindex*tupSize];
+        hDim_t temp3 = gindex*rts;
+        for(blockOffset = 0; blockOffset < lts; blockOffset++)
+        {
+            hDim_t temp2 = blockOffset*temp1 + temp3;
+            for(modOffset = 0; modOffset < rts; modOffset++)
+            {
+                idx = temp2 + modOffset;
+                CMPLX_IMUL(y[idx*tupSize],coeff);
+            }
+        }
+    }
+}
+
+void ppGCRTRq (void* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, void* gcoeffs, hInt_t* qs)
+{
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    
+#ifdef DEBUG_MODE
+    printf("gcoeffs for p=%" PRId32 ", e=%" PRId16 "\t[", pe.prime, pe.exponent);
+    int i;
+    for(i = 0; i < ((p-1)*ipow(p,e-1)); i++) {
+        printf("%" PRId64 ",", ((hInt_t*)gcoeffs)[i]);
+    }
+    printf("]\n");
+#endif
+    
+    if (p != 2)
+    {
+      for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+        gCRTRq (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p, ((hInt_t*)gcoeffs)+tupIdx, qs[tupIdx]);
+      }
+    }
+}
+
+void ppGCRTC (void* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, void* gcoeffs, hInt_t* qs)
+{
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    
+#ifdef DEBUG_MODE
+    printf("gcoeffs for p=%" PRId32 ", e=%" PRId16 "\t[", pe.prime, pe.exponent);
+    int i;
+    for(i = 0; i < ((p-1)*ipow(p,e-1)); i++) {
+        printf("(%f,%f),", ((complex_t*)gcoeffs)[i].real, ((complex_t*)gcoeffs)[i].imag);
+    }
+    printf("]\n");
+#endif
+    
+    if (p != 2)
+    {
+      for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+        gCRTC (((complex_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p, ((complex_t*)gcoeffs)+tupIdx);
+      }
+    }
+}
+
+void gInvDecR (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p)
+{
+  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;
+      hInt_t lastOut = 0;
+      for (i=1; i < p; ++i)
+      {
+        lastOut += i * y[(tensorOffset + (i-1)*rts)*tupSize];
+      }
+      hInt_t acc = lastOut / p;
+      ASSERT (acc * p == lastOut);  // this line asserts that lastOut % p == 0, without calling % operator
+      for (i = p-2; i > 0; --i)
+      {
+        hDim_t idx = tensorOffset + i*rts;
+        hInt_t tmp = acc;
+        acc -= y[idx*tupSize]; // we already divided acc by p, do not multiply y[idx] by p
+        y[idx*tupSize] = tmp;
+      }
+      y[tensorOffset*tupSize] = acc;
+    }
+  }
+}
+
+void gInvDecRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hInt_t q)
+{
+  hDim_t blockOffset;
+  hDim_t modOffset;
+  hDim_t i;
+  hDim_t tmp1 = rts*(p-1);
+  hInt_t reciprocalOfP = reciprocal (q,p);
+
+  for (blockOffset = 0; blockOffset < lts; ++blockOffset)
+  {
+    hDim_t tmp2 = blockOffset*tmp1;
+    for (modOffset = 0; modOffset < rts; ++modOffset)
+    {
+      hDim_t tensorOffset = tmp2 + modOffset;
+      hInt_t lastOut = 0;
+      for (i=1; i < p; ++i)
+      {
+        lastOut += (i * y[(tensorOffset + (i-1)*rts)*tupSize]);
+      }
+      //in the previous loop, |lastOut| <= p*p*q
+      lastOut = lastOut % q;
+      hInt_t acc = (lastOut * reciprocalOfP) % q;
+      // |acc| <= p*q
+      for (i = p-2; i > 0; --i)
+      {
+        hDim_t idx = tensorOffset + i*rts;
+        hInt_t tmp = acc;
+        acc = acc - y[idx*tupSize];
+        y[idx*tupSize] = tmp % q;
+      }
+      y[tensorOffset*tupSize] = acc % q;
+    }
+  }
+}
+
+void ppGInvDecR (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs)
+{
+#ifdef DEBUG_MODE
+  ASSERT (q==0);
+#endif
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if (p != 2)
+    {
+      for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+        gInvDecR (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);
+      }
+    }
+}
+
+void ppGInvDecRq (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs)
+{
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if (p != 2)
+    {
+      for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+        gInvDecRq (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p, qs[tupIdx]);
+      }
+    }
+}
+
+#ifdef STATS
+int gprCtr = 0;
+int gprqCtr = 0;
+int gdrCtr = 0;
+int gdrqCtr = 0;
+int giprCtr = 0;
+int giprqCtr = 0;
+int gidrCtr = 0;
+int gidrqCtr = 0;
+int gcrqCtr = 0;
+int gccCtr = 0;
+int gicrqCtr = 0;
+int giccCtr = 0;
+
+struct timespec gprTime = {0,0};
+struct timespec gprqTime = {0,0};
+struct timespec gdrTime = {0,0};
+struct timespec gdrqTime = {0,0};
+struct timespec giprTime = {0,0};
+struct timespec giprqTime = {0,0};
+struct timespec gidrTime = {0,0};
+struct timespec gidrqTime = {0,0};
+struct timespec gcrqTime = {0,0};
+struct timespec gccTime = {0,0};
+#endif
+
+void tensorGPowR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE)
+{
+#ifdef STATS
+    gprCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppGPowR, totm, peArr, sizeOfPE, (hInt_t*)0);
+
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    gprTime = tsAdd(gprTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorGPowRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs)
+{
+#ifdef STATS
+    gprqCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppGPowRq, totm, peArr, sizeOfPE, qs);
+
+  hDim_t j;
+  for(int tupIdx = 0; tupIdx<tupSize; tupIdx++) {
+    hInt_t q = qs[tupIdx];
+    for(j = 0; j < totm; j++)
+    {
+        if(y[j*tupSize+tupIdx]<0)
+        {
+            y[j*tupSize+tupIdx]+=q;
+        }
+    }
+  }
+
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    gprqTime = tsAdd(gprqTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorGPowC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE)
+{
+#ifdef STATS
+    gpcCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppGPowC, totm, peArr, sizeOfPE, (hInt_t*)0);
+
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    gpcTime = tsAdd(gpcTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorGDecR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE)
+{
+#ifdef STATS
+    gdrCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppGDecR, totm, peArr, sizeOfPE, (hInt_t*)0);
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    gdrTime = tsAdd(gdrTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorGDecRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs)
+{
+#ifdef STATS
+    gdrqCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppGDecRq, totm, peArr, sizeOfPE, qs);
+
+  hDim_t j;
+  for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+    hInt_t q = qs[tupIdx];
+    for(j = 0; j < totm; j++)
+    {
+        if(y[j*tupSize+tupIdx]<0)
+        {
+            y[j*tupSize+tupIdx]+=q;
+        }
+    }
+  }
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    gdrqTime = tsAdd(gdrqTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorGInvPowR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE)
+{
+#ifdef STATS
+    giprCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppGInvPowR, totm, peArr, sizeOfPE, (hInt_t*)0);
+
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    giprTime = tsAdd(giprTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorGInvPowRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs)
+{
+#ifdef STATS
+    giprqCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppGInvPowRq, totm, peArr, sizeOfPE, qs);
+
+  hDim_t j;
+  for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+    hInt_t q = qs[tupIdx];
+    for(j = 0; j < totm; j++)
+    {
+        if(y[j*tupSize+tupIdx]<0)
+        {
+            y[j*tupSize+tupIdx]+=q;
+        }
+    }
+  }
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    giprqTime = tsAdd(giprqTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorGInvPowC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE)
+{
+#ifdef STATS
+    gipcCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppGInvPowC, totm, peArr, sizeOfPE, (hInt_t*)0);
+
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    gipcTime = tsAdd(gipcTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorGInvDecR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE)
+{
+#ifdef STATS
+    gidrCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppGInvDecR, totm, peArr, sizeOfPE, (hInt_t*)0);
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    gidrTime = tsAdd(gidrTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorGInvDecRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs)
+{
+#ifdef STATS
+    gidrqCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+    tensorFuser (y, tupSize, ppGInvDecRq, totm, peArr, sizeOfPE, qs);
+
+  hDim_t j;
+  for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+    hInt_t q = qs[tupIdx];
+    for(j = 0; j < totm; j++)
+    {
+        if(y[j*tupSize+tupIdx]<0)
+        {
+            y[j*tupSize+tupIdx]+=q;
+        }
+    }
+  }
+
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    gidrqTime = tsAdd(gidrqTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorGCRTRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t** gcoeffs, hInt_t* qs)
+{
+#ifdef STATS
+    gcrqCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+#ifdef DEBUG_MODE
+    printf("\n\nEntered tensorGCRTRq\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\tq=%" PRId64 "\n[", totm, sizeOfPE, q);
+    hDim_t j;
+    for(j = 0; j < totm; j++) {
+        printf("%" PRId64 ",", y[j]);
+    }
+    printf("]\n[");
+    for(j = 0; j < sizeOfPE; j++) {
+        printf("(%" PRId32 ",%" PRId16 "),", peArr[j].prime, peArr[j].exponent);
+    }
+    printf("]\n");
+#endif
+    void** vgcoeffs = (void**)malloc(sizeOfPE*sizeof(void*));
+    hDim_t i;
+    for(i = 0; i < sizeOfPE; i++)
+    {
+        vgcoeffs[i] = (void*) (gcoeffs[i]);
+    }
+
+    tensorFuserCRT (y, tupSize, ppGCRTRq, totm, peArr, sizeOfPE, vgcoeffs, qs);
+
+#ifdef DEBUG_MODE
+    for(j = 0; j < totm; j++)
+  {
+      if(y[j]<0)
+      {
+          printf("tensorGCRTRq\n");
+      }
+  }
+#endif
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    gcrqTime = tsAdd(gcrqTime, tsSubtract(t1,s1));
+#endif
+}
+void tensorGCRTC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, complex_t** gcoeffs)
+{
+#ifdef STATS
+    gccCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+#ifdef DEBUG_MODE
+    printf("\n\nEntered tensorGCRTC\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\n[", totm, sizeOfPE);
+    hDim_t j;
+    for(j = 0; j < totm; j++) {
+        printf("(%f,%f),", (y[j]).real, (y[j]).imag);
+    }
+    printf("]\n[");
+    for(j = 0; j < sizeOfPE; j++) {
+        printf("(%" PRId32 ",%" PRId16 "),", peArr[j].prime, peArr[j].exponent);
+    }
+    printf("]\n");
+#endif
+    void** vgcoeffs = (void**)malloc(sizeOfPE*sizeof(void*));
+    hDim_t i;
+    for(i = 0; i < sizeOfPE; i++)
+    {
+        vgcoeffs[i] = (void*) (gcoeffs[i]);
+    }
+
+    tensorFuserCRT (y, tupSize, ppGCRTC, totm, peArr, sizeOfPE, vgcoeffs, (hInt_t*)0);
+
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    gccTime = tsAdd(gccTime, tsSubtract(t1,s1));
+#endif
+}
+void tensorGInvCRTRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t** gcoeffs, hInt_t* qs)
+{
+#ifdef STATS
+    gicrqCtr++;
+#endif
+    tensorGCRTRq (tupSize, y, totm, peArr, sizeOfPE, gcoeffs, qs); //output is already shifted
+}
+void tensorGInvCRTC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, complex_t** gcoeffs)
+{
+#ifdef STATS
+    giccCtr++;
+#endif
+    tensorGCRTC (tupSize, y, totm, peArr, sizeOfPE, gcoeffs);
+}
+

@@ -0,0 +1,305 @@+#include "tensorTypes.h"
+
+hDim_t ipow(hDim_t base, hShort_t exp)
+{
+#ifdef DEBUG_MODE
+    ASSERT(exp >= 0);
+#endif
+	hDim_t result = 1;
+    while (exp)
+    {
+        if (exp & 1)
+        {
+            result *= base;
+        }
+        exp >>= 1;
+        base *= base;
+    }
+    return result;
+}
+
+complex_t cmplxpow(complex_t base, hShort_t exp)
+{
+	complex_t result = (complex_t){1,0};
+    while (exp)
+    {
+        if (exp & 1)
+        {
+            CMPLX_IMUL(result,base);
+        }
+        exp >>= 1;
+        CMPLX_IMUL(base,base);
+    }
+    return result;
+}
+
+hInt_t qpow(hInt_t base, hShort_t exp, hInt_t q)
+{
+	hInt_t result = 1;
+    while (exp)
+    {
+        if (exp & 1)
+        {
+            result = (result*base)%q;
+        }
+        exp >>= 1;
+        base = (base*base)%q;
+    }
+    return result;
+}
+
+// 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;
+	}
+	ASSERT (a==1);  // if this one fails, then b is not invertible mod a
+
+	// this actually returns EITHER the reciprocal OR reciprocal + fieldSize
+	hInt_t res = lasty + fieldSize;
+#ifdef DEBUG_MODE
+	ASSERT (0);
+	ASSERT ((res >= 0) && (res < fieldSize + fieldSize));
+	hInt_t test = res * b % fieldSize;
+	ASSERT (test == 1);
+#endif
+	return res;
+
+}
+
+//for square transforms
+void tensorFuser (void* y, hShort_t tupSize, funcPtr 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;
+        (*f) (y, tupSize, pe, lts, rts, qs);
+        rts  *= dim;
+    }
+}
+
+void tensorFuserCRT (void* y, hShort_t tupSize, crtFuncPtr f, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, void** ru, hInt_t* q)
+{
+    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;
+        (*f) (y, tupSize, lts, rts, pe, ru[i], q);
+        rts  *= dim;
+    }
+}
+
+struct  timespec  tsSubtract (struct  timespec  time1, struct  timespec  time2)
+{    /* Local variables. */
+    struct  timespec  result ;
+
+/* Subtract the second time from the first. */
+
+    if ((time1.tv_sec < time2.tv_sec) ||
+        ((time1.tv_sec == time2.tv_sec) &&
+         (time1.tv_nsec <= time2.tv_nsec))) {		/* TIME1 <= TIME2? */
+        result.tv_sec = result.tv_nsec = 0 ;
+    } else {						/* TIME1 > TIME2 */
+        result.tv_sec = time1.tv_sec - time2.tv_sec ;
+        if (time1.tv_nsec < time2.tv_nsec) {
+            result.tv_nsec = time1.tv_nsec + 1000000000L - time2.tv_nsec ;
+            result.tv_sec-- ;				/* Borrow a second. */
+        } else {
+            result.tv_nsec = time1.tv_nsec - time2.tv_nsec ;
+        }
+    }
+
+    return (result) ;
+}
+
+struct  timespec  tsAdd (struct  timespec  time1, struct  timespec  time2)
+{    /* Local variables. */
+    struct  timespec  result ;
+
+/* Add the two times together. */
+
+    result.tv_sec = time1.tv_sec + time2.tv_sec ;
+    result.tv_nsec = time1.tv_nsec + time2.tv_nsec ;
+    if (result.tv_nsec >= 1000000000L) {		/* Carry? */
+        result.tv_sec++ ;  result.tv_nsec = result.tv_nsec - 1000000000L ;
+    }
+
+    return (result) ;
+}
+
+const  char  *tsShow (struct  timespec  binaryTime, bool  inLocal, const  char  *format)
+{    /* Local variables. */
+    struct  tm  calendarTime ;
+#define  MAX_TIMES  4
+    static  char  asciiTime[MAX_TIMES][64] ;
+    static  int  current = 0 ;
+
+/* Convert the TIMESPEC to calendar time: year, month, day, etc. */
+
+#ifdef VXWORKS
+    if (inLocal)
+        localtime_r ((time_t *) &binaryTime.tv_sec, &calendarTime) ;
+    else
+        gmtime_r ((time_t *) &binaryTime.tv_sec, &calendarTime) ;
+#else
+    if (inLocal)
+        calendarTime = *(localtime ((time_t *) &binaryTime.tv_sec)) ;
+    else
+        calendarTime = *(gmtime ((time_t *) &binaryTime.tv_sec)) ;
+#endif
+
+/* Format the time in ASCII. */
+
+    current = (current + 1) % MAX_TIMES ;
+
+    if (format == NULL) {
+        strftime (asciiTime[current], 64, "%Y-%j-%H:%M:%S", &calendarTime) ;
+        sprintf (asciiTime[current] + strlen (asciiTime[current]),
+                 ".%06ld", (binaryTime.tv_nsec % 1000000000L) / 1000L) ;
+    } else {
+        strftime (asciiTime[current], 64, format, &calendarTime) ;
+        sprintf (asciiTime[current] + strlen (asciiTime[current]),
+                 ".%06ld", (binaryTime.tv_nsec % 1000000000L) / 1000L) ;
+    }
+
+    return (asciiTime[current]);
+}
+
+
+
+const char* timeformat = "%M:%S";
+
+void getStats() { 
+
+#ifdef STATS
+    struct timespec total;
+    printf("CRT Stats:\n");
+    printf("CRT_Rq times: Real:%s\tMono:%s\tProc:%s\tThread:%s\n", tsShow(crttime1, false, timeformat),tsShow(crttime2, false, timeformat),tsShow(crttime3, false, timeformat),tsShow(crttime4, false, timeformat));
+    printf("CTR_Rq: %d\t%s\t%d\t%s\n", crtRqCtr, tsShow(crttime1, false, timeformat), crtInvRqCtr, tsShow(crtInvRqTime, false, timeformat));
+    printf("CTR_C: %d\t%s\t%d\t%s\n", crtCCtr, tsShow(crtCTime, false, timeformat), crtInvCCtr, tsShow(crtInvCTime, false, timeformat));
+    
+    printf("\nG Stats:\n");
+    printf("GPow_R: %d\t%s\t%d\t%s\n", gprCtr, tsShow(gprTime, false, timeformat), giprCtr, tsShow(giprTime, false, timeformat));
+    printf("GPow_Rq: %d\t%s\t%d\t%s\n", gprqCtr, tsShow(gprqTime, false, timeformat), giprqCtr, tsShow(giprqTime, false, timeformat));
+    printf("GPow_C: %d\t%s\t%d\t%s\n", gpcCtr, tsShow(gpcTime, false, timeformat), gipcCtr, tsShow(gipcTime, false, timeformat));
+    printf("GDec_R: %d\t%s\t%d\t%s\n", gdrCtr, tsShow(gdrTime, false, timeformat), gidrCtr, tsShow(gidrTime, false, timeformat));
+    printf("GDec_Rq: %d\t%s\t%d\t%s\n", gdrqCtr, tsShow(gdrqTime, false, timeformat), gidrqCtr, tsShow(gidrqTime, false, timeformat));
+    printf("GCRT_Rq: %d\t%d\t%s\n", gcrqCtr, gicrqCtr, tsShow(gcrqTime, false, timeformat));
+    printf("GCRT_C: %d\t%d\t%s\n", gccCtr, giccCtr, tsShow(gccTime, false, timeformat));
+
+    printf("\nL Stats:\n");
+    printf("L_R: %d\t%s\t%d\t%s\n", lrCtr, tsShow(lrTime, false, timeformat), lirCtr, tsShow(lirTime, false, timeformat));
+    printf("L_Rq: %d\t%s\t%d\t%s\n", lrqCtr, tsShow(lrqTime, false, timeformat), lirqCtr, tsShow(lirqTime, false, timeformat));
+    printf("L_D: %d\t%s\t%d\t%s\n", ldCtr, tsShow(ldTime, false, timeformat), lidCtr, tsShow(lidTime, false, timeformat));
+    printf("L_C: %d\t%s\t%d\t%s\n", lcCtr, tsShow(lcTime, false, timeformat), licCtr, tsShow(licTime, false, timeformat));
+
+    printf("\nNorm Stats:\n");
+    printf("NormSq_R %d\t%s\n", normrCtr, tsShow(normrTime, false, timeformat));
+
+    printf("\nBasic Stats:\n");
+    printf("Mul: %d\t%s\n", mulCtr, tsShow(mulTime, false, timeformat));
+    printf("Add: %d\t%s\n", addCtr, tsShow(addTime, false, timeformat));
+
+    total = tsAdd(norrTime, tsAdd(crttime1, tsAdd(crtInvRqTime, tsAdd(crtCTime, tsAdd(crtInvCTime, tsAdd(gprTime, tsAdd(giprTime, tsAdd(gdrTime, tsAdd(gidrTime, tsAdd(gprqTime, tsAdd(giprqTime, tsAdd(gdrqTime, tsAdd(gidrqTime, tsAdd(gcrqTime, tsAdd(gccTime, tsAdd(lrTime, tsAdd(lirTime, tsAdd(lrqTime, tsAdd(lirqTime, tsAdd(ldTime, tsAdd(lidTime, tsAdd(lcTime, tsAdd(licTime, tsAdd(mulTime,addTime))))))))))))))))))))))));
+
+    printf("\nTotal C Time: %s\n\n", tsShow(total, false, timeformat));
+
+    crtRqCtr = 0;
+    crtInvRqCtr = 0;
+    crtCCtr = 0;
+    crtInvCCtr = 0;
+
+    gprCtr = 0;
+    gpcCtr = 0;
+    gprqCtr = 0;
+    gdrCtr = 0;
+    gdrqCtr = 0;
+    giprCtr = 0;
+    gipcCtr = 0;
+    giprqCtr = 0;
+    gidrCtr = 0;
+    gidrqCtr = 0;
+    gcrqCtr = 0;
+    gccCtr = 0;
+    gicrqCtr = 0;
+    giccCtr = 0;
+
+    lrqCtr = 0;
+    lrCtr = 0;
+    ldCtr = 0;
+    lcCtr = 0;
+    lirqCtr = 0;
+    lirCtr = 0;
+    lidCtr = 0;
+    licCtr = 0;
+
+    normrCtr = 0;
+
+    mulCtr = 0;
+    addCtr = 0;
+
+    mulTime = (struct timespec){0,0};
+    addTime = (struct timespec){0,0};
+
+    lrqTime = (struct timespec){0,0};
+    lrTime = (struct timespec){0,0};
+    ldTime = (struct timespec){0,0};
+    lcTime = (struct timespec){0,0};
+    lirqTime = (struct timespec){0,0};
+    lirTime = (struct timespec){0,0};
+    lidTime = (struct timespec){0,0};
+    licTime = (struct timespec){0,0};
+
+    normrTime = (struct timespec){0,0};
+
+    gprTime = (struct timespec){0,0};
+    gpcTime = (struct timespec){0,0};
+    gprqTime = (struct timespec){0,0};
+    gdrTime = (struct timespec){0,0};
+    gdrqTime = (struct timespec){0,0};
+    giprTime = (struct timespec){0,0};
+    gipcTime = (struct timespec){0,0};
+    giprqTime = (struct timespec){0,0};
+    gidrTime = (struct timespec){0,0};
+    gidrqTime = (struct timespec){0,0};
+    gcrqTime = (struct timespec){0,0};
+    gccTime = (struct timespec){0,0};
+
+    crttime1 = (struct timespec){0,0};
+    crttime2 = (struct timespec){0,0};
+    crttime3 = (struct timespec){0,0};
+    crttime4 = (struct timespec){0,0};
+
+    crtInvRqTime = (struct timespec){0,0};
+    crtCTime = (struct timespec){0,0};
+    crtInvCTime = (struct timespec){0,0};
+#endif
+    fflush(stdout);
+}
+
+

@@ -0,0 +1,434 @@+#include "tensorTypes.h"
+
+void lpRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hInt_t q) {
+  hDim_t blockOffset;
+  hDim_t modOffset;
+  int 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 idx = tmp2 + modOffset + rts;
+      for (i = 1; i < p-1; ++i) {
+        hInt_t temp = y[(idx-rts)*tupSize] + y[idx*tupSize];
+        if (temp >= q) y[idx*tupSize]=temp-q;
+        else y[idx*tupSize] = temp;
+        idx += rts;
+      }
+    }
+  }
+}
+
+void lpR (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) {
+  hDim_t blockOffset;
+  hDim_t modOffset;
+  int 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 idx = tmp2 + modOffset + rts;
+      for (i = 1; i < p-1; ++i) {
+        y[idx*tupSize] += y[(idx-rts)*tupSize];
+        idx += rts;
+      }
+    }
+  }
+}
+
+void lpDouble (double* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) {
+  hDim_t blockOffset;
+  hDim_t modOffset;
+  int 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 idx = tmp2 + modOffset + rts;
+      for (i = 1; i < p-1; ++i) {
+        y[idx*tupSize] += y[(idx-rts)*tupSize];
+        idx += rts;
+      }
+    }
+  }
+}
+
+void lpC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) {
+  hDim_t blockOffset;
+  hDim_t modOffset;
+  int 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 idx = tmp2 + modOffset + rts;
+      for (i = 1; i < p-1; ++i) {
+        CMPLX_IADD (y[idx*tupSize], y[(idx-rts)*tupSize]);
+        idx += rts;
+      }
+    }
+  }
+}
+
+void lpInvRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hInt_t q) {
+  hDim_t blockOffset;
+  hDim_t modOffset;
+  int 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;
+      hDim_t idx = tensorOffset + (p-2) * rts;
+      for (i = p-2; i != 0; --i) {
+        hInt_t temp = y[idx*tupSize] - y[(idx-rts)*tupSize] + q;
+        if (temp >= q) y[idx*tupSize]=temp-q;
+        else y[idx*tupSize] = temp;
+        idx -= rts;
+      }
+    }
+  }
+}
+
+void lpInvR (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) {
+  hDim_t blockOffset;
+  hDim_t modOffset;
+  int 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;
+      hDim_t idx = tensorOffset + (p-2) * rts;
+      for (i = p-2; i != 0; --i) {
+        y[idx*tupSize] -= y[(idx-rts)*tupSize] ;
+        idx -= rts;
+      }
+    }
+  }
+}
+
+void lpInvDouble (double* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) {
+  hDim_t blockOffset;
+  hDim_t modOffset;
+  int 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;
+      hDim_t idx = tensorOffset + (p-2) * rts;
+      for (i = p-2; i != 0; --i) {
+        y[idx*tupSize] -= y[(idx-rts)*tupSize] ;
+        idx -= rts;
+      }
+    }
+  }
+}
+
+void lpInvC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p) {
+  hDim_t blockOffset;
+  hDim_t modOffset;
+  int 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;
+      hDim_t idx = tensorOffset + (p-2) * rts;
+      for (i = p-2; i != 0; --i) {
+        CMPLX_ISUB (y[idx*tupSize], y[(idx-rts)*tupSize]);
+        idx -= rts;
+      }
+    }
+  }
+}
+
+void ppLRq (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs) {
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if(p == 2) return;
+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+      lpRq (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p, qs[tupIdx]);
+    }
+}
+
+void ppLR (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs) {
+#ifdef DEBUG_MODE
+  ASSERT (q==0);
+#endif
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if(p == 2) return;
+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+      lpR (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);
+    }
+}
+
+void ppLDouble (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs) {
+#ifdef DEBUG_MODE
+  ASSERT (q==0);
+#endif
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if(p == 2) return;
+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+      lpDouble (((double*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);
+    }
+}
+
+void ppLC (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs) {
+#ifdef DEBUG_MODE
+  ASSERT (q==0);
+#endif
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if(p == 2) return;
+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+      lpC (((complex_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);
+    }
+}
+
+
+void ppLInvRq (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs) {
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if(p == 2) return;
+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+      lpInvRq (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p, qs[tupIdx]);
+    }
+}
+
+void ppLInvR (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs) {
+#ifdef DEBUG_MODE
+  ASSERT (q==0);
+#endif
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if(p == 2) return;
+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+      lpInvR (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);
+    }
+}
+
+void ppLInvDouble (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs) {
+#ifdef DEBUG_MODE
+  ASSERT (q==0);
+#endif
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if(p == 2) return;
+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+      lpInvDouble (((double*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);
+    }
+}
+
+void ppLInvC (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs) {
+#ifdef DEBUG_MODE
+  ASSERT (q==0);
+#endif
+    hDim_t p = pe.prime;
+    hShort_t e = pe.exponent;
+    if(p == 2) return;
+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
+      lpInvC (((complex_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);
+    }
+}
+
+#ifdef STATS
+int lrqCtr = 0;
+int lrCtr = 0;
+int ldCtr = 0;
+int lcCtr = 0;
+int lirqCtr = 0;
+int lirCtr = 0;
+int lidCtr = 0;
+int licCtr = 0;
+
+struct timespec lrqTime = {0,0};
+struct timespec lrTime = {0,0};
+struct timespec ldTime = {0,0};
+struct timespec lcTime = {0,0};
+struct timespec lirqTime = {0,0};
+struct timespec lirTime = {0,0};
+struct timespec lidTime = {0,0};
+struct timespec licTime = {0,0};
+#endif
+
+
+void tensorLRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs) {
+#ifdef STATS
+    lrqCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+#ifdef DEBUG_MODE
+    hDim_t i;
+    printf("\n\nEntered tensorLRq\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\tq=%" PRId64 "\n[", totm, sizeOfPE, q);
+    /*for(i = 0; i < totm; i++) {
+        printf("%" PRId64 ",", y[i]);
+    }*/
+    printf("]\n[");
+    for(i = 0; i < sizeOfPE; i++) {
+        printf("(%" PRId32 ",%" PRId16 "),", peArr[i].prime, peArr[i].exponent);
+    }
+    printf("]\n");
+#endif
+  tensorFuser (y, tupSize, ppLRq, totm, peArr, sizeOfPE, qs); // don't need to shift here
+#ifdef DEBUG_MODE
+  for(i = 0; i < totm*tupSize; i++) {
+      if(y[i]<0) {
+          printf("tensorLRq\n");
+      }
+  }
+#endif
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    lrqTime = tsAdd(lrqTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorLR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) {
+#ifdef STATS
+    lrCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+#ifdef DEBUG_MODE
+    printf("\n\nEntered tensorLR\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\n[", totm, sizeOfPE);
+    hDim_t i;
+    for(i = 0; i < totm; i++) {
+        printf("%" PRId64 ",", y[i]);
+    }
+    printf("]\n[");
+    for(i = 0; i < sizeOfPE; i++) {
+        printf("(%" PRId32 ",%" PRId16 "),", peArr[i].prime, peArr[i].exponent);
+    }
+    printf("]\n");
+#endif
+  tensorFuser (y, tupSize, ppLR, totm, peArr, sizeOfPE, (hInt_t*)0);
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    lrTime = tsAdd(lrTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorLDouble (hShort_t tupSize, double* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) {
+#ifdef STATS
+    ldCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+#ifdef DEBUG_MODE
+    printf("\n\nEntered tensorLDouble\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\n[", totm, sizeOfPE);
+    hDim_t i;
+    for(i = 0; i < totm; i++) {
+        printf("%f,", y[i]);
+    }
+    printf("]\n[");
+    for(i = 0; i < sizeOfPE; i++) {
+        printf("(%" PRId32 ",%" PRId16 "),", peArr[i].prime, peArr[i].exponent);
+    }
+    printf("]\n");
+#endif
+  tensorFuser (y, tupSize, ppLDouble, totm, peArr, sizeOfPE, (hInt_t*)0);
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    ldTime = tsAdd(ldTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorLC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) {
+#ifdef STATS
+    lcCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+#ifdef DEBUG_MODE
+    printf("\n\nEntered tensorLC\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\n[", totm, sizeOfPE);
+    hDim_t i;
+    for(i = 0; i < totm; i++) {
+        printf("(%f,%f),", y[i].real, y[i].imag);
+    }
+    printf("]\n[");
+    for(i = 0; i < sizeOfPE; i++) {
+        printf("(%" PRId32 ",%" PRId16 "),", peArr[i].prime, peArr[i].exponent);
+    }
+    printf("]\n");
+#endif
+  tensorFuser (y, tupSize, ppLC, totm, peArr, sizeOfPE, (hInt_t*)0);
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    lcTime = tsAdd(lcTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorLInvRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs) {
+#ifdef STATS
+    lirqCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppLInvRq, totm, peArr, sizeOfPE, qs);
+#ifdef DEBUG_MODE
+  hDim_t i;
+  for(i = 0; i < totm*tupSize; i++)
+  {
+      if(y[i]<0)
+      {
+          printf("tensorLInvRq\n");
+      }
+  }
+#endif
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    lirqTime = tsAdd(lirqTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorLInvR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) {
+#ifdef STATS
+    lirCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppLInvR, totm, peArr, sizeOfPE, (hInt_t*)0);
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    lirTime = tsAdd(lirTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorLInvDouble (hShort_t tupSize, double* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) {
+#ifdef STATS
+    lidCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppLInvDouble, totm, peArr, sizeOfPE, (hInt_t*)0);
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    lidTime = tsAdd(lidTime, tsSubtract(t1,s1));
+#endif
+}
+
+void tensorLInvC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) {
+#ifdef STATS
+    licCtr++;
+    struct timespec s1,t1;
+    clock_gettime(CLOCK_REALTIME, &s1);
+#endif
+  tensorFuser (y, tupSize, ppLInvC, totm, peArr, sizeOfPE, (hInt_t*)0);
+#ifdef STATS
+    clock_gettime(CLOCK_REALTIME, &t1);
+    licTime = tsAdd(licTime, tsSubtract(t1,s1));
+#endif
+}
+

@@ -0,0 +1,86 @@+#include "tensorTypes.h"++#ifdef STATS+int normrCtr = 0;+struct timespec normrTime = {0,0};+#endif++void pNormSqR (hInt_t* 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;+      hInt_t 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;+      }+    }+  }+}++void ppNormSqR (void* y, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs) {+#ifdef DEBUG_MODE+  ASSERT (q==0);+#endif+    hDim_t p = pe.prime;+    hShort_t e = pe.exponent;+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {+      pNormSqR (((hInt_t*)y)+tupIdx, tupSize, lts*ipow(p,e-1), rts, p);+    }+}++void tensorNormSqR (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE) {+#ifdef STATS+    normrCtr++;+    struct timespec s1,t1;+    clock_gettime(CLOCK_REALTIME, &s1);+#endif+#ifdef DEBUG_MODE+    printf("\n\nEntered tensorNormSqR\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\n[", totm, sizeOfPE);+    hDim_t i;+    for(i = 0; i < totm; i++) {+        printf("%" PRId64 ",", y[i]);+    }+    printf("]\n[");+    for(i = 0; i < sizeOfPE; i++) {+        printf("(%" PRId32 ",%" PRId16 "),", peArr[i].prime, peArr[i].exponent);+    }+    printf("]\n");+#endif++  hInt_t* tempSpace = (hInt_t*)malloc(totm*tupSize*sizeof(hInt_t));+  for(hDim_t i = 0; i < totm*tupSize; i++) {+    tempSpace[i]=y[i];+  }++  tensorFuser(y, tupSize, ppNormSqR, 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);++#ifdef STATS+    clock_gettime(CLOCK_REALTIME, &t1);+    lrTime = tsAdd(normrTime, tsSubtract(t1,s1));+#endif+}

@@ -0,0 +1,74 @@++#include <math.h>+#include <stdlib.h>+#include "tensorTypes.h"++// this function takes *inverse* RUs, so no negation is needed on the indexing+// I had been negating the ru-idx, but this was causing a *negative* mod, resulting in a hard-to-find bug+void primeD (double *y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hDim_t rustride, complex_t* ruinv)+{+	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 * ruinv[((row*col) % p)*rustride*tupSize].real * y[(tensorOffset+rts*(col-1))*tupSize];+        }+        for(col = (p>>1)+1; col <= p-1; col++)+        {+          acc += 2 * ruinv[((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 (void *y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, void *ruinv, hInt_t* qs)+{+    hDim_t p = pe.prime;+    hDim_t e = pe.exponent;+#ifdef DEBUG_MODE+    ASSERT(e != 0);+#endif+    hDim_t mprime = ipow(p,e-1);+    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {+      primeD (((double*)y)+tupIdx, tupSize, lts*mprime, rts, p, mprime, ((complex_t*)ruinv)+tupIdx);+    }+}++//the contents of y will be destroyed, but should be initialized in Haskell-land to independent Guassians over the reals+void tensorGaussianDec (hShort_t tupSize, double *y, hDim_t totm, PrimeExponent *peArr, hShort_t sizeOfPE, complex_t** ruinv)+{+  void** ruinvs = (void**)malloc(sizeOfPE*sizeof(void*));+  hShort_t i;+  for(i = 0; i < sizeOfPE; i++)+  {+      ruinvs[i] = (void*) (ruinv[i]);+  }+    +	tensorFuserCRT (y, tupSize, ppD, totm, peArr, sizeOfPE, ruinvs, (hInt_t*)0);+	+	free(ruinvs);+}

@@ -0,0 +1,231 @@+
+#ifndef TENSORTYPES_H_
+#define TENSORTYPES_H_
+
+
+// remove next line for more efficient code
+//#define DEBUG_MODE
+
+
+#include <stdbool.h>
+#include <inttypes.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+
+
+#define ASSERT(EXP) { \
+	if (!(EXP)) { \
+		fprintf (stderr, "Assertion in file '%s' line %d : " #EXP "  is false\n", __FILE__, __LINE__); \
+		exit(-1); \
+	} \
+}
+
+
+//timers and counters
+#ifdef STATS
+extern int crtRqCtr;
+extern int crtInvRqCtr;
+extern int crtCCtr;
+extern int crtInvCCtr;
+
+extern int gprCtr;
+extern int gpcCtr;
+extern int gprqCtr;
+extern int gdrCtr;
+extern int gdrqCtr;
+extern int giprCtr;
+extern int gipcCtr;
+extern int giprqCtr;
+extern int gidrCtr;
+extern int gidrqCtr;
+extern int gcrqCtr;
+extern int gccCtr;
+extern int gicrqCtr;
+extern int giccCtr;
+
+extern int lrqCtr;
+extern int lrCtr;
+extern int ldCtr;
+extern int lcCtr;
+extern int lirqCtr;
+extern int lirCtr;
+extern int lidCtr;
+extern int licCtr;
+
+extern int normrCtr;
+
+extern int mulCtr;
+extern struct timespec mulTime;
+extern int addCtr;
+extern struct timespec addTime;
+
+extern struct timespec lrqTime;
+extern struct timespec lrTime;
+extern struct timespec ldTime;
+extern struct timespec lcTime;
+extern struct timespec lirqTime;
+extern struct timespec lirTime;
+extern struct timespec lidTime;
+extern struct timespec licTime;
+
+extern struct timespec normrTime;
+
+extern struct timespec crttime1;
+extern struct timespec crttime2;
+extern struct timespec crttime3;
+extern struct timespec crttime4;
+extern struct timespec crtInvRqTime;
+extern struct timespec crtCTime;
+extern struct timespec crtInvCTime;
+
+extern struct timespec gprTime;
+extern struct timespec gpcTime;
+extern struct timespec gprqTime;
+extern struct timespec gdrTime;
+extern struct timespec gdrqTime;
+extern struct timespec giprTime;
+extern struct timespec gipcTime;
+extern struct timespec giprqTime;
+extern struct timespec gidrTime;
+extern struct timespec gidrqTime;
+extern struct timespec gcrqTime;
+extern struct timespec gccTime;
+#endif
+
+typedef int64_t hInt_t ;
+typedef int32_t hDim_t ;
+typedef int16_t hShort_t ;
+typedef int8_t hByte_t ;
+
+typedef struct
+{
+	hDim_t prime;
+	hShort_t exponent;
+}  PrimeExponent;
+
+
+typedef struct
+{
+	double real;
+	double imag;
+} complex_t;
+
+#define CMPLX_ADD(a,b)  ((complex_t){((a).real + (b).real), ((a).imag + (b).imag)})
+#define CMPLX_ADD3(a,b,c)  ((complex_t){((a).real + (b).real + (c).real), ((a).imag + (b).imag + (c).imag)})
+#define CMPLX_ADD4(a,b,c,d)  ((complex_t){((a).real + (b).real + (c).real + (d).real), ((a).imag + (b).imag + (c).imag + (d).imag)})
+#define CMPLX_ADD5(a,b,c,d,e)  ((complex_t){((a).real + (b).real + (c).real + (d).real + (e).real), ((a).imag + (b).imag + (c).imag + (d).imag + (e).imag)})
+#define CMPLX_ADD6(a,b,c,d,e,f)  ((complex_t){((a).real + (b).real + (c).real + (d).real + (e).real + (f).real), ((a).imag + (b).imag + (c).imag + (d).imag + (e).imag + (f).imag)})
+#define CMPLX_ADD7(a,b,c,d,e,f,g)  ((complex_t){((a).real + (b).real + (c).real + (d).real + (e).real + (f).real + (g).real), ((a).imag + (b).imag + (c).imag + (d).imag + (e).imag + (f).imag + (g).imag)})
+
+#define CMPLX_SUB(a,b)  ((complex_t){((a).real - (b).real), ((a).imag - (b).imag)})
+#define CMPLX_MUL(a,b)  ((complex_t){((a).real*(b).real - (a).imag*(b).imag), \
+	                                  (a).real*(b).imag + (a).imag*(b).real})
+#define CMPLX_DIV(a,b)  ((complex_t){((a).real*(b).real + (a).imag*(b).imag)/((b).real*(b).real+(b).imag*(b).imag), \
+                                     ((a).imag*(b).real - (a).real*(b).imag)/((b).real*(b).real+(b).imag*(b).imag)})
+
+// 'inside' operators
+#define CMPLX_IADD(a,b)  { (a).real += (b).real;  (a).imag += (b).imag; }
+#define CMPLX_ISUB(a,b)  { (a).real -= (b).real;  (a).imag -= (b).imag; }
+#define CMPLX_IMUL(a,b)  { double temp = ((a).real*(b).real - (a).imag*(b).imag); \
+	                       (a).imag = ((a).real*(b).imag + (a).imag*(b).real); \
+	                       (a).real = temp; }
+
+// calculates base ** exp
+hDim_t ipow(hDim_t base, hShort_t exp);
+complex_t cmplxpow(complex_t base, hShort_t exp);
+hInt_t qpow(hInt_t base, hShort_t exp, hInt_t q);
+
+hInt_t reciprocal (hInt_t a, hInt_t b);
+
+struct  timespec  tsSubtract (struct  timespec  time1, struct  timespec  time2);
+struct  timespec  tsAdd (struct  timespec  time1, struct  timespec  time2);
+const  char  *tsShow (struct  timespec  binaryTime, bool  inLocal, const  char  *format);
+
+void getStats();
+
+void mulRq (hShort_t tupSize, hInt_t* a, hInt_t* b, hDim_t totm, hInt_t* qs);
+//void mulMq (hInt_t* a, const hInt_t* b, const hDim_t totm, const hByte_t logr, const hInt_t k, const hInt_t q);
+void mulC (hShort_t tupSize, complex_t* a, complex_t* b, hDim_t totm);
+
+void addR (hShort_t tupSize, hInt_t* a, hInt_t* b, hDim_t totm);
+void addRq (hShort_t tupSize, hInt_t* a, const hInt_t* b, const hDim_t totm, const hInt_t* qs);
+//void addMq (hInt_t* a, const hInt_t* b, const hDim_t totm, const hByte_t logr, const hInt_t k, const hInt_t q);
+void addC (hShort_t tupSize, complex_t* a, complex_t* b, hDim_t totm);
+void addD (hShort_t tupSize, double* a, double* b, hDim_t totm);
+
+typedef void (*funcPtr) (void* outputVec, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* qs);
+void tensorFuser (void* y, hShort_t tupSize, funcPtr f, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs);
+
+typedef void (*normFuncPtr) (void* outputVec, hShort_t tupSize, PrimeExponent pe, hDim_t lts, hDim_t rts, hInt_t* output, hInt_t* qs);
+void tensorFuserNorm (void* y, hShort_t tupSize, normFuncPtr f, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* output, hInt_t q);
+
+typedef void (*crtFuncPtr) (void* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, void* ru, hInt_t* q);
+void tensorFuserCRT (void* y, hShort_t tupSize, crtFuncPtr f, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, void** ru, hInt_t* q);
+
+void tensorGPowR (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+void tensorGPowRq (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs);
+
+void tensorGPowC (hShort_t tupSize, complex_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+void tensorGDecR (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+void tensorGDecRq (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs);
+
+void tensorGInvPowR (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+void tensorGInvPowRq (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs);
+
+void tensorGInvPowC (hShort_t tupSize, complex_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+void tensorGInvDecR (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+void tensorGInvDecRq (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs);
+
+void tensorGCRTRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t** gcoeffs, hInt_t* qs);
+
+void tensorGInvCRTRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t** gcoeffs, hInt_t* qs);
+
+void tensorGCRTC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, complex_t** gcoeffs);
+
+void tensorGInvCRTC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, complex_t** gcoeffs);
+
+
+
+void tensorLRq (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs);
+
+void tensorLR (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+void tensorLDouble (hShort_t tupSize, double* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+void tensorLC (hShort_t tupSize, complex_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+void tensorLInvRq (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t* qs);
+
+void tensorLInvR (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+void tensorLInvDouble (hShort_t tupSize, double* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+void tensorLInvC (hShort_t tupSize, complex_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+
+
+void tensorCRTRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t** ru, hInt_t* qs);
+
+void tensorCRTC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, complex_t** ru);
+
+void tensorCRTInvRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t** ru, hInt_t* minv, hInt_t* qs);
+
+void tensorCRTInvC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, complex_t** ru, complex_t* minv);
+
+
+
+void tensorGaussianDec (hShort_t tupSize, double* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, complex_t** ru);
+
+void tensorNormSqR (hShort_t tupSize, hInt_t* x, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE);
+
+
+#endif /* TENSORTYPES_H_ */
+

@@ -0,0 +1,258 @@+{-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts,+             FlexibleInstances, GADTs, MultiParamTypeClasses,+             NoImplicitPrelude, PolyKinds, RebindableSyntax,+             RoleAnnotations, ScopedTypeVariables, StandaloneDeriving,+             TypeFamilies, TypeOperators, UndecidableInstances #-}++-- | A pure, repa-based implementation of the Tensor interface.++module Crypto.Lol.Cyclotomic.Tensor.RepaTensor+( RT ) where++import Crypto.Lol.Cyclotomic.Tensor                      as T+import Crypto.Lol.Cyclotomic.Tensor.RepaTensor.CRT+import Crypto.Lol.Cyclotomic.Tensor.RepaTensor.Dec+import Crypto.Lol.Cyclotomic.Tensor.RepaTensor.Extension+import Crypto.Lol.Cyclotomic.Tensor.RepaTensor.GL+import Crypto.Lol.Cyclotomic.Tensor.RepaTensor.RTCommon  as RT+import Crypto.Lol.LatticePrelude                         as LP hiding+                                                                ((!!))+import Crypto.Lol.Reflects+import Crypto.Lol.Types.FiniteField                      as FF+import Crypto.Lol.Types.IZipVector++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        hiding (arr)+import Control.DeepSeq      (NFData (rnf))+import Control.Monad.Random+import Data.Coerce+import Data.Constraint      hiding ((***))+import Data.Foldable        as F+import Data.Maybe+import Data.Traversable     as T+import Data.Vector          as V hiding (force)+import Data.Vector.Unboxed  as U hiding (force)+import Test.QuickCheck++-- | An implementation of 'Tensor' backed by repa.+data RT (m :: Factored) r where+  RT :: Unbox r => !(Arr m r) -> RT m r+  ZV :: IZipVector m r -> RT m r++deriving instance Show r => Show (RT m r)++instance Eq r => Eq (RT m r) where+  (ZV a) == (ZV b) = a == b+  (RT a) == (RT b) = a == b+  a@(RT _) == b = a == toRT b+  a == b@(RT _) = toRT a == b+  {-# INLINABLE (==) #-}++zvToArr :: Unbox r => IZipVector m r -> Arr m r+zvToArr v = let vec = convert $ unIZipVector v+            in Arr $ fromUnboxed (Z :. U.length vec) vec++-- converts to RT constructor+toRT :: Unbox r => RT m r -> RT m r+toRT v@(RT _) = v+toRT (ZV v) = RT $ zvToArr v++toZV :: Fact m => RT m r -> RT m r+toZV (RT (Arr v)) = ZV $ fromMaybe (error "toZV: internal error") $+                    iZipVector $ convert $ toUnboxed v+toZV v@(ZV _) = v++{-# INLINABLE wrap #-}+wrap :: Unbox r => (Arr l r -> Arr m r) -> RT l r -> RT m r+wrap f (RT v) = RT $ f v+wrap f (ZV v) = RT $ f $ zvToArr v++{-# INLINABLE wrapM #-}+wrapM :: (Unbox r, Monad mon) => (Arr l r -> mon (Arr m r))+         -> RT l r -> mon (RT m r)+wrapM f (RT v) = RT <$> f v+wrapM f (ZV v) = RT <$> f (zvToArr v)++instance Tensor RT where++  type TElt RT r = (Unbox r, Elt 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 = RT . scalarPow'++  l = wrap fL+  lInv = wrap fLInv++  mulGPow = wrap fGPow+  mulGDec = wrap fGDec++  divGPow = wrapM fGInvPow+  divGDec = wrapM fGInvDec++  crtFuncs = (,,,,) <$>+             ((RT .) <$> scalarCRT') <*>+             (wrap <$> mulGCRT') <*>+             (wrap <$> divGCRT') <*>+             (wrap <$> fCRT) <*>+             (wrap <$> fCRTInv)++  tGaussianDec = fmap RT . tGaussianDec'++  gSqNormDec (RT e) = gSqNormDec' e+  gSqNormDec e = gSqNormDec $ toRT e++  twacePowDec = wrap twacePowDec'++  embedPow = wrap embedPow'+  embedDec = wrap embedDec'++  crtExtFuncs = (,) <$> (wrap <$> twaceCRT') <*> (wrap <$> embedCRT')++  coeffs = wrapM coeffs'++  powBasisPow = (RT <$>) <$> powBasisPow'++  crtSetDec = (RT <$>) <$> crtSetDec'++  fmapT f (RT v) = RT $ (coerce $ force . RT.map f) v+  fmapT f v@(ZV _) = fmapT f $ toRT v++  -- Repa arrays don't have mapM, so apply to underlying Unboxed+  -- vector instead+  fmapTM f (RT (Arr arr)) = (RT . Arr . fromUnboxed (extent arr)) <$>+                            U.mapM f (toUnboxed arr)+  fmapTM f v = fmapTM f $ toRT v++  zipWithT f (RT (Arr a1)) (RT (Arr a2)) = RT $ Arr $ force $ RT.zipWith f a1 a2+  zipWithT f v1 v2 = zipWithT f (toRT v1) (toRT v2)++  unzipTElt (RT (Arr arr)) = (RT . Arr . fromUnboxed (extent arr)) ***+                             (RT . Arr . fromUnboxed (extent arr)) $+                             U.unzip $ toUnboxed arr+  unzipTElt v = unzipTElt $ toRT v++  unzipT v@(RT _) = unzipT $ toZV v+  unzipT (ZV v) = ZV *** ZV $ unzipIZV v++  {-# 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 #-}+  {-# INLINABLE crtExtFuncs #-}+  {-# INLINABLE coeffs #-}+  {-# INLINABLE powBasisPow #-}+  {-# INLINABLE crtSetDec #-}+  {-# INLINABLE fmapT #-}+  {-# INLINABLE fmapTM #-}+  {-# INLINABLE zipWithT #-}+  {-# INLINABLE unzipTElt #-}+  {-# INLINABLE unzipT #-}+++---------- Category-theoretic instances ----------++instance Fact m => Functor (RT m) where+  -- Functor instance is implied by Applicative+  fmap f x = pure f <*> x++instance Fact m => Applicative (RT m) where+  pure = ZV . pure++  -- RT can never hold an a -> b+  (ZV f) <*> (ZV a) = ZV (f <*> a)+  f@(ZV _) <*> v@(RT _) = f <*> toZV v++instance Fact m => Foldable (RT m) where+  -- Foldable instance is implied by Traversable+  foldMap = foldMapDefault++instance Fact m => Traversable (RT m) where+  traverse f r@(RT _) = T.traverse f $ toZV r+  traverse f (ZV v) = ZV <$> T.traverse f v+++---------- Numeric Prelude instances ----------++--CJP: Additive, Ring are not necessary when we use zipWithT+--EAC: But we need an Additive instance for the Module instance++instance (Unbox r, Additive (Arr m r)) => Additive.C (RT m r) where+  zero = RT zero++  (RT a) + (RT b) = RT $ a + b+  a + b = toRT a + toRT b++  (RT a) - (RT b) = RT $ a - b+  a - b = toRT a - toRT b++  negate (RT a) = RT $ negate a+  negate a = negate $ toRT a++  {-# INLINABLE (+) #-}+  {-# INLINABLE (-) #-}+  {-# INLINABLE zero #-}+  {-# INLINABLE negate #-}++{-+instance (Unbox r, Ring (Arr m r)) => Ring.C (RT m r) where+  (RT a) * (RT b) = RT $ a * b+  a * b = toRT a * toRT b++  fromInteger = RT . fromInteger+  {-# INLINABLE (*) #-}+  {-# INLINABLE fromInteger #-}+-}++instance (ZeroTestable (Arr m r), ZeroTestable (IZipVector m r))+    => ZeroTestable.C (RT m r) where+  isZero (RT a) = isZero a+  isZero (ZV v) = isZero v+  {-# INLINABLE isZero #-}++instance (GFCtx fp d, Fact m, Additive (RT m fp))+    => Module.C (GF fp d) (RT m fp) where++  r *> v = case v of+    RT (Arr arr) -> RT $ Arr $ RT.fromList (extent arr) $ unCoeffs $ r *> Coeffs $ RT.toList arr+    ZV zv -> ZV $ fromJust $ iZipVector $ V.fromList $ unCoeffs $ r *> Coeffs $ V.toList $ unIZipVector zv++---------- Miscellaneous instances ----------++instance (Unbox r, Random (Arr m r)) => Random (RT m r) where+  random = runRand $ RT <$> (liftRand random)++  randomR = error "randomR nonsensical for RT"++instance (Unbox r, Arbitrary (Arr m r)) => Arbitrary (RT m r) where+  arbitrary = RT <$> arbitrary++instance (NFData r) => NFData (RT m r) where+  rnf (RT v) = rnf v+  rnf (ZV v) = rnf v

@@ -0,0 +1,215 @@+{-# LANGUAGE ConstraintKinds, FlexibleContexts, GADTs, NoImplicitPrelude,+             ScopedTypeVariables #-}++-- | Functions to support the chinese remainder transform on Repa arrays++module Crypto.Lol.Cyclotomic.Tensor.RepaTensor.CRT+( scalarCRT'+, fCRT, fCRTInv+, mulGCRT', divGCRT'+, gCRT, gInvCRT+) where++import Crypto.Lol.CRTrans+import Crypto.Lol.Cyclotomic.Tensor+import Crypto.Lol.Cyclotomic.Tensor.RepaTensor.GL+import Crypto.Lol.Cyclotomic.Tensor.RepaTensor.RTCommon as RT+import Crypto.Lol.LatticePrelude                        as LP++import Control.Applicative+import Data.Coerce+import Data.Singletons.Prelude++-- | Embeds a scalar into the CRT basis (when it exists).+scalarCRT' :: forall m r . (Fact m, CRTrans r, Unbox r)+              => Maybe (r -> Arr m r)+{-# INLINABLE scalarCRT' #-}+scalarCRT'+  = let pps = proxy ppsFact (Proxy::Proxy m)+        sz = Z :. totientPPs pps+    in pure $ Arr . force . fromFunction sz . const++-- | Multiply by @g_m@ in the CRT basis (when it exists).+mulGCRT' :: forall m r . (Fact m, CRTrans r, Unbox r, Elt r)+            => Maybe (Arr m r -> Arr  m r)+{-# INLINABLE mulGCRT' #-}+mulGCRT' = (coerce (\x -> force . RT.zipWith (*) x) `asTypeOf` asTypeOf) <$> gCRT++-- | Divide by @g@ in the CRT basis (when it exists).+divGCRT' :: (Fact m, CRTrans r, IntegralDomain r, ZeroTestable r,+             Unbox r, Elt r) => Maybe (Arr m r -> Arr m r)+{-# INLINABLE divGCRT' #-}+divGCRT' =  (coerce (\x -> force . RT.zipWith (*) x) `asTypeOf` asTypeOf) <$> gInvCRT++-- | The representation of @g@ in the CRT basis (when it exists).+gCRT :: (Fact m, CRTrans r, Unbox r, Elt r) => Maybe (Arr m r)+{-# INLINABLE gCRT #-}+gCRT = fCRT <*> pure (fGPow $ scalarPow' LP.one)++-- | The representation of @g^{ -1 }@ in the CRT basis (when it exists).+gInvCRT:: (Fact m, CRTrans r, IntegralDomain r,+           ZeroTestable r, Unbox r, Elt r)+          => Maybe (Arr m r)+{-# INLINABLE gInvCRT #-}+gInvCRT = fCRT <*> fGInvPow (scalarPow' LP.one)++fCRT, fCRTInv ::+  forall m r . (Fact m, CRTrans r, Unbox r, Elt r)+  => Maybe (Arr m r -> Arr m r)++{-# INLINABLE fCRT #-}+{-# INLINABLE fCRTInv #-}++-- | The chinese remainder transform.+-- Exists if and only if crt exists for all prime powers.+fCRT = evalM $ fTensor ppCRT++-- divide by mhat after doing crtInv'+-- | The inverse chinese remainder transform.+-- Exists if and only if crt exists for all prime powers.+fCRTInv = do -- in Maybe+  (_, mhatInv) :: (CRTInfo r) <- proxyT crtInfoFact (Proxy :: Proxy m)+  let totm = proxy totientFact (Proxy :: Proxy m)+      divMhat = trans totm $ RT.map (*mhatInv)+  evalM $ (divMhat .*) <$> fTensor ppCRTInv'++ppDFT, ppDFTInv', ppCRT, ppCRTInv' ::+  forall pp r . (PPow pp, CRTrans r, Unbox r, Elt r)+  => TaggedT pp Maybe (Trans r)++{-# INLINABLE ppDFT #-}+{-# INLINABLE ppDFTInv' #-}+{-# INLINABLE ppCRT #-}+{-# INLINABLE ppCRTInv' #-}++ppDFT = case (sing :: SPrimePower pp) of+  (SPP (STuple2 sp SO)) -> tagT $ withWitnessT pDFT sp+  spp@(SPP (STuple2 sp (SS se'))) ->+    tagT $ do+      let spp' = SPP (STuple2 sp se')+      pp'dft <- withWitnessT ppDFT spp'+      pptwid <- withWitnessT (ppTwid False) spp+      pdft <- withWitnessT pDFT sp+      return $ (pp'dft @* Id (dim pdft)) .* pptwid .* (Id (dim pp'dft) @* pdft)++ppDFTInv' = case (sing :: SPrimePower pp) of+  (SPP (STuple2 sp SO)) -> tagT $ withWitnessT pDFTInv' sp+  spp@(SPP (STuple2 sp (SS se'))) ->+    tagT $ do+      let spp' = SPP (STuple2 sp se')+      pp'dftInv' <- withWitnessT ppDFTInv' spp'+      pptwidInv <- withWitnessT (ppTwid True) spp+      pdftInv' <- withWitnessT pDFTInv' sp+      return $ (Id (dim pp'dftInv') @* pdftInv') .* pptwidInv .*+                 (pp'dftInv' @* Id (dim pdftInv'))++ppCRT = case (sing :: SPrimePower pp) of+  (SPP (STuple2 sp SO)) -> tagT $ withWitnessT pCRT sp+  spp@(SPP (STuple2 sp (SS se'))) ->+    tagT $ do+      let spp' = SPP (STuple2 sp se')+      pp'dft <- withWitnessT ppDFT spp'+      pptwid <- withWitnessT (ppTwidHat False) spp+      pcrt <- withWitnessT pCRT sp+      return $+        (pp'dft @* Id (dim pcrt)) .* pptwid .*+        -- save some work when p=2+        (if dim pcrt > 1 then Id (dim pp'dft) @* pcrt else Id (dim pp'dft))++ppCRTInv' = case (sing :: SPrimePower pp) of+  (SPP (STuple2 sp SO)) -> tagT $ withWitnessT pCRTInv' sp+  spp@(SPP (STuple2 sp (SS se'))) ->+    tagT $ do+      let spp' = SPP (STuple2 sp se')+      pp'dftInv' <- withWitnessT ppDFTInv' spp'+      pptwidInv <- withWitnessT (ppTwidHat True) spp+      pcrtInv' <- withWitnessT pCRTInv' sp+      return $+        (Id (dim pp'dftInv') @* pcrtInv') .* pptwidInv .*+        (pp'dftInv' @* Id (dim pcrtInv'))++butterfly :: (Additive r, Unbox r, Elt r) => Trans r+butterfly = trans 2 $ \arr ->+            fromFunction (extent arr) $+                             \(sh:.j) -> case j of+                                           0 -> arr ! (sh:.0) ++                                                arr ! (sh:.1)+                                           1 -> arr ! (sh:.0) -+                                                arr ! (sh:.1)++-- DFT_p, CRT_p, scaled DFT_p^{ -1 } and CRT_p^{ -1 }+pDFT, pDFTInv', pCRT, pCRTInv' ::+  forall p r . (Prim p, CRTrans r, Unbox r, Elt r)+  => TaggedT p Maybe (Trans r)++{-# INLINABLE pDFT #-}+{-# INLINABLE pDFTInv' #-}+{-# INLINABLE pCRT #-}+{-# INLINABLE pCRTInv' #-}++pDFT = let pval = proxy valuePrime (Proxy::Proxy p)+       in if pval == 2+          then return butterfly+          else do (omegaPPow, _) <- crtInfoPrime+                  return $ trans pval $ mulMat $ force $+                         fromFunction (Z :. pval :. pval)+                                          (\(Z:.i:.j) -> omegaPPow (i*j))++pDFTInv' = let pval = proxy valuePrime (Proxy::Proxy p)+           in if pval == 2+              then return butterfly+              else do (omegaPPow, _) <- crtInfoPrime+                      return $ trans pval $ mulMat $ force $+                             fromFunction (Z :. pval :. pval)+                                              (\(Z:.i:.j) -> omegaPPow (-i*j))++pCRT = let pval = proxy valuePrime (Proxy::Proxy p)+       in if pval == 2+          then return $ Id 1+          else do (omegaPPow, _) <- crtInfoPrime+                  return $ trans (pval-1) $ mulMat $ force $+                         fromFunction (Z :. pval-1 :. pval-1)+                                          (\(Z:.i:.j) -> omegaPPow ((i+1)*j))++-- crt_p * this = \hat{p}*I, for all prime p.+pCRTInv' =+  let pval = proxy valuePrime (Proxy::Proxy p)+  in if pval == 2 then return $ Id 1+     else do+       (omegaPPow, _) <- crtInfoPrime+       return $ trans (pval-1) $  mulMat $ force $+              fromFunction (Z :. pval-1 :. pval-1)+                               (\(Z:.i:.j) -> omegaPPow (negate i*(j+1)) -+                                              omegaPPow (j+1))++-- twiddle factors for DFT_pp and CRT_pp decompositions+ppTwid, ppTwidHat ::+  forall pp r . (PPow pp, CRTrans r, Unbox r, Elt r)+  => Bool -> TaggedT pp Maybe (Trans r)++{-# INLINABLE ppTwid #-}+{-# INLINABLE ppTwidHat #-}++ppTwid inv =+  let pp@(p,e) = proxy ppPPow (Proxy :: Proxy pp)+      ppval = valuePP pp+  in do+    (omegaPPPow, _) <- crtInfoPPow+    return $ trans ppval $ mulDiag $ force $+                           fromFunction (Z :. ppval)+                           (\(Z:.i) -> let (iq,ir) = i `divMod` p+                                           pow = (if inv then negate else id)+                                                 ir * digitRev (p,e-1) iq+                                       in omegaPPPow pow)++ppTwidHat inv =+  let pp@(p,e) = proxy ppPPow (Proxy :: Proxy pp)+      pptot = totientPP pp+  in do+    (omegaPPPow, _) <- crtInfoPPow+    return $ trans pptot $ mulDiag $ force $+                           fromFunction (Z :. pptot)+                           (\(Z:.i) -> let (iq,ir) = i `divMod` (p-1)+                                           pow = (if inv then negate else id)+                                                 (ir+1) * digitRev (p,e-1) iq+                                       in omegaPPPow pow)