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()
