Discrete Weighted Transform - A Primitive Nth Root of Unity

[Sane]

Hello, I got pointed to this forum by OfficeShreder. I have a question I've been puzzling myself over for a while.

I am currently trying to implement the "Discrete Weighted Transform". I have reached a step where I need to determine "a primitive Nth root of unity in the appropriate domain".

I have been reading up on what roots of unity are, but I do not know what the "appropriate domain" is, nor how to determine 'g' according to this.

The research paper can be found here: http://faginfamily.net/barry/Papers/Discrete Weighted Transforms.pdf

The following formula is given:

$$X_{k} = \sum_{j=0}^{N-1} a_{j}x_{j}g^{-jk}$$

Along with that formula, in section (2.5), variable 'g', is a primitive Nth root of unity in the appropriate domain.

This is all part of the variant of Algorithm W, as referenced in section 6.
Algorithm W can be found at the end of section 3.

If someone could help me figure out how to determine variable 'g', I would be very appreciative. Thanks in advance.

 

[Sane]

I'm trying to avoid complaining, because I understand that I have no right to receive your help; it is merely my privilege. But I'm a little stressed out and frustrated as I've posted this on two separate forums and have received very little help (and in this case, none at all). Also, this is fairly time sensitive, and therefore adamently require some assistance.

If you don't know the answer that's fine, but would you be able to point me to another forum or person who may be able to help me?

Sorry if I'm sounding like a jerk in any way. >_>

 

[shmoe]

Are you familiar with the usual discrete Fourier transform? This weighted version appears to be a modification the same thing. "appropriate domain" can be whatever ring you are considering your x_i's to be a part of. You could take the complex numbers for example, a primitive nth root of unity would then be e^(2*pi*i/n). You could take the integers modulo something, they talk about number theoretic transform version in section 8 there. What domain you want to take will depend on what you are trying to do. It would probably be fine for you to just work with the complex numbers to get the general idea of what's going on.

You also might have better luck getting help in a computer science forum (sorry don't know any) where people are more likely to know about algorithms like this (I can only claim a passing familiarity). Perhaps the mersenneforum.org people would be a good place to go, they're interested in multiplying large numbers quickly, so will be better acquainted with methods like this.

 

[Sane]

You mean the usual "fast Fourier transform"? Then, no, I have not implemented the FFT before.

My x_i's are determined from the equation in section (6.5):

$$x = \sum_{j=0}^{N-1} x_{j}2^{Ceiling(qj/N)}$$

I don't know what "ring" this is part of. But what I am trying to do is what's explained in section 6, according to a variant of Algorithm W from the end of section 3 (as I've stated in the first post).

$$n = x^{2} (\bmod{p})$$

I've tried posting this on a computer science algorithm forum, but still have no replies (three forums I've posted this on now). If I keep running out of luck, I'll try the fourth forum you've pointed me to.

 

[chroot]

The fast Fourier transform is a numerical algorithm. shmoe is not talking about that.

Shmoe is talking about the discrete Fourier transform, or DFT, which is essentially the same thing as the normal continuous Fourier transform, with sums in place of integrals.

- Warren

 

[Sane]

Ahh, okay. Well I haven't implemented the DFT either. I'm only working on the DWT.

The reason that confused me is because the DWT is supposedly analogous to the FFT.

 

[shmoe]

the FFT is a way to compute the DFT with less operations than the 'naive' way of perform the transformation (just computing the terms in the sum as they appear). the DFT is a special case of the DWT, taking the weights to be 1's.

For their algorithm in section 6, use a primitive Nth complex root of unity. Your weights are irrational, so you can consider your "appropriate domain" to be the complex numbers. They give some details of a numerical example in section 6. Try to fill in the steps with exp^(i*2*Pi/N)=i, as N=4 in this example.

 

[Sane]

Where is the "numerical example in section 6" that demonstrates the step involving finding a primitive Nth root of unity?

I see an example of the entire algorithm at the end, but it skips steps 2, 3 and 4 in Algorithm W. Those steps being the ones that require a primitive Nth root of unity.

It might help slightly to show you where I'm at, with a Pythonic implementation of the algorithm.

import math

"""
A Python Implementation of the IBDWT algorithm, for transferring to C.

    Implementation - Aaron Voelker                  - 28/10/06
    Algorithm      - Richard Crandall & Barry Fagin - 1994

Algorithm W variant: Weighted convolution algorithm for x^2 modulo (2^q-1).

    1) Choose run length N >= q and establish bit-sizes b_j for digits according to (6.6).
       Represent x = {x_j}, according to (6.5), zero-padding only to N digits inclusive.
       Compute the components of the weight signal a according to (6.7).
    2) Compute X = DWT (N, a)x.
    3) Compute Z = X^2.
    4) Compute z = DWT^-1 (N, a)Z.  (This is the weighted convolution x*^a x.)
    5) Set z = Round(z), if noninteger (e.g., floating-point) FFTs were used.
    6) Adjust the digits {z_n} to the digit representation of choice.

Ex:

    Perform x^2 (mod p).

    Given,
        x = 78314567209,    q = 37,    p = 2^37-1,    N = 4

    Calculate,
        b = {10, 9, 9, 9}
        x = {553, 93, 381, 291}
        a = {1.0, 1.681792830507429, 1.4142135623730951, 1.189207115002721}
        z = {704383., 324600., 523365.0000000001, 463577.9999999999}

    Therefore,
        x^2 (mod p) = 704383 + 324600 * 2^10 + 523365 * 2^19 + 463578 * 2^28

http://faginfamily.net/barry/Papers/Discrete%20Weighted%20Transforms.pdf
"""

def ceiling (z):
    """Round z up to the nearest integer."""
    # See section (3.8).
    # The paper describes the "ceiling" function notation as int (z + 1).
    # However, it's actually supposed to be (z + n), where (n | lim -> 1).
    # Fortunately, the math.ceil function handles this accordingly.
    return int (math.ceil (z))

class algorithm_w:
    
    def __init__ (self, q, N):
        """Initiate the handle for performing a Discrete Weighted Transform as according to Algorithm W."""
        # See section (6) and end of section (3).
        # This handle will be used to execute Algorithm W, where (x^2 modulo (2^q-1)) is calculated.
        
        # q denotes that the result will be modulo (2^q - 1).
        # N is the run time of x^2 modulo (2^q - 1).
        self.q    = q
        self.intN = N        # Store as an integer for use as a counter.
        self.N    = float(N) # Store as a float for sufficient accuracy.


    '''Debugging functions.'''
    

    def sigma (self, x_j):
        """Calculate the value of x based on its base representation, x_j."""
        # See section (6.5).
        # Take the sum of each part times its base.
        # I swear, it's only for debugging purposes, really.
        return sum (x_j[j] * self.x_j_base (j) for j in range (self.intN))


    '''Worker functions.'''


    def qjN (self, j):
        """This little formula is ugly, we'll just push it off to the side."""
        return j * self.q / self.N
        
    def x_j_base (self, j):
        """Calculate the base of the j'th number in the x series."""
        # See section (6.5).
        return 2 ** ceiling (self.qjN (j))


    '''Algorithm W exclusive functions. End of section (3).'''


    def b_j (self):
        """Calculate the series b_j: The number of bits allocated for digit x_(j-1)."""
        # See section (6.6).
        # Step 1 of Algorithm W.
        # Ugh, why the **** do we need the bit sizes again?
        return [ceiling (self.qjN (j)) - ceiling (self.qjN (j-1)) for j in range (1, self.intN + 1)]

    def a_j (self):
        """Calculate the series a_j: The constant signal, required for a properly weighted convolution."""
        # See section (6.7).
        # Step 1 of Algorithm W.
        return [2 ** (ceiling (self.qjN (j)) - self.qjN (j)) for j in range (self.intN)]

    def x_j (self, x):
        """Calculate the series x_j: The base representation of x, it being the sum of its parts."""
        # Where's the section for this god damnit? Why did I have to figure it out myself?
        # Step 1 of Algorithm W.
            
        # Divide and modulate in the same manner of how we'd jump from one number system. to another
        # Do this by calculating the base for each part, effectively reversing the sigma.

        # ****, this is probably pretty slow for large values of x. Is there a faster, synthetic way of determining the base representation of x?
        # There has got to be, or else this defeats the whole purpose.

        x_j = []
        for j in range (self.intN):       # Start at the end, since we are reversing the sum.
            base = self.x_j_base (self.N - j - 1)
            x_j  = [int (x / base)] + x_j # Convert from long to integer, and add to the result list.
            x   %= base                   # Determine the newly diminished sum of the sigma notation.
        return x_j

    def dwt (self, x_j, a_j):
        """Calculate X_k = DWT(N, a_j)x_j: The Discrete Weighted Transform, X_k, of x_j according to weight signal a_j, comprised of N non-zero constants."""
        # See section (2.5).
        # Step 2 of Algorithm W.

        pass
        


def run_test ():
    x1   = 78314567209          # this number is supposedly represented as
    x_j1 = [553, 93, 381, 291]  # these four values

    q = 37                      # when these values are true
    N = 4                       # how the hell did we get this number?

    dwt_h = algorithm_w (q, N)

    bits = dwt_h.b_j ()         # calculate the bit sizes for x_j
    print bits

    signals = dwt_h.a_j ()      # calculate the constant signal a
    print signals

    x2 = dwt_h.sigma (x_j1)     # x_j1 produces the x value: 78314567209
    print x1
    print x2
    print x1 == x2              # True

    x_j2 = dwt_h.x_j (x1)       # x1 produces the x_j values: [553, 93, 381, 291]
    print x_j1
    print x_j2
    print x_j1 == x_j2          # True


if __name__ == '__main__':
    run_test()

 

[shmoe]

Where is the "numerical example in section 6" that demonstrates the step involving finding a primitive Nth root of unity?

I see an example of the entire algorithm at the end, but it skips steps 2, 3 and 4 in Algorithm W. Those steps being the ones that require a primitive Nth root of unity.

That's the example I was talking about, I didn't say they showed you what primitive root to use or how to find it. I said you can fill in their steps using the primitive root I suggested, i.e. take the data from (6.10) and you can reproduce all the numbers in their example when squaring 78314567209 (mod p) by taking exp^(i*2*Pi/4)=i to be your primitive 4th root of unity. It's been a couple of weeks but I recall it worked out (can't seem to find the scrap paper anymore).

All you were missing was the primitive root, no? Or are there other parts of the algorithm you don't follow? (I haven't read your python code)