{-# LANGUAGE TypeFamilies #-} {- | Functions for generating arrays of complex unit numbers, needed as twiddle factors and for chirps. -} module Numeric.BLAS.Rank1.Unit where import qualified Numeric.FFTW.Extra.Shape as CShape import Numeric.FFTW.Extra.Utility (takeShape, takeSubshape, assembleInterval1) import Numeric.FFTW.Extra.NumberTheory (sqrtToPowerOfTwo) import qualified Numeric.FFTW.Rank1 as Rank1 import qualified Numeric.BLAS.Matrix.RowMajor as Matrix import qualified Numeric.BLAS.Vector.Chunk as Chunk import qualified Numeric.BLAS.Vector.Symbolic as VectorSymb import qualified Numeric.BLAS.Vector.Slice as VectorSlice import qualified Numeric.BLAS.Vector as Vector import qualified Numeric.BLAS.Subobject.Layout as Layout import qualified Numeric.BLAS.Slice as Slice import qualified Numeric.Netlib.Class as Class -- import Numeric.Netlib.Modifier (Conjugation(NonConjugated)) import qualified Data.Array.Comfort.Storable as Array import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable (Array) import qualified Data.List as List import qualified Data.Complex as Complex import Data.Complex (Complex) import Control.Applicative (liftA2) {- $setup >>> :set -XTypeFamilies >>> import Numeric.BLAS.Rank1.Unit >>> import Test.Numeric.FFTW.Extra.Utility >>> import DocTest.NumberType_.Numeric.FFTW.Extra.Rank1.Convolution (complex_) >>> import qualified Numeric.FFTW.Rank1 as Rank1 >>> import qualified Numeric.BLAS.Matrix.RowMajor as Matrix >>> import qualified Numeric.BLAS.Vector.Slice as VectorSlice >>> >>> import qualified Data.Array.Comfort.Shape as Shape >>> >>> import qualified Test.QuickCheck as QC -} _twiddleVectorChecked :: (Shape.Indexed shape, Shape.Index shape ~ n, Integral n, Class.Real a) => a -> n -> shape -> n -> Array shape (Complex a) _twiddleVectorChecked c n shape i = Array.sample shape (\j -> if i*j < n then Complex.cis (c * fromIntegral (i*j)) else error "twiddle index too big") {- ToDo: Fast cis computation. Store an array 'precomp' of 128 sin values from 0 to pi/4. Compute precomp ! intPart x * let y = fracPart x in (1-y^2/2) :+ (y-y^3/3). This would be precise enough for Float. We might require the user to choose this optimization via a Cabal flag. *Numeric.BLAS.Matrix.RowMajor> pi/(4*128) :: Float 6.1359233e-3 *Numeric.BLAS.Matrix.RowMajor> 2^^(-23::Int) :: Float 1.1920929e-7 *Numeric.BLAS.Matrix.RowMajor> (pi/(4*128))^(3::Int)/6::Float 3.8502463e-8 That is, for those small numbers, we can substitute sin x by x in Float. *Numeric.BLAS.Matrix.RowMajor> 2^^(-53::Int) :: Double 1.1102230246251565e-16 *Numeric.BLAS.Matrix.RowMajor> (pi/4/512)^(5::Int)/120::Double 7.078127084107336e-17 For those numbers we would approximate sin x = x-x^3/6 in Double. Maybe 'cis' is already doing this? We could setup precomputed data in arrays of fixed-point values. This is independent from floating point format. Maybe one quadrant of sine in Arrayn 1024 Int64 and another quadrant of sine in Array 225 Int64. 225 = 3*3*5*5. The maximum number of odd prime factors found in results of ceilingEfficient. -} twiddleVectorWrap :: (Shape.Indexed shape, Shape.Index shape ~ n, Integral n, Class.Real a) => a -> n -> shape -> n -> Array shape (Complex a) twiddleVectorWrap c n shape i = Array.sample shape (\j -> Complex.cis (c * fromIntegral (mod (i*j) n))) divUp :: (Integral n) => n -> n -> n divUp n m = - div (-n) m {- | prop> :{ QC.forAll genVectorShape $ \shape -> let n = Shape.zeroBasedSize shape in let c = 2*pi / fromIntegral n in QC.forAll (chooseLogarithmic n) $ \m i -> VectorSlice.toVector (twiddleVectorCascaded c n shape m i) =~= (twiddleVectorWrap c n shape i &:: complex_) :} -} twiddleVectorCascaded :: (Integral n, Shape.ZeroBased n ~ shape, Class.Real a) => a -> n -> shape -> n -> n -> VectorSlice.T Layout.Subvector shape (Complex a) twiddleVectorCascaded c nn shape@(Shape.ZeroBased n) k i = let m = divUp n k in takeSubshape shape $ Matrix.tensorProduct (twiddleVectorWrap c nn (Shape.ZeroBased m) (k*i)) (twiddleVectorWrap c nn (Shape.ZeroBased k) i) {- twiddleVectorAuto_ :: (Integral n, Shape.ZeroBased n ~ shape, Class.Real a) => a -> n -> shape -> n -> VectorSlice.T Layout.Subvector shape (Complex a) twiddleVectorAuto_ c nn shape@(Shape.ZeroBased n) i = if n < 16 then VectorSlice.fromVector $ twiddleVectorWrap c nn shape i else twiddleVectorCascaded c nn shape (sqrtToPowerOfTwo n) i -} twiddleVectorAuto :: (Integral n, Shape.ZeroBased n ~ shape, Class.Real a) => a -> n -> shape -> n -> VectorSlice.T Layout.Subvector shape (Complex a) twiddleVectorAuto c nn shape@(Shape.ZeroBased n) i = if n < 16 then twiddleVectorCascaded c nn shape (sqrtToPowerOfTwo n) i -- 32 was the fastest in my measurements else twiddleVectorAutoLoop 32 c nn shape i twiddleVectorAutoLoop :: (Integral n, Shape.ZeroBased n ~ shape, Class.Real a) => n -> a -> n -> shape -> n -> VectorSlice.T Layout.Subvector shape (Complex a) twiddleVectorAutoLoop d c nn = let go shape@(Shape.ZeroBased n) i = if n < 16 then twiddleVectorCascaded c nn shape (sqrtToPowerOfTwo n) i else takeSubshape shape $ Matrix.tensorProduct (go (Shape.ZeroBased $ divUp n d) (d*i)) (twiddleVectorWrap c nn (Shape.ZeroBased d) i) in go {- 49*1024 - 223*225 == 1 49/225 - 223/1024 == 1/(225*1024) *Main Key> Key.minimum (\(a,b,c) -> max a(max b c)) $ map (\k -> (k,mod(49*k)225, mod(223*k)1024)) [1..(225*1024::Integer)] (23,2,9) *Main> mapM_ print $ take 10 $ filter (\(a,b,c) -> max b c < 100) $ map (\k -> (k,mod(49*k)225, mod(223*k)1024)) [1..(100::Integer)] Embed binary data in a library: https://hackage.haskell.org/package/file-embed/ https://hackage.haskell.org/package/data-embed/ -} {-# SPECIALIZE precomputed1024 :: Array (Shape.ZeroBased Int) Float #-} {-# SPECIALIZE precomputed1024 :: Array (Shape.ZeroBased Int) Double #-} {-# SPECIALIZE precomputed225 :: Array (Shape.ZeroBased Int) Float #-} {-# SPECIALIZE precomputed225 :: Array (Shape.ZeroBased Int) Double #-} precomputed1024, precomputed225 :: (Class.Real a) => Array (Shape.ZeroBased Int) a precomputed1024 = let n = 1024 in let c = pi/2 / fromIntegral n in Array.sample (Shape.ZeroBased n) (\j -> sin (c * fromIntegral j)) precomputed225 = let n = 225 in let c = pi/2 / fromIntegral n in Array.sample (Shape.ZeroBased n) (\j -> sin (c * fromIntegral j)) twiddleVector :: (Shape.Indexed shape, Shape.Index shape ~ n, Integral n, Class.Real a) => a -> shape -> n -> Array shape (Complex a) twiddleVector c shape i = Array.sample shape (\j -> Complex.cis (c * fromIntegral (i*j))) {- | prop> :{ QC.forAll (QC.choose (-2*pi,2*pi)) $ \c -> QC.forAll genVectorShape $ \shape -> let n = Shape.zeroBasedSize shape in QC.forAll (chooseLogarithmic n) $ \m -> QC.forAll (QC.choose (1, div n m + 1)) $ \i -> VectorSlice.toVector (twiddleVectorFast c shape m i) =~= (twiddleVector c shape i &:: complex_) :} -} -- ToDo: we need VectorSlice.mapSourceShape in order to hide (shape, shape) structure twiddleVectorFast :: (Integral n, Shape.ZeroBased n ~ shape, Class.Real a) => a -> shape -> n -> n -> VectorSlice.T Layout.Slice shape (Complex a) twiddleVectorFast c shape@(Shape.ZeroBased n) k i = takeShape shape $ Matrix.tensorProduct (twiddleVector c (Shape.ZeroBased $ divUp n k) (k*i)) (twiddleVector c (Shape.ZeroBased k) i) {- | prop> :{ QC.forAll (QC.choose (-pi,pi)) $ \c -> QC.forAll genVectorShape $ \(Shape.ZeroBased n) -> QC.forAll (QC.choose (-2*n,n)) $ \start -> let shape = Shape.Shifted start n in QC.forAll (chooseLogarithmic n) $ \m -> QC.forAll (QC.choose (1, div n m + 1)) $ \i -> VectorSlice.toVector (twiddleVectorShiftedFast c shape m i) =~= (twiddleVector c shape i &:: complex_) :} -} twiddleVectorShiftedFast :: (Integral n, Shape.Shifted n ~ shape, Class.Real a) => a -> shape -> n -> n -> VectorSlice.T Layout.Slice shape (Complex a) twiddleVectorShiftedFast c shape@(Shape.Shifted from n) k i = takeShape shape $ Matrix.tensorProduct (twiddleVector c (Shape.ZeroBased $ divUp n k) (k*i)) (twiddleVector c (Shape.Shifted from k) i) {- | The product i*j must not exceed n. This is true as found for twiddle matrices in the regular Cooley-Tukey transform. It is not true when abusing twiddleMatrix to compute the Fourier matrix. In those cases we need to wrap-around indices with (mod (i*j) n) for more floating point accuracy. -} twiddleMatrix :: (Shape.Indexed height, Shape.Index height ~ n, Shape.Indexed width, Shape.Index width ~ n, Integral n, Class.Real a) => a -> (height, width) -> Array (height, width) (Complex a) twiddleMatrix c shape = Array.sample shape (\(j,k) -> Complex.cis (c * fromIntegral (j*k))) -- ToDo: use twiddleVectorFast twiddleMatrixRows :: (Shape.IntegralRange height, Shape.Index height ~ n, Shape.Indexed width, Shape.Index width ~ n, Integral n, Class.Real a) => a -> n -> (height, width) -> [Array width (Complex a)] twiddleMatrixRows c step (height, width) = map (twiddleVector c width) $ take (Shape.size height) $ iterate (step+) $ Shape.shiftedOffset $ Shape.toShifted height -- iterate (step+) (Shape.intervalFrom $ Shape.toInterval height) {- | prop> :{ QC.forAll (QC.choose (-2*pi,2*pi)) $ \c -> QC.forAll genVectorShape $ \(Shape.ZeroBased nk) -> QC.forAll (chooseLogarithmic nk) $ \k -> let n = div nk k in let width = Shape.ZeroBased k in QC.forAll (QC.choose (-2*n,n)) $ \start -> let height = Shape.intervalFromShifted $ Shape.Shifted start n in QC.forAll (chooseLogarithmic n) $ \m -> let tm = twiddleMatrix c (height, width) &:: complex_ in tm =~= (Matrix.reshapeHeight height $ Matrix.fromRows width $ take n $ twiddleMatrixRowsCascaded c m (height,width)) :} -} twiddleMatrixRowsCascaded :: (Shape.IntegralRange height, Shape.Index height ~ n, Shape.Indexed width, Shape.Index width ~ n, Eq width, Integral n, Class.Real a) => a -> n -> (height, width) -> [Array width (Complex a)] twiddleMatrixRowsCascaded c k (height, width) = let Shape.Shifted from n = Shape.toShifted height in liftA2 Vector.mul (twiddleMatrixRows c k (Shape.ZeroBased (divUp n k), width)) (twiddleMatrixRows c 1 (Shape.Shifted from k, width)) twiddleMatrixRowsAuto :: (Shape.IntegralRange height, Shape.Index height ~ n, Shape.Indexed width, Shape.Index width ~ n, Eq width, Integral n, Class.Real a) => a -> (height, width) -> [Array width (Complex a)] twiddleMatrixRowsAuto c (height, width) = let n = Shape.shiftedSize (Shape.toShifted height) in if n < 16 then twiddleMatrixRows c 1 (height,width) else twiddleMatrixRowsCascaded c (sqrtToPowerOfTwo n) (height,width) {- (k+j)^2/2 - k^2/2 - j^2/2 = k*j -} {- | prop> :{ QC.forAll genVectorShape $ \width sign -> let n = Shape.zeroBasedSize width in let cu = 2*pi / fromIntegral n in let c = case sign of Rank1.Forward -> cu; Rank1.Backward -> -cu in QC.forAll (chooseLogarithmic n) $ \m -> let height = Shape.ZeroBased m in (twiddleMatrix c (height, width) &:: complex_) =~= (Matrix.reshapeHeight height $ Matrix.fromRows width $ twiddleMatrixViaChirp sign n (height,width)) :} -} twiddleMatrixViaChirp :: (Shape.ZeroBased n ~ height, Shape.ZeroBased n ~ width, Integral n, Class.Real a) => Rank1.Sign -> n -> (height, width) -> [Array width (Complex a)] twiddleMatrixViaChirp sign n (height, width) = let h = Shape.zeroBasedSize height in let w = Shape.zeroBasedSize width in let hw = h+w-1 in let chirp = nonNegChirp sign n hw in let chirpPrefix = takeShape width chirp in map (\(j,c) -> Vector.scale (Complex.conjugate c) $ VectorSlice.mul (VectorSymb.conjugate chirpPrefix) $ VectorSlice.sliceVector (Slice.reshape width . Slice.sub j w) chirp) $ Array.toAssociations $ Vector.take h chirp {- | prop> :{ QC.forAll genVectorShape $ \widthZB fromW sign -> let n = Shape.zeroBasedSize widthZB in let width = Shape.Shifted fromW n in let cu = 2*pi / fromIntegral n in let c = case sign of Rank1.Forward -> cu; Rank1.Backward -> -cu in QC.forAll (chooseLogarithmic n) $ \m fromH -> let height = Shape.Shifted fromH m in (twiddleMatrix c (height, width) &:: complex_) =~= (Matrix.reshapeHeight height $ Matrix.fromRows width $ twiddleMatrixViaChirpShifted sign n (height,width)) :} -} {- Numeric.FFTW.Extra.Rank1.Transform:333: *** Failed! Falsified (after 76 tests and 5 shrinks): ZeroBased {zeroBasedSize = 1} 47 Forward 1 74 StorableArray.fromList (Shifted {shiftedOffset = 74, shiftedSize = 1},Shifted {shiftedOffset = 47, shiftedSize = 1}) [0.999999 :+ 1.423504e-3] StorableArray.fromList (Shifted {shiftedOffset = 74, shiftedSize = 1},Shifted {shiftedOffset = 47, shiftedSize = 1}) [1.0 :+ 0.0] ((74,47),(-1.013279e-6) :+ 1.423504e-3) chirpForShape already contains 'mod'. Thus we cannot easily avoid this inaccuracy. -} {- This is more general than twiddleMatrixViaChirp but also requires to compute three independent chirps. -} -- ToDo: alternative Shape.toShifted :: (Shape.IntegralRange sh) => sh -> Shape.Shifted (Shape.Index sh) twiddleMatrixViaChirpShifted :: (Shape.Shifted n ~ height, Shape.Shifted n ~ width, Integral n, Class.Real a) => Rank1.Sign -> n -> (height, width) -> [Array width (Complex a)] twiddleMatrixViaChirpShifted sign n (height, width@(Shape.Shifted _ w)) = let shapeUnshift = Shape.ZeroBased . Shape.shiftedSize in let chirp = Array.mapShape shapeUnshift $ chirpForShape sign n $ CShape.convolve height width in let chirpH = Array.mapShape shapeUnshift $ chirpForShape sign n height in let chirpW = chirpForShape sign n width in map (\(j,c) -> -- ToDo: merge mul and scale Vector.scale (Complex.conjugate c) $ VectorSlice.mul (VectorSymb.conjugate chirpW) $ -- ToDo: Slice.subShifted would be nice to have -- or even: Slice.subShape start shape VectorSlice.sliceVector (Slice.reshape width . Slice.sub j w) chirp) $ Array.toAssociations chirpH -- ToDo: the chirp is 2*n-periodic, could use Spectrum.Symmetric shape chirpForShape :: (Shape.Indexed sh, Shape.Index sh ~ n, Integral n, Class.Real a) => Rank1.Sign -> n -> sh -> Array sh (Complex a) chirpForShape sign n shape = let cu = pi / fromIntegral n in let c = case sign of Rank1.Forward -> cu; Rank1.Backward -> -cu in Array.sample shape (\k -> Complex.cis (c * fromIntegral (mod (k*k) (2*n)))) nonNegChirp :: (Integral n, Class.Real a) => Rank1.Sign -> n -> n -> Array (Shape.ZeroBased n) (Complex a) nonNegChirp sign n m = chirpForShape sign n (Shape.ZeroBased m) {- | prop> :{ QC.forAll genVectorShape $ \shape (QC.Positive n) sign -> chirpForShape sign n shape =~= (chirpForShapeScaled sign n 1 shape &:: complex_) :} -} chirpForShapeScaled :: (Shape.Indexed sh, Shape.Index sh ~ n, Integral n, Class.Real a) => Rank1.Sign -> n -> n -> sh -> Array sh (Complex a) chirpForShapeScaled sign n j shape = let cu = pi / fromIntegral n in let c = case sign of Rank1.Forward -> cu; Rank1.Backward -> -cu in Array.sample shape (\k -> Complex.cis (c * fromIntegral (mod (j*k*k) (2*n)))) {- | prop> :{ QC.forAll genVectorShape $ \(Shape.ZeroBased m) (QC.Positive n) sign -> QC.forAll (chooseLogarithmic m) $ \k -> nonNegChirp sign n m =~= (nonNegChirpFast sign n m k &:: complex_) :} -} nonNegChirpFast :: (Integral n, Class.Real a) => Rank1.Sign -> n -> n -> n -> Array (Shape.ZeroBased n) (Complex a) nonNegChirpFast sign n m k = let cu = 2*pi*fromIntegral k / fromIntegral n in let c = case sign of Rank1.Forward -> cu; Rank1.Backward -> -cu in let mk = divUp m k in VectorSlice.toVector $ takeShape (Shape.ZeroBased m) $ -- Vector.mul (twiddleMatrix c (Shape.ZeroBased mk, Shape.ZeroBased k)) $ Vector.mul (Matrix.reshapeHeight (Shape.ZeroBased mk) $ Matrix.fromRows (Shape.ZeroBased k) $ List.genericTake mk $ twiddleMatrixRowsAuto c (Shape.ZeroBased mk, Shape.ZeroBased k)) $ Matrix.tensorProduct (chirpForShapeScaled sign n (k*k) (Shape.ZeroBased mk)) (chirpForShape sign n (Shape.ZeroBased k)) {- (i+j*k)^2 = i^2 + 2*i*j*k + (j*k)^2 k*(i+j)^2 = k*i^2 + 2*i*j*k + k*j^2 (i+j*k)^2 = i^2 + (j*k)^2 + k*(i+j)^2 - k*i^2 - k*j^2 = (1-k)*i^2 + (k-1)*k*j^2 + k*(i+j)^2 -} {- | prop> :{ QC.forAll genVectorShape $ \(Shape.ZeroBased m) (QC.Positive n) sign -> QC.forAll (chooseLogarithmic m) $ \k -> nonNegChirp sign n m =~= (nonNegChirpChirpy sign n m k &:: complex_) :} -} nonNegChirpChirpy :: (Integral n, Class.Real a) => Rank1.Sign -> n -> n -> n -> Array (Shape.ZeroBased n) (Complex a) nonNegChirpChirpy sign n m k = let mk = divUp m k in let longChirp = chirpForShapeScaled sign n k (Shape.ZeroBased (k+mk-1)) in VectorSlice.toVector $ takeShape (Shape.ZeroBased m) $ Vector.mul (Matrix.reshapeHeight (Shape.ZeroBased mk) $ Matrix.fromRows (Shape.ZeroBased k) $ map (\j -> VectorSlice.sliceVector (Slice.sub j k) longChirp) $ List.genericTake mk $ iterate (1+) 0) $ Matrix.tensorProduct (chirpForShapeScaled sign n ((k-1)*k) (Shape.ZeroBased mk)) (chirpForShapeScaled sign n (1-k) (Shape.ZeroBased k)) {- | prop> :{ QC.forAll genVectorShape $ \(Shape.ZeroBased m) from (QC.Positive n) sign -> QC.forAll (chooseLogarithmic m) $ \k -> let shape = Shape.Shifted from m in chirpForShape sign n shape =~= (chirpForShapeChirpy sign n shape k &:: complex_) :} -} -- ToDo: This subsumes nonNegChirpChirpy and has not much overhead. chirpForShapeChirpy :: (Shape.IntegralRange sh, Shape.Index sh ~ n, Integral n) => (Class.Real a) => Rank1.Sign -> n -> sh -> n -> Array sh (Complex a) chirpForShapeChirpy sign n shape k = let Shape.Shifted from m = Shape.toShifted shape in let mk = divUp m k in let longLen = k+mk-1 in let longChirp = Array.reshape (Shape.ZeroBased longLen) $ chirpForShapeScaled sign n k (Shape.Shifted from longLen) in VectorSlice.toVector $ takeShape shape $ Vector.mul (Array.reshape (Shape.ZeroBased mk, Shape.Shifted from k) $ Matrix.fromRows (Shape.ZeroBased k) $ map (\j -> VectorSlice.sliceVector (Slice.sub j k) longChirp) $ List.genericTake mk $ iterate (1+) 0) $ Matrix.tensorProduct (chirpForShapeScaled sign n ((k-1)*k) (Shape.ZeroBased mk)) (chirpForShapeScaled sign n (1-k) (Shape.Shifted from k)) appendRevToInterval1 :: (VectorSlice.C v, VectorSlice.C w, Integral n, Class.Floating a) => v (Shape.ZeroBased n) a -> w (Shape.ZeroBased n) a -> Array (Shape.Interval n) a appendRevToInterval1 x y = Array.reshape (assembleInterval1 (VectorSlice.shape x) (VectorSlice.shape y)) $ Chunk.toVector (Chunk.reversed x <> Chunk.simple y) {- | prop> :{ QC.forAll genVectorShape $ \shape sign -> let n = Shape.zeroBasedSize shape in forChoose (0,n) $ \to -> forChoose (-n,to) $ \from -> let rng = Shape.Interval from to in VectorSlice.toVector (interval2FromChirp rng (nonNegChirp sign n (n+1))) =~= chirpForShape sign n rng &:: complex_ :} -} interval2FromChirp :: (Integral n, Class.Floating a) => Shape.Interval n -> Array (Shape.ZeroBased n) a -> VectorSlice.T Layout.Slice (Shape.Interval n) a interval2FromChirp rng@(Shape.Interval from to) chirp = if from >= 0 then VectorSlice.sliceVector (Slice.subInterval rng) $ Array.mapShape Shape.intervalFromZeroBased chirp else VectorSlice.fromVector $ appendRevToInterval1 (takeShape (Shape.ZeroBased (1-from)) chirp) (VectorSlice.sliceVector (Slice.sub 1 to) chirp) {- | prop> :{ QC.forAll genVectorShape $ \shape sign from (QC.NonNegative offset) -> let n = Shape.zeroBasedSize shape in let to = from+n+offset in let rng = Shape.Interval from to in interval3FromChirp rng (nonNegChirp sign n (n+1)) =~= chirpForShape sign n rng &:: complex_ :} -} interval3FromChirp :: (Integral n, Class.Floating a) => Shape.Interval n -> Array (Shape.ZeroBased n) a -> Array (Shape.Interval n) a interval3FromChirp rng@(Shape.Interval from to) chirp = let n1 = Shape.zeroBasedSize $ Array.shape chirp in let n = n1-1 in Array.reshape rng $ Chunk.toVector $ let go rev chunkStart = if to <= chunkStart+n then if rev then Chunk.reversed $ VectorSlice.sliceVector (Slice.drop (chunkStart+n-to)) chirp else Chunk.simple $ takeShape (Shape.ZeroBased (to+1-chunkStart)) chirp else (if rev then Chunk.reversed $ VectorSlice.sliceVector (Slice.drop 1) chirp else Chunk.simple $ takeShape (Shape.ZeroBased n) chirp) <> go (not rev) (chunkStart+n) in let fromR = mod from (2*n) in if fromR < n then Chunk.simple (VectorSlice.sliceVector (Slice.sub fromR (n-fromR)) chirp) <> go True (from-fromR+n) else Chunk.reversed (VectorSlice.sliceVector (Slice.sub 1 (2*n-fromR)) chirp) <> go False (from-fromR+2*n)