mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Programming

Reply
 
Thread Tools
Old 2021-08-19, 16:23   #1
bentonsar
 
"Sarah B"
Aug 2021
Washington, Seattle

22 Posts
Question Implementing the IBDWT in Python?

Hi fellow crunchers,

I am trying to implement the IBDWT 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
bentonsar is offline   Reply With Quote
Old 2021-08-19, 17:37   #2
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

3,863 Posts
Default

See the posts here. It should be relatively easy to take the expressive Pari/GP code and convert into Python,

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.

Last fiddled with by paulunderwood on 2021-08-19 at 18:59
paulunderwood is offline   Reply With Quote
Old 2021-08-26, 21:52   #3
bentonsar
 
"Sarah B"
Aug 2021
Washington, Seattle

22 Posts
Default

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
bentonsar is offline   Reply With Quote
Old 2021-08-28, 10:33   #4
axn
 
axn's Avatar
 
Jun 2003

514910 Posts
Default

This thread contains some description for balanced digit representation

Also wikipedia

Finally, have a look at cudaLucas source code 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

Last fiddled with by axn on 2021-08-28 at 10:33
axn is online now   Reply With Quote
Old 2021-08-28, 18:04   #5
bentonsar
 
"Sarah B"
Aug 2021
Washington, Seattle

410 Posts
Default

Hi Axn,

Thx for the links.

Take the example from page 498 in Prime Numbers: A Computational Exploration 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
bentonsar is offline   Reply With Quote
Old 2021-09-07, 00:38   #6
ewmayer
2ω=0
 
ewmayer's Avatar
 
Sep 2002
República de California

101101100010102 Posts
Default

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:
Originally Posted by bentonsar View Post
Hi Axn,

Thx for the links.

Take the example from page 498 in Prime Numbers: A Computational Exploration 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 is offline   Reply With Quote
Old 2021-09-07, 02:51   #7
chalsall
If I May
 
chalsall's Avatar
 
"Chris Halsall"
Sep 2002
Barbados

37×269 Posts
Default

Quote:
Originally Posted by ewmayer View Post
Assuming your inputs are properly normalized, i.e. x,y < (2^p-1).
What if that isn't a safe assumption?
chalsall is offline   Reply With Quote
Old 2021-09-07, 20:50   #8
ewmayer
2ω=0
 
ewmayer's Avatar
 
Sep 2002
República de California

1165810 Posts
Default

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[0] = 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[0] = 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.

Last fiddled with by ewmayer on 2021-09-07 at 21:00
ewmayer is offline   Reply With Quote
Old 2021-09-09, 03:04   #9
bentonsar
 
"Sarah B"
Aug 2021
Washington, Seattle

1002 Posts
Default Getting Closer

Quote:
Originally Posted by ewmayer View Post
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.
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[0] = 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()


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

Last fiddled with by paulunderwood on 2021-09-09 at 06:07
bentonsar is offline   Reply With Quote
Old 2021-09-09, 05:27   #10
ewmayer
2ω=0
 
ewmayer's Avatar
 
Sep 2002
República de California

2·3·29·67 Posts
Default

Quote:
Originally Posted by bentonsar View Post
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
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.
ewmayer is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
Is there utility in implementing an LLT-like sequence for GIMPS? kriesel Math 11 2019-01-21 18:16
Implementing Factoring Algorithms Raman Hobbies 45 2009-05-11 05:11
Implementing MPQS: SOS! smoking81 Factoring 10 2007-10-02 12:30
Implementing algorithms, did I do this right? ShiningArcanine Programming 18 2005-12-29 21:47
improvement needed in SSE2 for Pentium 4 for non-IBDWT case for LLR/PRP wfgarnett3 Hardware 3 2004-06-29 03:03

All times are UTC. The time now is 07:14.


Thu Oct 21 07:14:58 UTC 2021 up 90 days, 1:43, 1 user, load averages: 1.58, 1.42, 1.28

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.