Category Archives: math

N-ary m-sequence generator in Python

Posted on by
Reply

Maximum length sequences (or m-sequences) are bit sequences generated using maximal LFSRs (Linear Feedback Shift Registers) and are so called because they are periodic and reproduce every binary sequence that can be represented by the shift registers (i.e., for length-m registers they produce a sequence of length 2^{m}-1). Although originally defined for binary sequences, the concept of m-sequence can be extended to sequences over GF(n), n>2.

Despite being deterministically generated, m-sequences have nice statistical properties, which resemble those of white noise (they are spectrally flat, with the exception of a near-zero constant term). For this reason, m-sequences are often referred to as pseudo random noise code (PN code). Unfortunately though, unpredictability is not among those properties (see below).

The configuration of the feedback taps in a binary LFSR can be expressed in finite field arithmetic as a polynomial in GF(2^{n}). The LFSR is maximal-length if and only if the corresponding feedback polynomial is primitive. The linear complexity (LC) of a given periodic sequence is the number of cells in the shortest LFSR that can generate that sequence. The correlation between two sequences is the complex inner product of the first sequence with a shifted version of the second sequence. It is called an autocorrelation if the two sequences are the same or a cross-correlation if they are distinct. Particularly interesting for spread spectrum application are sequences (not necessarily m-sequences) that exhibit low autocorrelation and low cross-correlation (the ideal correlation function would then be the Kronecker’s delta).

Apart from spread spectrum communications, m-sequences are used in a variety of areas including digital watermarking and cryptography. As m-sequences are derived from linear recurrence relations — which lead to fairly easy cryptanalysis — they must be processed by nonlinear functions in order to be used in cryptographic applications. A good example is Trivium, a 80-bit stream cipher that consists of three shift registers of different lengths.

#!/usr/bin/python3
#
# Author: Joao H de A Franco (jhafranco@acm.org)
#
# Description: N-ary sequence generator in Python
#
# Version: 1
#
# Date: 2014-02-14
#
# License: Attribution-NonCommercial-ShareAlike 3.0 Unported
#          (CC BY-NC-SA 3.0)
#===========================================================

from functools import reduce
from fractions import Fraction
from cmath import rect,pi

def flatList(s,q):
    """Converts (list of) polynomials into (list of) elements
       of finite field mapped as integers.
         s: list of polynomials
         q: order of finite field"""
    if type(s[0]) is int:
        return s
    elif type(s[0]) is list:
        return [reduce(lambda i,j:(i*q)+j,e) for e in s]
    else:
        raise TypeError

def printListOfNumbers(S,L,*X):
    """Print list of numbers according to their particular
       format (integer, float or complex)"""
    print(S+" [",end="")
    for i,e in enumerate(L):
        if type(e) is complex:
            print("{:5.3f}{:+5.3f}i".format(e.real,e.imag),end="")
        elif type(e) is float:
            print("{:5.3f}".format(e),end="")
        elif type(e) is int:
            print("{:d}".format(e),end="")
        elif type(e) is Fraction:
            print("{:s}".format(e),end="")
        else:
            raise TypeError
        if i<len(L)-1:
            print(", ",end="")
    print("]",*X)

def xac(s,q=None,d=None,n=None):
    """Evaluates the complex auto-correlation of a periodic
       sequence with increasing delays within the period.
       Input length must correspond to a full period.
         s: q-ary periodic sequence
         q: order of finite field
         d: maximum denominator allowed in fractions
         n: number of decimal places in floating point numbers"""
    return xcc(s,s,q,d,n)
    
def xcc(s1,s2,q=None,d=None,n=None):
    """Evaluates the complex cross-correlation between two
       equally periodic q-ary sequences with increasing delays
       within the period. Input length must correspond to a full
       period.
         s: q-ary periodic sequence
         q: order of finite field
         d: maximum denominator allowed in fractions
         n: number of decimal places in floating point numbers"""
    
    def cc(s1,s2):
        """Evaluates the complex correlation between two equally
           periodic numerical sequences.        
             s1,s2: q-ary periodic sequences"""      
        assert type(s1[0]) == type(s2[0])
        if type(s1[0]) is list:
            s3 = [[(j-i)%q for i,j in zip(u,v)] for u,v in zip(s1,s2)]
            s4 = [reduce(lambda x,y:(x*q)+y,e) for e in s3]
            z = sum(rect(1,2*pi*i/q) for i in s4)/len(s1)
        elif type(s1[0]) is int:
            z = sum(rect(1,2*pi*(j-i)/q) for i,j in zip(s1,s2))/len(s1)
        else:
            raise TypeError
        zr,zi = round(z.real,n),round(z.imag,n)
        if abs(zi%1)<10**-n:
            if abs(zr-round(zr))<10**-n:
                return int(zr)
            elif Fraction(z.real).limit_denominator().denominator<=d:
                return Fraction(z.real).limit_denominator()
            else:
                return zr
        else:
            return complex(zr,zi)

    q = 2 if q is None else q
    d = 30 if d is None else d
    n = 3 if n is None else n
    assert len(s1) == len(s2)
    return [cc(s1,s2[i:]+s2[:i]) for i in range(len(s1))]

def LFSR(P,S,M,N,K):
    """Outputs K-ary sequence with N elements. Each element is
       derived from M successive values of the LFSR sequence
       generated by polynomial P and initial state S. Polynomial
       P is represented by a list of coefficients in decreasing
       power order."""
    def LFSR2():
        """Generates linear K-ary sequence according to polynomial
           P and initial state S. If P is primitive, sequence length
           is exactly one period."""
        seq,st = [S[-1]],S
        for j in range(K**len(S)-2):
            st0 = sum([i*j for i,j in zip(st,P[1:])])%K
            st = [st0]+st[:-1]
            seq += [st[-1]]
        return seq

    assert len(P) > 1 and len(P)-1 == len(S)
    s = LFSR2()
    L = len(s)
    assert M <= L
    return [s[i%L] if M == 1 else (s[i%L:]+s[:i%L])[:M] for i in range(N)] 

if __name__ == "__main__":
       
    # Binary sequence over GF(2) generated by the LFSR
    # defined by non primitive polynomial x⁴+x³+x+1.
    # This is not a m-sequence since its period is 6
    # (<15 = 2⁴-1)
    s01 = LFSR([1,1,0,1,1],[0,0,0,1],1,15,2)
    print("s01 =",s01,"\n")

s01 = [1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0]

    # Sequence over GF(2⁴) derived from sequence "s01".
    # The polynomial elements were also mapped to integers
    # in Z16, in order to show the period. 
    s04 = LFSR([1,1,0,1,1],[0,0,0,1],4,15,2)
    print("s04 =",s04,"\n")
    print("s04 =",flatList(s04,2),"\n")

s04 = [[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 1], [0, 1, 1, 1],
       [1, 1, 1, 0], [1, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1],
       [0, 0, 1, 1], [0, 1, 1, 1], [1, 1, 1, 0], [1, 1, 0, 0],
       [1, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]]

s04 = [8, 1, 3, 7, 14, 12, 8, 1, 3, 7, 14, 12, 9, 2, 4] 

    # m-sequence over GF(2) generated by the LFSR defined
    # by the primitive polynomial x⁴+x+1. Its period is
    # equal to 2⁴-1=15
    s11 = LFSR([1,1,0,0,1],[0,0,0,1],1,15,2)
    print("s11 =",s11,"\n")

s11 = [1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0] 

    # Autocorrelation of m-sequence "s11". This function
    # is two-valued and has the same period as sequence
    # "s11"
    printListOfNumbers("ac(s11) =",xac(s11),"\n")

ac(s11) = [1, -1/15, -1/15, -1/15, -1/15, -1/15, -1/15, -1/15,
           -1/15, -1/15, -1/15, -1/15, -1/15, -1/15, -1/15]

    # Sequence over GF(2²) derived from m-sequence "s11"
    s12 = LFSR([1,0,0,1,1],[0,0,0,1],2,15,2)
    print("s12 =",s12,"\n")

s12 = [[1, 0], [0, 0], [0, 0], [0, 1], [1, 0], [0, 0], [0, 1],
       [1, 1], [1, 0], [0, 1], [1, 0], [0, 1], [1, 1], [1, 1], [1, 1]] 

    
    # Autocorrelation of sequence "s12"
    printListOfNumbers("ac(s12) =",xac(s12,4),"\n")

ac(s12) = [1, 7/15, 7/15, 7/15, 7/15, 7/15, 7/15, 7/15, 7/15,
           7/15, 7/15, 7/15, 7/15, 7/15, 7/15]

    # Sequence over GF(2³) derived from m-sequence "s11"
    s13 = LFSR([1,0,0,1,1],[0,0,0,1],3,15,2)
    print("s13 =",s13,"\n")

s13 = [[1, 0, 0], [0, 0, 0], [0, 0, 1], [0, 1, 0], [1, 0, 0],
       [0, 0, 1], [0, 1, 1], [1, 1, 0], [1, 0, 1], [0, 1, 0],
       [1, 0, 1], [0, 1, 1], [1, 1, 1], [1, 1, 1], [1, 1, 0]] 

    # Autocorrelation of sequence "s13"
    printListOfNumbers("ac(s13) =",xac(s13,8,d=1000),"\n")

ac(s13) = [1, 0.844, 0.844, 0.844, 0.844, 0.844, 0.844, 0.844,
           0.844, 0.844, 0.844, 0.844, 0.844, 0.844, 0.844] 

    # Sequence over GF(2⁴) derived from the m-sequence "s11"
    s14 = LFSR([1,0,0,1,1],[0,0,0,1],4,15,2)
    print("s14 =",s14,"\n")

s14 = [[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0],
       [1, 0, 0, 1], [0, 0, 1, 1], [0, 1, 1, 0], [1, 1, 0, 1],
       [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 1], [0, 1, 1, 1],
       [1, 1, 1, 1], [1, 1, 1, 0], [1, 1, 0, 0]] 

    # Autocorrelation of sequence "s14" 
    printListOfNumbers("ac(s14) =",xac(s14,16,d=1000),"\n")

ac(s14) = [1, 0.959, 0.959, 0.959, 0.959, 0.959, 0.959, 0.959,
           0.959, 0.959, 0.959, 0.959, 0.959, 0.959, 0.959] 

    # Binary sequence over GF(2) generated by the LFSR
    # defined by the primitive polynomial (x⁴+x³+1).
    # This is also a m-sequence with period equal to
    # 15 = 2⁴-1
    s21 = LFSR([1,0,0,1,1],[0,0,0,1],1,15,2)
    print("s21 =",s21,"\n")

s21 = [1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1]

    # Autocorrelation of sequence "s21"       
    printListOfNumbers("ac(s21) =",xac(s21),"\n")

ac(s21) = [1, -1/15, -1/15, -1/15, -1/15, -1/15, -1/15,
           -1/15, -1/15, -1/15, -1/15, -1/15, -1/15, -1/15, -1/15]

    # Cross-correlation between m-sequences "s11" and "s21".
    # This function is four-valued and has the same period
    # as sequences "s11" and "s21"
    printListOfNumbers("cc(s11,s21) =",xcc(s11,s21),"\n")

cc(s11,s21) = [-1/15, 1/5, -1/15, 7/15, 1/5, -1/3, -1/15,
               1/5, 7/15, -1/3, 1/5, -1/3, -1/3, -1/15, -1/15] 

    # m-sequence over GF(3) generated by the LFSR defined
    # by the primitive polynomial x²+2x+1. Its period is
    # equal to 3²-1=8 
    s31 = LFSR([1,2,1],[1,0],1,8,3)
    print("s31 =",s31,"\n")

s31 = [0, 1, 2, 2, 0, 2, 1, 1]

    # Autocorrelation of sequence "s31"
    printListOfNumbers("ac(s31) =",xac(s31,3),"\n")

ac(s31) = [1, -1/8, -1/8, -1/8, -1/8, -1/8, -1/8, -1/8]

    # Sequence over GF(3²) derived from m-sequence "s31"
    s32 = LFSR([1,2,1],[1,0],2,8,3)
    print("s32 =",s32,"\n")

s32 = [[0, 1], [1, 2], [2, 2], [2, 0], [0, 2], [2, 1], [1, 1], [1, 0]]

    # Autocorrelation of sequence "s32"
    printListOfNumbers("ac(s32) =",xac(s32,3),"\n")

ac(s32) = [1, -1/8, -1/8, -1/8, -1/8, -1/8, -1/8, -1/8]

    # Sequence over GF(3³) derived from m-sequence "s31"
    s33 = LFSR([1,2,1],[1,0],3,8,3)
    print("s33 =",s33,"\n")

s33 = [[0, 1, 2], [1, 2, 2], [2, 2, 0], [2, 0, 2],
       [0, 2, 1], [2, 1, 1], [1, 1, 0], [1, 0, 1]]

    # Autocorrelation of sequence "s33"
    printListOfNumbers("ac(s33) =",xac(s33,3),"\n")

ac(s33) = [1, -1/8, -1/8, -1/8, -1/8, -1/8, -1/8, -1/8]

    # m-sequence over GF(3) generated by the LFSR defined
    # by the primitive polynomial x³+x²+2. Its period is
    # equal to 3³-1=26 
    s41 = LFSR([1,1,0,2],[1,0,0],1,26,3)
    print("s41 =",s41,"\n")

s41 = [0, 0, 1, 1, 1, 0, 2, 1, 1, 2, 1, 0, 1, 0, 0, 2, 2, 2, 0, 1,
       2, 2, 1, 2, 0, 2] 

    # Autocorrelation of m-sequence "s41". This function is
    # two-valued and has the same period as sequence "s41"    
    printListOfNumbers("ac(s41) =",xac(s41,3),"\n")

ac(s41) = [1, -1/26, -1/26, -1/26, -1/26, -1/26, -1/26, -1/26, -1/26, 
           -1/26, -1/26, -1/26, -1/26, -1/26, -1/26, -1/26, -1/26, -1/26,
           -1/26, -1/26, -1/26, -1/26, -1/26, -1/26, -1/26, -1/26]

    # Sequence over GF(3²) derived from m-sequence "s41"
    s42 = LFSR([1,1,0,2],[1,0,0],2,26,3)
    print("s42 =",s42,"\n")

s42 = [[0, 0], [0, 1], [1, 1], [1, 1], [1, 0], [0, 2],
       [2, 1], [1, 1], [1, 2], [2, 1], [1, 0], [0, 1],
       [1, 0], [0, 0], [0, 2], [2, 2], [2, 2], [2, 0],
       [0, 1], [1, 2], [2, 2], [2, 1], [1, 2], [2, 0],
       [0, 2], [2, 0]]

    # Autocorrelation of sequence "s42"
    printListOfNumbers("ac(s42) =",xac(s42,9),"\n")

ac(s42) = [1, 0.701, 0.701, 0.701, 0.701, 0.701, 0.701,
           0.701, 0.701, 0.701, 0.701, 0.701, 0.701, 0.838,
           0.701, 0.701, 0.701, 0.701, 0.701, 0.701, 0.701,
           0.701, 0.701, 0.701, 0.701, 0.701]

    # Sequence over GF(3³) derived from m-sequence "s31"
    s43 = LFSR([1,1,0,2],[1,0,0],3,26,3)
    print("s43 =",s43,"\n")

s43 = [[0, 0, 1], [0, 1, 1], [1, 1, 1], [1, 1, 0],
       [1, 0, 2], [0, 2, 1], [2, 1, 1], [1, 1, 2],
       [1, 2, 1], [2, 1, 0], [1, 0, 1], [0, 1, 0],
       [1, 0, 0], [0, 0, 2], [0, 2, 2], [2, 2, 2],
       [2, 2, 0], [2, 0, 1], [0, 1, 2], [1, 2, 2],
       [2, 2, 1], [2, 1, 2], [1, 2, 0], [2, 0, 2],
       [0, 2, 0], [2, 0, 0]]

    # Autocorrelation of sequence "s43"
    printListOfNumbers("ac(s43) =",xac(s43,27),"\n")

ac(s43) = [1, 0.963, 0.963, 0.963, 0.963, 0.963, 0.963,
           0.963, 0.963, 0.963, 0.963, 0.963, 0.963, 0.981,
           0.963, 0.963, 0.963, 0.963, 0.963, 0.963, 0.963,
           0.963, 0.963, 0.963, 0.963, 0.963]

    # Sequence over Z20 with ideal autocorrelation function
    # (Kronecker's delta)
    s51 = [t**2%20 for t in range(10)]
    print("s51 =",s51,"\n")

s51 = [0, 1, 4, 9, 16, 5, 16, 9, 4, 1]

    # Autocorrelation of sequence "s51"  
    printListOfNumbers("xac(s51) =",xac(s51,20),"\n")

xac(s51) = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0] 

    # Sequence over Z7 with ideal autocorrelation function
    # (Kronecker's delta)
    s61 = [(t*(t+1)//2)%7 for t in range(7)]
    print("s61 =",s61,"\n")

s61 = [0, 1, 3, 6, 3, 1, 0]

    # Autocorrelation of sequence "s61"  
    printListOfNumbers("xac(s61) =",xac(s61,7),"\n")

xac(s61) = [1, 0, 0, 0, 0, 0, 0]

Fast Fourier Transform in Haskell

Posted on by
Reply

The Discrete Fourier Transform (DFT) converts a finite list of equally-spaced samples of a function into the list of coefficients of a finite combination of complex sinusoids, ordered by their frequencies, that has those same sample values.

The DFT is used to perform Fourier analysis in different applications like digital signal processing, polynomial multiplication, partial differential equations solving, large integers multiplication etc. It has also been used in statistical tests of random and pseudorandom number generators for cryptographic applications (see NIST Special Publication 800-22).

The Haskell function dft (shown below) evaluates the DFT for any number of points. Since the vectors [e^{2{\pi}kn/N} | n = 0,1..N-1] constitute an orthogonal basis over the set of N-dimensional complex vectors, the inverse DFT can be calculated by the same algorithm using the complex conjugate vectors (see function idft). In order to make the forward and the inverse transforms unitary, a normalization factor equal to \sqrt{1/N} is used in both functions. 

-------------------------------------------------------------
--
-- Author: Joao H de A Franco (jhafranco@acm.org)
--
-- Description: DFT implementation in Haskell
--
-- Date: 2013-Mar-28
--
-- License: Attribution-NonCommercial-ShareAlike 3.0 Unported
--          (CC BY-NC-SA 3.0)
--
-------------------------------------------------------------
module Aux.DFT (dft,idft) where
import Data.Complex (Complex((:+)))

dft, idft :: [Complex Float] -> [Complex Float]
dft = dft' 1
idft = dft' (-1)
           
dft' :: Float -> [Complex Float] -> [Complex Float]
dft' e x = let t = length x
               y = 1 / fromIntegral t
               dft'' k = let w = fromIntegral k * 2 * e * pi * y
                         in (sqrt y :+ 0) *
                             sum (zipWith (\n z -> z * exp (0 :+ n * w)) [0..] x)                         
           in map dft'' [0..t-1]

The DFT can be seen as a linear transformation where an real-valued N-tuple is multiplied by an N x N matrix returning a complex-valued N-tuple. This matrix operation takes N² scalar multiplications, which makes it very expensive to deal with large data sets.

Fast Fourier Transform is the generic name given to algorithms that compute the DFT (or its inverse) faster than O(N²). The most commonly used FFT algorithm is due to Cooley & Tukey (1965), a “divide and conquer” procedure that recursively breaks down a DFT of any composite size N = N1 * N2 into many smaller DFTs of sizes N1 and N2. The Cooley-Tukey and other FFTs schemes reduce the algorithm complexity from O(N²) to O(N log N).

The advent of FFT algorithms has revolutionized the use of the Discrete Fourier Transform, making it possible to use it on a variety of areas: applied mechanics, biomedical engineering, communications, electromagnetics, finance, instrumentation, numerical methods, radar, signal processing, sonics and acoustics etc.

It’s interesting to note, however, that the general idea behind the FFT algorithms was invented by Carl Friedrich Gauss (“The Prince of Mathematicians”) around 1805, when he was calculating the orbit of an asteroid. Amazing!

The Cooley–Tukey algorithm is usually employed to recursively divide a N-point transform into two N/2-point transforms and for this reason it is limited to cases where N is a power of two. The Haskell code below implements this technique, called decimation-in-time radix-2 (fft implements the forward transform and ifft the inverse one).

-------------------------------------------------------------
--
-- Author: Joao H de A Franco (jhafranco@acm.org)
--
-- Description: FFT implementation in Haskell
--
-- Date: 2013-Mar-28
--
-- License: Attribution-NonCommercial-ShareAlike 3.0 Unported
--          (CC BY-NC-SA 3.0)
--
-------------------------------------------------------------
module Aux.FFT (fft,ifft) where
import Data.Complex(Complex((:+)))

fft, ifft :: [Complex Float] -> [Complex Float]
fft = fft' 1
ifft = fft' (-1)

fft' :: Float -> [Complex Float] -> [Complex Float]
fft' e xs = if pow2 (length xs) 
            then let z = sqrt (1 / fromIntegral (length xs))
                  in map ((z :+ 0) *) $ fft'' e xs
            else error "Number of points is not a power of 2"
            where pow2 n
                     | n == 1 || n == 2 = True
                     | otherwise = n `mod` 2 == 0 && pow2 (n `div` 2)

fft'' :: Float -> [Complex Float] -> [Complex Float]
fft'' _ [] = []
fft'' _ [x] = [x]
fft'' e xs = fft'' e (evens xs) <+> t (fft'' e (odds xs))
  where (<+>) r s = zipWith (+) (r ++ r) (s ++ map negate s)
        evens [] = []
        evens [u] = [u]
        evens (v:_:vs) = v:evens vs
        odds = evens . drop 1
        n = 2 * pi / fromIntegral (length xs)
        t = zipWith (\k z -> z * exp (0 :+ k * n * e)) ([0..] :: [Float])

Besides classic references on this subject such as Brigham’s and Bracewell’s books, I recommend the excellent lectures on The Fourier Transforms and Its Applications (EE261) by Prof. Brad Osgood from Stanford University. As a matter of fact, I’ve watched only lectures 20 to 22, as they cover DFT and the basics of FFT algorithms, but most probably the remaining lectures should be equally good. Extensive course materials (syllabus, handouts, past exams, exercises & solutions) are available.

The box below shows some simple tests to check the correctness of the above implementations.

-------------------------------------------------------------
--
-- Author: Joao H de A Franco (jhafranco@acm.org)
--
-- Description: Test module for FFT implementation in Haskell
--
-- Date: 2013-Mar-28
--
-- License: Attribution-NonCommercial-ShareAlike 3.0 Unported
--          (CC BY-NC-SA 3.0)
--
-------------------------------------------------------------
module Main where
import Data.Complex (Complex((:+)),realPart,imagPart)
import Aux.DFT (dft,idft)
import Aux.FFT (fft,ifft)

main :: IO ()
main = print (if test1 && test2 && test3 then "ok" else "nok")

-- 7-point DFT (prime number) compared to Mathematica's evaluation
test1 :: Bool
test1 = let mathematica1 = [3.40168 :+ 0, (-0.566947) :+ 0.564962,
                            (-0.566947) :+ (-1.41899), (-0.566947) :+ 0.645967,
                            (-0.566947) :+ (-0.645967), (-0.566947) :+ 1.41899,
                            (-0.566947) :+ (-0.564962)]
        in map (roundComplex 3 0.001) (dft [0,1,2,3,0,1,2]) == 
           map (roundComplex 3 0.001) mathematica1

-- 8-point FFT (power of 2) compared to Mathematica's evaluation
test2 :: Bool
test2 = let mathematica2 = [4.24264 :+ 0, 0 :+ 0,
                            (-1.41421) :+ (-1.41421), 0 :+ 0, 
                            (-1.41421) :+ 0, 0 :+ 0,
                            (-1.41421) :+ 1.41421, 0 :+ 0]
        in map (roundComplex 3 0.001) (fft [0,1,2,3,0,1,2,3]) == 
           map (roundComplex 3 0.001) mathematica2

-- 512-point DFT compared to FFT
test3 :: Bool
test3 = let list = take 512 $ cycle ([0,1,2,3] :: [Complex Float])
            z1 = map (roundComplex 2 0.001) $ dft list
            z2 = map (roundComplex 2 0.001) $ fft list
         in z1 == z2

-- 128-point time-shifted impulse forward & inverse DFT (auto-test)
test4 :: Bool
test4 = all (==True) t4
        where t4 = do
              n <- [0..127]
              let x = impulse n 128 :: [Complex Float]
              return $ map (roundComplex 3 0.001) (idft (dft x)) == x

-- 1024-point time-shifted impulse forward & inverse FFT (auto-test)         
test5 :: Bool
test5 = all (==True) t5
        where t4 = do
              n <- [0..1023]
              let x = impulse n 1024 :: [Complex Float]
              return $ map (roundComplex 3 0.001) (ifft (fft x)) == x

-- Discrete Dirac delta generator
impulse :: Num a => Int -> Int -> [a]
impulse n m = replicate n 0 ++ 1:replicate (m - n - 1) 0

-- Rounding function
roundComplex :: Int -> Float -> Complex Float -> Complex Float
roundComplex n e x = let zeroFloat e' f = if abs f < e' then 0 else f
                         trim = roundFloat n . zeroFloat e
                         roundFloat m y = fromIntegral (round $ y * 10 ^ m :: Int) / 10 ^^ m
                     in trim (realPart x) :+ trim (imagPart x)

You can run the above FFT script to evaluate the forward and inverse FFT for a small data set using the facilities provided by Codepad.org, an online compiler/interpreter for Haskell and other languages.

A quick benchmark shows the speed-up achieved by the radix-2 FFT over the classic DFT: 18 ms vs. 2.5 s for 4,096 points (two orders of magnitude!). FFTW, a C subroutine library for computing the DFT in one or more dimensions, is much faster, as it takes only 0.7 ms for processing 4,096 points on the same Linux system.

Multiplication over the binary finite field GF(2^m)

Posted on by
Reply

The multGF2() function shown in the Python script below implements the element (polynomial) multiplication over a binary finite field GF(2^{m}). The second function, setGF2(), sets the three constants needed for its colleague to perform its multiplication task: "mask1" and "mask2" (used in “and” operations) and "polyred", a polynomial constant used in the polynomial reduction of the product. These constants are defined from the parameters passed to setGF2(): "degree" (the extension degree of the binary field GF(2^{m})) and "irPoly" (the irreducible polynomial used for the product reduction). The version of multGF2() in the post AES implementation in Python, dubbed “mult()“, was configured specifically for the AES binary finite field used in the MixColumns operation, GF(2^{8})/x^{8} + x^{4} + x^{3} + x + 1.

#!/usr/bin/python3
#
# Author: Joao H de A Franco (jhafranco@acm.org)
#
# Description: Binary finite field multiplication in Python 3
#
# Date: 2012-02-16
#
# License: Attribution-NonCommercial-ShareAlike 3.0 Unported
#          (CC BY-NC-SA 3.0)
#===========================================================
from functools import reduce

# constants used in the multGF2 function
mask1 = mask2 = polyred = None

def setGF2(degree, irPoly):
    """Define parameters of binary finite field GF(2^m)/g(x)
       - degree: extension degree of binary field
       - irPoly: coefficients of irreducible polynomial g(x)
    """
    global mask1, mask2, polyred
    mask1 = mask2 = 1 << degree
    mask2 -= 1
    if sum(irPoly) <= len(irPoly):
        polyred = reduce(lambda x, y: (x << 1) + y, irPoly[1:])    
    else:
        polyred = poly2Int(irPoly[1:]) 
        
def multGF2(p1, p2):
    """Multiply two polynomials in GF(2^m)/g(x)"""
    p = 0
    while p2:
        if p2 & 1:
            p ^= p1
        p1 <<= 1
        if p1 & mask1:
            p1 ^= polyred
        p2 >>= 1
    return p & mask2

#=============================================================================
#                        Auxiliary formatting functions
#=============================================================================
def int2Poly(bInt):
    """Convert a "big" integer into a "high-degree" polynomial"""
    exp = 0
    poly = []
    while bInt:
        if bInt & 1:
            poly.append(exp)
        exp += 1
        bInt >>= 1
    return poly[::-1]

def poly2Int(hdPoly):
    """Convert a "high-degree" polynomial into a "big" integer"""
    bigInt = 0
    for exp in hdPoly:
        bigInt += 1 << exp
    return bigInt

def i2P(sInt):
    """Convert a "small" integer into a "low-degree" polynomial"""
    return [(sInt >> i) & 1 for i in reversed(range(sInt.bit_length()))]

def p2I(ldPoly):
    """Convert a "low-degree" polynomial into a "small" integer"""
    return reduce(lambda x, y: (x << 1) + y, ldPoly)

def ldMultGF2(p1, p2):
    """Multiply two "low-degree" polynomials in GF(2^n)/g(x)"""
    return multGF2(p2I(p1), p2I(p2))

def hdMultGF2(p1, p2):
    """Multiply two "high-degree" polynomials in GF(2^n)/g(x)"""
    return multGF2(poly2Int(p1), poly2Int(p2))

if __name__ == "__main__":
  
    # Define binary field GF(2^3)/x^3 + x + 1
    setGF2(3, [1, 0, 1, 1])
    
    # Alternative way to define GF(2^3)/x^3 + x + 1
    setGF2(3, i2P(0b1011))    
    
    # Check if (x + 1)(x^2 + 1) == x^2
    assert ldMultGF2([1, 1], [1, 0, 1]) == p2I([1, 0, 0])

    # Check if (x^2 + x + 1)(x^2 + 1) == x^2 + x
    assert ldMultGF2([1, 1, 1], [1, 0, 1]) == p2I([1, 1, 0])
    
    # Define binary field GF(2^8)/x^8 + x^4 + x^3 + x + 1
    setGF2(8, [1, 0, 0, 0, 1, 1, 0, 1, 1])
    
    # Alternative way to define GF(2^8)/x^8 + x^4 + x^3 + x + 1
    setGF2(8, i2P(0b100011011))
    
    # Check if (x)(x^7 + x^2 + x + 1) == x^4 + x^2 + 1
    assert ldMultGF2([1, 0], [1, 0, 0, 0, 0, 1, 1, 1]) == p2I([1, 0, 1, 0, 1])
    
    # Check if (x + 1)(x^6 + x^5 + x^3 + x^2 + x) == x^7 + x^5 + x^4 + x
    assert ldMultGF2([1, 1], [1, 1, 0, 1, 1, 1, 0]) == p2I([1, 0, 1, 1, 0, 0, 1, 0])

    # Define binary field GF(2^571)/x^571 + x^10 + x^5 + x^2 + x
    setGF2(571, [571, 10, 5, 2, 1])

    # Calculate the product of two polynomials in GF(2^571)/x^571 + x^10 + x^5 + x^2 + x,
    # x^518 + x^447 + x^320 + x^209 + x^119 + x + 1 and x^287 + x^145 + x^82 + + x^44
    print(int2Poly(hdMultGF2([518, 447, 320, 209, 119, 1, 0], [287, 145, 82, 44])))

The output of the above script is:

[562, 529, 496, 491, 465, 406, 402, 364, 354, 291, 288, 287, 264, 253, 244,
239, 235, 201, 173, 168, 165, 164, 163, 146, 145, 102, 97, 94, 93, 83, 82,
46, 45, 44, 41, 39, 38, 37, 34, 30, 26, 23, 22]

This means that the product of the polynomials x^{518} + x^{447} + x^{320} + x^{209} + x^{119} + x + 1 and x^{287} + x^{145} + x^{82} + x^{44} is:

x^562 + x^529 + x^496 + x^491 + x^465 + x^406 + x^402 + x^364 + x^354 +
x^291 + x^288 + x^287 + x^264 + x^253 + x^244 + x^239 + x^236 + x^235 +
x^201 + x^173 + x^168 + x^165 + x^164 + x^163 + x^146 + x^145 + x^102 +
x^97 + x^94 + x^93 + x^83 + x^82 + x^46 + x^45 + x^44 + x^41 + x^39 +
x^38 + x^37 + x^34 + x^30 + x^26 + x^23 + x^22

Generating Permutations From Factoradic Numbers

Posted on by
Reply

This is the title of an interesting tutorial by Prof. Joel Guilherme S. Filho. It describes how to generate permutations of a set of n elements from a given decimal index based on its representation as a “factoradic number”, i.e. a number represented in the mixed base [(n-1)!, (n-2)!, … , 2!, 1!, 0!].

Thanks to the list manipulation facilities and the support for arbitrary long integers provided by the Python programming language, it’s easy to write a short script (shown below) that generates permutations even for large indexes. Please read his tutorial for additional details.

#!/usr/bin/python3
#
# Author: Joao H de A Franco (jhafranco@acm.org)
#
# Description: Generating Permutations From Factoradic Numbers
#
# Date: 2012-02-06
#
# License: Attribution-NonCommercial-ShareAlike 3.0 Unported
#          (CC BY-NC-SA 3.0)
#=============================================================
from sys import argv
from math import factorial

N = int(argv[1])

n = 0
while N > factorial(n):
    n += 1

S = [i for i in range(n)]
Factoradic, Permutation = [None] * n, [None] * n

v, s = N, S[:]
for i in range(n):
    r, v = divmod(v, factorial(n - i - 1))
    Factoradic[i], Permutation[i] = r, s.pop(r)

print("For N = {:d}\n".format(N))
print("n = {:d}\n".format(n))
print("Permutation[0] = {:s}\n".format(S))
print("Factoradic[N] = {:s}\n".format(Factoradic))
print("Permutation[N] = {:s}\n".format(Permutation))

Executing this script with 51 as the command-line argument, we get the output below:

For N = 51

n = 5

Perm[0] = [0, 1, 2, 3, 4]

Factoradic[N] = [2, 0, 1, 1, 0]

Permutation[N] = [2, 0, 3, 4, 1]

Trying now a big number like the Mersenne number \small 2^{607} - 1 (which happens to be the Mersenne prime # 14), we get the following results:

For N = 53113799281676709868958820655246862732959311772703192319944413 ... 29486246501015346579337652707239409519978766587351943831270835393219031728127

n = 113

Permutation[0] = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112]

Factoradic[N] = [2, 77, 30, 88, 71, 102, 97, 70, 4, 12, 75, 71, 27, 9, 91, 43, 75, 58, 30, 73, 10, 11, 68, 33, 77, 79, 50, 80, 41, 77, 48, 5, 32, 78, 25, 74, 0, 9, 49, 3, 43, 14, 6, 20, 2, 6, 61, 14, 29, 7, 37, 41, 5, 15, 30, 54, 26, 0, 41, 19, 29, 50, 6, 6, 2, 25, 6, 8, 44, 10, 26, 13, 16, 0, 14, 4, 5, 24, 24, 30, 22, 6, 6, 29, 10, 24, 12, 23, 3, 20, 15, 20, 16, 11, 12, 7, 2, 15, 10, 5, 5, 8, 6, 7, 2, 0, 0, 1, 0, 1, 0, 1, 0]

Permutation[N] = [2, 78, 31, 91, 73, 107, 102, 72, 5, 14, 82, 77, 30, 11, 104, 49, 87, 65, 36, 88, 13, 16, 84, 42, 98, 101, 61, 106, 52, 103, 60, 7, 43, 111, 34, 108, 0, 17, 69, 6, 63, 24, 12, 35, 4, 18, 99, 28, 53, 20, 67, 75, 15, 37, 58, 105, 54, 1, 85, 45, 64, 110, 22, 23, 9, 62, 26, 32, 112, 39, 74, 46, 51, 3, 50, 25, 29, 83, 86, 96, 80, 38, 40, 109, 55, 94, 59, 95, 21, 92, 76, 97, 81, 66, 70, 47, 19, 100, 71, 44, 48, 79, 57, 89, 27, 8, 10, 41, 33, 68, 56, 93, 90]

Math poetry: prime numbers and Nemo

Posted on by
Reply

“How many computing operations have been performed by all machines across all of world history?” … To put it roughly, right around the turn of the century (2000 AD), about one mole — that is, the Avogadro number 6 \times 10^{23} of chemistry, call it 10^{24} — is the total operation count for all machines for all of history. In spite of the usual mystery and awe that surrounds the notion of industrial and government super-computing, it is the huge collection of personal computers that allows this 10^{24}, this mole. The relevance is that a task such as trial dividing an integer N \simeq 10^{50} directly for prime factors is hopeless in the sense that one would essentially have to replicate the machine effort of all time. To convey an idea of scale, a typical instance of the deepest factoring or primality-proving runs of the modern era involves perhaps 10^{16} to 10^{18} machine operations. Similarly, a full-length graphically rendered synthetic movie of today—for example, the 2003 Pixar/Disney movie Finding Nemo—involves operation counts in the 10^{18} range. It is amusing that for this kind of Herculean machine effort one may either obtain a single answer (a factor, maybe even a single “prime/composite” decision bit) or create a full-length animated feature whose character is as culturally separate from a one-bit answer as can be. It is interesting that a computational task of say 10^{18} operations is about one ten-millionth of the overall historical computing effort by all Earth-bound machinery.

[excerpt from Crandall & Pomerance, “Prime Numbers – A Computational Perspective”, 2nd edition, Springer (2005), p. 5]

Math curiosity: expansion of 1/998001

Posted on by
Reply 
1/998001 = 0.0000010020030040050060070080090100110120130140150160170180190200210
22023024025026027028029030031032033034035036037038039040041042043044045046047048
04905005105205305405505605705805906006106206306406506606706806907007107207307407
50760770780790800810820830840850860870880890900910920930940950960970980991001011
02103104105106107108109110111112113114115116117118119120121122123124125126127128
12913013113213313413513613713813914014114214314414514614714814915015115215315415
51561571581591601611621631641651661671681691701711721731741751761771781791801811
82183184185186187188189190191192193194195196197198199200201202203204205206207208
20921021121221321421521621721821922022122222322422522622722822923023123223323423
52362372382392402412422432442452462472482492502512522532542552562572582592602612
62263264265266267268269270271272273274275276277278279280281282283284285286287288
28929029129229329429529629729829930030130230330430530630730830931031131231331431
53163173183193203213223233243253263273283293303313323333343353363373383393403413
42343344345346347348349350351352353354355356357358359360361362363364365366367368
36937037137237337437537637737837938038138238338438538638738838939039139239339439
53963973983994004014024034044054064074084094104114124134144154164174184194204214
22423424425426427428429430431432433434435436437438439440441442443444445446447448
44945045145245345445545645745845946046146246346446546646746846947047147247347447
54764774784794804814824834844854864874884894904914924934944954964974984995005015
02503504505506507508509510511512513514515516517518519520521522523524525526527528
52953053153253353453553653753853954054154254354454554654754854955055155255355455
55565575585595605615625635645655665675685695705715725735745755765775785795805815
82583584585586587588589590591592593594595596597598599600601602603604605606607608
60961061161261361461561661761861962062162262362462562662762862963063163263363463
56366376386396406416426436446456466476486496506516526536546556566576586596606616
62663664665666667668669670671672673674675676677678679680681682683684685686687688
68969069169269369469569669769869970070170270370470570670770870971071171271371471
57167177187197207217227237247257267277287297307317327337347357367377387397407417
42743744745746747748749750751752753754755756757758759760761762763764765766767768
76977077177277377477577677777877978078178278378478578678778878979079179279379479
57967977987998008018028038048058068078088098108118128138148158168178188198208218
22823824825826827828829830831832833834835836837838839840841842843844845846847848
84985085185285385485585685785885986086186286386486586686786886987087187287387487
58768778788798808818828838848858868878888898908918928938948958968978988999009019
02903904905906907908909910911912913914915916917918919920921922923924925926927928
92993093193293393493593693793893994094194294394494594694794894995095195295395495
595695795895996096196296396496596696796896997097197297397497

Note that the result contains the sequence “000”, “111”, “222” up to “999” (with exception of “998”), after which it starts again. This expansion was obtained evaluating the expression 1 + (1/999^2).n(digits = 3007) with the help of Sage Online (sagenb.org). A “1” was added to the expression to show the first zeroes, otherwise Sage uses exponential notation and omits them.