mersenneforum.org (https://www.mersenneforum.org/index.php)
-   Programming (https://www.mersenneforum.org/forumdisplay.php?f=29)
-   -   Implementing the IBDWT in Python? (https://www.mersenneforum.org/showthread.php?t=27089)

 bentonsar 2021-08-19 16:23

Implementing the IBDWT in Python?

Hi fellow crunchers,

I am trying to implement the [URL="http://www.faginfamily.net/barry/Papers/Discrete%20Weighted%20Transforms.pdf"]IBDWT[/URL] in Python (for clarity not for speed) and I am having some trouble. I read through the relevant code in gwnum, gpuowl, gpulucas, and Prime Numbers: a computational perspective and I still can't seem to get it running correctly.

I am using Numpy right now for the FFTs but I am fine changing that if need be.

Has somebody already implemented one in python?

Any help would be greatly appreciated.

-Sarah

 paulunderwood 2021-08-19 17:37

See the [URL="https://mersenneforum.org/showthread.php?p=531780#post531780"]posts here[/URL]. It should be relatively easy to take the expressive Pari/GP code and convert into Python, :smile:

The code there might be for Wagstaff numbers. I am sure someone here will point out the differences for it an that for Mersenne numbers.

 bentonsar 2021-08-26 21:52

Thanks Paul.

That helped a ton. The part I am stuck on now is the digit to balanced and then back again (after the FFT) part. Is there any pseudocode anywhere for that?

-Sarah

 axn 2021-08-28 10:33

Also [URL="https://en.wikipedia.org/wiki/Signed-digit_representation"]wikipedia[/URL]

Finally, have a look at [URL="https://sourceforge.net/p/cudalucas/code/HEAD/tree/trunk/"]cudaLucas source code[/URL] for a GPU based implementation of LL test. The core FFT is based on Nvidia cuFFT library, but all the weightings, balancing, carry etc are done there.

Ask if you have any specific questions/difficulties

 bentonsar 2021-08-28 18:04

Hi Axn,

Take the example from page 498 in [URL="http://thales.doa.fmph.uniba.sk/macaj/skola/teoriapoli/primes.pdf"]Prime Numbers: A Computational Exploration[/URL] where p = 2^521 - 1, q = 521, D=16, d = (33, 33, 32, 33, 32, 33, 32, 33, 33, 32, 33, 32, 33, 32, 33, 32), a = ...

Right now I can successfully generate d and a based on those input values and my FFTs seem to be working correctly. My question is, let's say we want to have the x and y integers be 10 and 10 (or some other integer, it doesn't really matter). How do you apply those weights that were created to the integers so that they are ready for the FFT (of size D)?

I read through the x = formula that is listed in 9.5.18 but I couldn't quite get it working correctly.

-Sarah

 ewmayer 2021-09-07 00:38

Assuming your inputs are properly normalized, i.e. x,y < (2^p-1), you break each into 32/33-bit pieces according to the per-word bitlengths in your list. (Balanced-digit representation is easy to achieve here, but not necessary for learning the basics.) Then multiply each of the D words of each so-discretized input by the corresponding fractional power of 2, fwdFFT each resulting weighted input vector, dyadic-multiply the resulting length-D forward transforms, invFFT the resulting vector, and multiply the resulting D discrete convolution outputs by the D inverse weights. The resulting outputs should be close enough to pure-integers that they can be confidently rounded-to-nearest and carries propagated. Check that the carry-propagated result matches x*y (mod 2^p-1) computed using exact multiprecision integer arithmetic.

[QUOTE=bentonsar;586739]Hi Axn,

Take the example from page 498 in [URL="http://thales.doa.fmph.uniba.sk/macaj/skola/teoriapoli/primes.pdf"]Prime Numbers: A Computational Exploration[/URL] where p = 2^521 - 1, q = 521, D=16, d = (33, 33, 32, 33, 32, 33, 32, 33, 33, 32, 33, 32, 33, 32, 33, 32), a = ...

Right now I can successfully generate d and a based on those input values and my FFTs seem to be working correctly. My question is, let's say we want to have the x and y integers be 10 and 10 (or some other integer, it doesn't really matter). How do you apply those weights that were created to the integers so that they are ready for the FFT (of size D)?

I read through the x = formula that is listed in 9.5.18 but I couldn't quite get it working correctly.

-Sarah[/QUOTE]

 chalsall 2021-09-07 02:51

[QUOTE=ewmayer;587417]Assuming your inputs are properly normalized, i.e. x,y < (2^p-1).[/QUOTE]

What if that isn't a safe assumption?

 ewmayer 2021-09-07 20:50

BTW, a useful identity for computing inverse weights for mod-M(p) transform is as follows - for length-N transform, start with the weights as per the original Crandall/Fagin paper (I dislike their p = 2^q - 1 notation, since "p" in this sort of context generally implies prime, but their q is prime and primality-or-not of p remains to be established ... I instead use M(p) = 2^p-1 below):

w[j] = 2^[ceiling(j*p/N) - j*p/N], j = 0,...,N-1 .[*]

This gives w = 1, followed by a sequence of nonrepeating fractional powers of 2. E.g. for p = 263 and N = 16 we have w[j] = 2^[0,9,2,11,4,13,6,15,8,1,10,3,12,5,14,7]/16 = 2^[j*sw (mod N)] for j = 0,...,N-1, where sw denotes the number of "smallwords" in our variable-base representation; if bw = p%N is the number of "bigwords", then sw = N - bw.

If we extend the formula[*] to j = N, we get w[N] = w = 1. Instead define w[N] := 2, then observe that

w[j] * w[N-j] = 2 for j = 0,...,N .

Thus w[j] = 2/w[N-j], and the jth inverse weight can be computed as 1/w[j] = w[N-j]/2 . The 1/N multiplier needed for inverse-transform outputs can be lumped into the denominator of the RHS, thus one can define "effective inverse weights" by winv[j] = w[N-j]/(2*N); this needs just a single reciprocal 1/2N to be computed.

 bentonsar 2021-09-09 03:04

Getting Closer

[QUOTE=ewmayer;587417]Assuming your inputs are properly normalized, i.e. x,y < (2^p-1), you break each into 32/33-bit pieces according to the per-word bitlengths in your list. (Balanced-digit representation is easy to achieve here, but not necessary for learning the basics.) Then multiply each of the D words of each so-discretized input by the corresponding fractional power of 2, fwdFFT each resulting weighted input vector, dyadic-multiply the resulting length-D forward transforms, invFFT the resulting vector, and multiply the resulting D discrete convolution outputs by the D inverse weights. The resulting outputs should be close enough to pure-integers that they can be confidently rounded-to-nearest and carries propagated. Check that the carry-propagated result matches x*y (mod 2^p-1) computed using exact multiprecision integer arithmetic.[/QUOTE]

Thx for that.

My current code is below. I think I am close. I can smell victory! I know I am missing the transfer back to the balanced int (to get the actual result of the operation), but I want to get the first iteration working (4x4 mod M) before tackling that.

[code]import numpy as np

import math

q = 521
#Q = 163

D = 16

d = np.array([np.ceil((q * i) / D) - np.ceil(q*(i-1) / D) for i in range(0, D)]) #creates the bit length array
print(d)
a = np.array([2**(np.ceil(q * j/D) - (q * j/D)) for j in range(0, D)]) #creates the weight signal array
print(a)
ainv = np.array([(1.0/a[i]/D) for i in range(0, len(a))]) #creates the inverse weight signal array
print("Inverted A: " + str(ainv))

i_hiBitArr = [int(0)] * D
dsignal = [int(0)] * D

llint_signal = [int(0)] * D
llint_signal = 4

for i in range(0, q-2):
ival = llint_signal[i]
print("Ival: " + str(ival))
if (i == 0):
ival += i_hiBitArr[D-1]
else:
ival += i_hiBitArr[i]
dsignal[i] = ival * a[i]
print(str(dsignal))
C = np.fft.fftn(dsignal, norm="forward")
print("First C: " + str(C))
C = C * C
print("Squared C: " + str(C))
DD = np.fft.irfft(C, norm="backward")
print("Inversed C value: " + str(DD))
if i == 0:
sig = DD[i] * ainv[i] - 2
else:
sig = DD[i] * ainv[i]
print(sig)
llint_signal[i] = np.around(sig)
dsignal[i] = sig
#using carry steps from Crandall Book
z = llint_signal
outputarr = [None] * len(z)
carry = 0
for n in range(0, len(z)-1):
B = 2**d[n+1]
v = z[n] + carry
outputarr[n] = math.remainder(v, B)
carry = np.floor(v//B)
print("Carry " + str(n) + ": " + str(carry))
print("V " + str(n) + ": " + str(v))
print("B " + str(n) + ": " + str(B))
print("outputarr " + str(n) + ": " + str(outputarr))
print(outputarr)
exit()[/Code]

Is this implemented correctly? If so, how do I convert that outputarr to the actual result of the sxs mod M operation in the lucas lehmer test?

-Sarah

 ewmayer 2021-09-09 05:27

[QUOTE=bentonsar;587541]Is this implemented correctly? If so, how do I convert that outputarr to the actual result of the sxs mod M operation in the lucas lehmer test?

-Sarah[/QUOTE]

It'll be easier for the rest of us to tell if it's implemented correctly if you would provide a sample input, say a random number of the desired bitlength p you are working (or pair of such if you're doing a general mod-M(p) multiply rather than a squaring), the resulting N mixed-base inputs and forward-weights, and the resulting outputs after inverse-weighting but before carry propagation, and then after you've done the carries.

 All times are UTC. The time now is 08:10.