cudadl/parrhoasm.cu

/*
 *  cudadl version 0.9: Compute discrete logs in smooth group orders
 *                      using CUDA
 *  Copyright (C) 2012 by Ryan Henry and Ian Goldberg
 *                        {rhenry,iang}@cs.uwaterloo.ca
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of version 3 of the GNU General Public License as
 *  published by the Free Software Foundation.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "cudadl.h"
#include <fstream>
#include <NTL/vec_ZZ.h>
#include <string.h>
#include <limits.h>
#include <sys/time.h>
#include <time.h>
#include <string>
#include <sstream>
#include <utility>
#include <map>

// #define CHECK_RESULTS

// #include "cuPrintf.cu"

#define BITS_PER_WORD (8 * sizeof(unsigned int))
#define TWO_32_DIV_3 1431655765U // floor( 2^32 / 3 )
#define TWO_32_DIV_3_X2 2863311530U // 2 * floor( 2^32 / 3 )
#define TWO_64_DIV_3 6148914691236517205UL // floor( 2^64 / 3 )

#define SIZE_G WORDS
#define SIZE_Y WORDS
#define SIZE_X (WORDS + 1)

#include "dpstream.cu"

// #define X(idx) x[idx<<5]
// #define Y(idx) (multtype == 0 ? c_g[idx] : multtype == 1 ? c_y[idx] : x[idx<<5])
// #define Y(idx) ((c_g[idx]&typemask0) | (c_y[idx]&typemask1) | (x[idx<<5]&typemask2))
// #define Z(idx) z[idx<<5]

NTL_CLIENT

#ifdef CHECK_RESULTS
static void dump(const char *prefix, const unsigned int *di, size_t words)
{
    size_t l = words * 4;
    const unsigned char *d = (const unsigned char *) di;
    d += l;
    printf("%s=", prefix);
    while(l)
    {
        --d;
        printf("%02X", *d);
        --l;
    }
    printf("\n");
}
#endif

static void checkCUDAError(const char *msg)
{
    cudaError_t err = cudaGetLastError();
    if (cudaSuccess != err) 
    {
        fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString(err));
        exit(EXIT_FAILURE);
    }                         
}

typedef struct
{
    unsigned int x[SIZE_X];
    unsigned int a[3], b[3];
} GlobalThreadState;


__device__ __constant__ unsigned int c_rho[WORDS + 1], c_r_inv[WORDS + 1], c_rho_prime;
__device__ __constant__ unsigned int c_y[SIZE_Y], c_g[SIZE_G];

/*
__device__ unsigned int getgy(int i)
{
    return c_rho[i] + c_y[i] + c_g[i];
}
*/

#include "cios.asm"

#if 0
__device__ void _sub(unsigned int *x)
// x <- x - y
{
    asm("sub.cc.u32  %0, %1, %2;" : "=r"(X(0))  : "r"(X(0)),  "r"(c_rho[0]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(1))  : "r"(X(1)),  "r"(c_rho[1]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(2))  : "r"(X(2)),  "r"(c_rho[2]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(3))  : "r"(X(3)),  "r"(c_rho[3]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(4))  : "r"(X(4)),  "r"(c_rho[4]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(5))  : "r"(X(5)),  "r"(c_rho[5]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(6))  : "r"(X(6)),  "r"(c_rho[6]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(7))  : "r"(X(7)),  "r"(c_rho[7]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(8))  : "r"(X(8)),  "r"(c_rho[8]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(9))  : "r"(X(9)),  "r"(c_rho[9]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(10)) : "r"(X(10)), "r"(c_rho[10]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(11)) : "r"(X(11)), "r"(c_rho[11]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(12)) : "r"(X(12)), "r"(c_rho[12]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(13)) : "r"(X(13)), "r"(c_rho[13]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(14)) : "r"(X(14)), "r"(c_rho[14]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(15)) : "r"(X(15)), "r"(c_rho[15]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(16)) : "r"(X(16)), "r"(c_rho[16]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(17)) : "r"(X(17)), "r"(c_rho[17]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(18)) : "r"(X(18)), "r"(c_rho[18]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(19)) : "r"(X(19)), "r"(c_rho[19]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(20)) : "r"(X(20)), "r"(c_rho[20]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(21)) : "r"(X(21)), "r"(c_rho[21]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(22)) : "r"(X(22)), "r"(c_rho[22]));
    asm("subc.cc.u32 %0, %1, %2;" : "=r"(X(23)) : "r"(X(23)), "r"(c_rho[23]));
    asm("subc.u32    %0, %1, %2;" : "=r"(X(24)) : "r"(X(24)), "r"(c_rho[24]));
}

__device__ inline bool _gt(const unsigned int * x)
// returns true iff x > y  ;  x is of length WORDS+1 and y is of length WORDS
{
    if (X(24)) return true;
    if (X(23) < c_rho[23]) return false;
    if (X(23) > c_rho[23]) return true;
    if (X(22) < c_rho[22]) return false;
    if (X(22) > c_rho[22]) return true;
    if (X(21) < c_rho[21]) return false;
    if (X(21) > c_rho[21]) return true;
    if (X(20) < c_rho[20]) return false;
    if (X(20) > c_rho[20]) return true;
    if (X(19) < c_rho[19]) return false;
    if (X(19) > c_rho[19]) return true;
    if (X(18) < c_rho[18]) return false;
    if (X(18) > c_rho[18]) return true;
    if (X(17) < c_rho[17]) return false;
    if (X(17) > c_rho[17]) return true;
    if (X(16) < c_rho[16]) return false;
    if (X(16) > c_rho[16]) return true;
    if (X(15) < c_rho[15]) return false;
    if (X(15) > c_rho[15]) return true;
    if (X(14) < c_rho[14]) return false;
    if (X(14) > c_rho[14]) return true;
    if (X(13) < c_rho[13]) return false;
    if (X(13) > c_rho[13]) return true;
    if (X(12) < c_rho[12]) return false;
    if (X(12) > c_rho[12]) return true;
    if (X(11) < c_rho[11]) return false;
    if (X(11) > c_rho[11]) return true;
    if (X(10) < c_rho[10]) return false;
    if (X(10) > c_rho[10]) return true;
    if (X(9)  < c_rho[9])  return false;
    if (X(9)  > c_rho[9])  return true;
    if (X(8)  < c_rho[8])  return false;
    if (X(8)  > c_rho[8])  return true;
    if (X(7)  < c_rho[7])  return false;
    if (X(7)  > c_rho[7])  return true;
    if (X(6)  < c_rho[6])  return false;
    if (X(6)  > c_rho[6])  return true;
    if (X(5)  < c_rho[5])  return false;
    if (X(5)  > c_rho[5])  return true;
    if (X(4)  < c_rho[4])  return false;
    if (X(4)  > c_rho[4])  return true;
    if (X(3)  < c_rho[3])  return false;
    if (X(3)  > c_rho[3])  return true;
    if (X(2)  < c_rho[2])  return false;
    if (X(2)  > c_rho[2])  return true;
    if (X(1)  < c_rho[1])  return false;
    if (X(1)  > c_rho[1])  return true;
    return (X(0) > c_rho[0]);
}
#endif

#define nmult 1000
#define nlaunch 2

#define threadDimx 32

__global__ void cudaMulmod(GlobalThreadState *global_ts,
	unsigned int order_0, unsigned int order_1, unsigned int order_2)
// x <- x * y mod rho
{
    // extern __shared__ unsigned int sharedmem[];
    register int i;
    register int tid = threadIdx.x + threadDimx * threadIdx.y;
    register int global_tid = tid + blockDim.x * blockDim.y * blockIdx.x;
    unsigned int a_0, a_1, a_2;
    unsigned int b_0, b_1, b_2;
    unsigned int *global_x_base = global_ts[global_tid].x;
    unsigned int *global_a_base = global_ts[global_tid].a;
    unsigned int *global_b_base = global_ts[global_tid].b;
    // unsigned int *shared_table_base = sharedmem;
    // unsigned int *shared_x_base = shared_table_base + threadIdx.y * 32 * (SIZE_X) + threadIdx.x;
    unsigned int xlow;

    // cuPrintf("d_z = %08X%08X%08X\n", global_x_base[2], global_x_base[1], global_x_base[0]);
    loadx(global_x_base);
    a_2 = global_a_base[2];
    a_1 = global_a_base[1];
    a_0 = global_a_base[0];
    b_2 = global_b_base[2];
    b_1 = global_b_base[1];
    b_0 = global_b_base[0];
    // cuPrintf("s_A = %08X\n", (unsigned int)a);
    // cuPrintf("s_B = %08X\n", (unsigned int)b);
    // cuPrintf("s_a = %08X%08X\n", global_a_base[1], global_a_base[0]);
    // cuPrintf("s_b = %08X%08X\n", global_b_base[1], global_b_base[0]);
    for (i = 0; i < nmult; ++i)
    {
        //memset(ts[tid].z, 0, (WORDS + 2) * sizeof(unsigned int));
        //__syncthreads();
        //cuPrintf("x    = %p %08X%08X%08X\n", ts[tid].x, ts[tid].x[2], ts[tid].x[1], ts[tid].x[0]);
        //cuPrintf("y    = %p %08X%08X\n", c_y, c_y[1], c_y[0]);
        //cuPrintf("z    = %p %08X%08X%08X%08X\n", ts[tid].z, ts[tid].z[3], ts[tid].z[2],  ts[tid].z[1], ts[tid].z[0]);
        int multtype;
	asm("mov.u32 %0, $xr0;" : "=r" (xlow));
        if (xlow < TWO_32_DIV_3)
        {
            multtype = 0;
	    asm("add.cc.u32 %0, %1, %2;" : "=r"(a_0) : "r"(a_0), "r"((unsigned int)1U));
	    asm("addc.cc.u32 %0, %1, %2;" : "=r"(a_1) : "r"(a_1), "r"((unsigned int)0U));
	    asm("addc.u32 %0, %1, %2;" : "=r"(a_2) : "r"(a_2), "r"((unsigned int)0U));
//	    a += 1;
	    // cuPrintf("inc a\n");
        }
        else if (xlow < TWO_32_DIV_3_X2)
        {
            multtype = 1;
	    asm("add.cc.u32 %0, %1, %2;" : "=r"(b_0) : "r"(b_0), "r"((unsigned int)1U));
	    asm("addc.cc.u32 %0, %1, %2;" : "=r"(b_1) : "r"(b_1), "r"((unsigned int)0U));
	    asm("addc.u32 %0, %1, %2;" : "=r"(b_2) : "r"(b_2), "r"((unsigned int)0U));
//	    b += 1;
	    // cuPrintf("inc b\n");
        }
        else
        {
            multtype = 2;
	    asm("add.cc.u32 %0, %1, %2;" : "=r"(a_0) : "r"(a_0), "r"(a_0));
	    asm("addc.cc.u32 %0, %1, %2;" : "=r"(a_1) : "r"(a_1), "r"(a_1));
	    asm("addc.u32 %0, %1, %2;" : "=r"(a_2) : "r"(a_2), "r"(a_2));
	    asm("add.cc.u32 %0, %1, %2;" : "=r"(b_0) : "r"(b_0), "r"(b_0));
	    asm("addc.cc.u32 %0, %1, %2;" : "=r"(b_1) : "r"(b_1), "r"(b_1));
	    asm("addc.u32 %0, %1, %2;" : "=r"(b_2) : "r"(b_2), "r"(b_2));
//	    a += a;
//	    b += b;
	    // cuPrintf("double\n");
        }
	if (a_2 > order_2 || ((a_2 == order_2) && (a_1 > order_1)) || (((a_2 == order_2) && (a_1 == order_1) && (a_0 > order_0)))) {
	    asm("sub.cc.u32 %0, %1, %2;" : "=r"(a_0) : "r"(a_0), "r"(order_0));
	    asm("subc.cc.u32 %0, %1, %2;" : "=r"(a_1) : "r"(a_1), "r"(order_1));
	    asm("subc.u32 %0, %1, %2;" : "=r"(a_2) : "r"(a_2), "r"(order_2));
	}
	if (b_2 > order_2 || ((b_2 == order_2) && (b_1 > order_1)) || (((b_2 == order_2) && (b_1 == order_1) && (b_0 > order_0)))) {
	    asm("sub.cc.u32 %0, %1, %2;" : "=r"(b_0) : "r"(b_0), "r"(order_0));
	    asm("subc.cc.u32 %0, %1, %2;" : "=r"(b_1) : "r"(b_1), "r"(order_1));
	    asm("subc.u32 %0, %1, %2;" : "=r"(b_2) : "r"(b_2), "r"(order_2));
	}
		
        modmul(multtype);
        //memcpy(ts[tid].x, ts[tid].z, (WORDS + 1) * sizeof(unsigned int));
	/*
        if (_gt(shared_x_base)) {
            _sub(shared_x_base);
	}
	*/
        //cuPrintf("y    = %p %08X%08X\n", c_y, c_y[1], c_y[0]);
        //cuPrintf("z    = %p %08X%08X\n", ts[tid].x, ts[tid].x[1], ts[tid].x[0]);
        //cuPrintf("y    = %p %08X%08X\n", c_y, c_y[1], c_y[0]);

	// Check for a distinguished point
	asm("mov.u32 %0, $xr0;" : "=r" (xlow));
	if ((xlow & 0x000003ff) == 0 || !(order_2|order_1) && order_0 < (1<<20)) {
	    DPstreamWrite(a_0,a_1,a_2,b_0,b_1,b_2);
	}
    }
    savex(global_x_base);
    global_a_base[0] = a_0;
    global_a_base[1] = a_1;
    global_a_base[2] = a_2;
    global_b_base[0] = b_0;
    global_b_base[1] = b_1;
    global_b_base[2] = b_2;
    // cuPrintf("d_z = %08X%08X%08X\n", global_x_base[2], global_x_base[1], global_x_base[0]);
    // cuPrintf("d_A = %08X\n", (unsigned int)a);
    // cuPrintf("d_B = %08X\n", (unsigned int)b);
    // cuPrintf("d_a = %08X%08X\n", global_a_base[1], global_a_base[0]);
    // cuPrintf("d_b = %08X%08X\n", global_b_base[1], global_b_base[0]);
    if (!(order_0|order_1|order_2)) evilhack();
}

int nthreads = 25600;
int nblocks = 50;

typedef map<std::string, pair<ZZ,ZZ> > DTable;

struct CBData {
    const ZZ_p &base;
    const ZZ_p &target;
    const ZZ &order;
    unsigned long long numdp;
    DTable dtable;
    bool found_collision;
    ZZ expon;

    CBData(const ZZ_p &_base, const ZZ_p &_target, const ZZ &_order) :
	base(_base), target(_target), order(_order), numdp(0),
	found_collision(false) {}
};

static void dpcallback(void *cbdata, unsigned short threadId,
    unsigned short blockId, string x, unsigned int a_0, unsigned int a_1,
    unsigned int a_2, unsigned int b_0, unsigned int b_1, unsigned int b_2)
{
    CBData *d = (CBData*)cbdata;
    // WARNING: this assumes
    //    sizeof(unsigned long) == sizeof(unsigned long long) !
    ZZ zz_a = to_ZZ(a_2);
    zz_a <<= 32;
    zz_a += a_1;
    zz_a <<= 32;
    zz_a += a_0;
    ZZ zz_b = to_ZZ(b_2);
    zz_b <<= 32;
    zz_b += b_1;
    zz_b <<= 32;
    zz_b += b_0;

    //ZZ_p dp = power(d->base, zz_a) * power(d->target, zz_b);

    // cerr << "DP " << ++(d->numdp) << "\r";
    pair<ZZ,ZZ> ab(zz_a,zz_b);
    pair<DTable::iterator, bool> res;
    res = d->dtable.insert(DTable::value_type(x, ab));
    if (!res.second) {
	// Collision!
	ZZ adiff = to_ZZ(res.first->second.first) - zz_a;
	ZZ bdiff = zz_b - to_ZZ(res.first->second.second);
	if (bdiff < 0) bdiff += d->order;
	ZZ binv;
	if (InvModStatus(binv, bdiff, d->order) == 0) {
	    d->expon = MulMod(binv, adiff, d->order);
	    d->found_collision = true;
	} else {
	    if (d->order > (1<<20)) {
		cerr << "Unhelpful collision\n";
	    }
	}
    }
}

ZZ cuda_dl(const ZZ_p &base, const ZZ_p &target, const ZZ &order,
		const ZZ &modulus)
{
    unsigned long long totmicros = 0;
    ZZ_pBak pbak;
    pbak.save();

    if (NumBits(modulus) <= ((WORDS-1)*BITS_PER_WORD) ||
	    NumBits(modulus) > WORDS*BITS_PER_WORD) {
	cerr << "modulus is not " << WORDS << " words long.\n";
	exit(1);
    }

    long orderbits = NumBits(order);
    if (orderbits > 92) {
	cerr << "order is larger than 92 bits.\n";
	exit(1);
    }

	unsigned int order_2 = trunc_long(order >> 64, 32);
	unsigned int order_1 = trunc_long(order >> 32, 32);
	unsigned int order_0 = trunc_long(order, 32);
    /*
    cerr << "order   = " << order << "\n";
    cerr << "orderll = " << orderll << "\n";
    */

    unsigned int *rho, rho_prime, *r, *r_inv;
    unsigned int l_y[SIZE_Y], l_g[SIZE_G];

//cout << "rho = " << zz_rho << "\n";
    ZZ zz_r = (to_ZZ(1) << (WORDS * BITS_PER_WORD));
    ZZ_p::init(zz_r);
    ZZ_p zz_p_rho_prime = -to_ZZ_p(1) / to_ZZ_p(modulus);
//cout << "rho_prime = " << zz_p_rho_prime << "\n";
    ZZ_p::init(modulus);
    ZZ_p zz_p_r, zz_p_r_inv;
    zz_p_r = to_ZZ_p(zz_r);
    zz_p_r_inv = to_ZZ_p(1) / zz_p_r;

    CBData cbdata(base, target, order);
//cout << "r = " << zz_r << "\n";
//cout << "r mod rho = " << rep(zz_p_r) << "\n";
//cout << "r_inv = " << rep(zz_p_r_inv) << "\n";
//cout << "r * r_inv = " << zz_p_r * zz_p_r_inv << "\n";
//cout << "r * r_inv - rho * rho_prime = " << zz_r * rep(zz_p_r_inv) - zz_rho * rep(zz_p_rho_prime) << "\n";
    rho = (unsigned int*) calloc((WORDS + 1), sizeof(unsigned int));
    r = (unsigned int*) calloc((WORDS + 1), sizeof(unsigned int));
    r_inv = (unsigned int*) calloc((WORDS + 1), sizeof(unsigned int));
    BytesFromZZ((unsigned char *) rho, modulus, WORDS * sizeof(unsigned int));
    BytesFromZZ((unsigned char *) &rho_prime, rep(zz_p_rho_prime), sizeof(unsigned int));
    BytesFromZZ((unsigned char *) r,  rep(zz_p_r),  (WORDS + 1) * sizeof(unsigned int));
    BytesFromZZ((unsigned char *) r_inv,  rep(zz_p_r_inv),  (WORDS + 1) * sizeof(unsigned int));
    cudaMemcpyToSymbol(c_rho, rho, (WORDS + 1) * sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
    cudaMemcpyToSymbol(c_rho_prime, &rho_prime, sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
    cudaMemcpyToSymbol(c_r_inv, r_inv, (WORDS + 1) * sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
    checkCUDAError("memcpytosymbol");

/*
    if (nthreads % threadDimx) {
	nthreads /= threadDimx;
	nthreads *= threadDimx;
	cerr << "Truncating to " << nthreads << " threads.\n";
    }
*/

    if (nthreads % nblocks) {
	cerr << "Error: " << nthreads << " not a multiple of " << nblocks << "\n";
	exit(1);
    }
    const int threadsPerBlock = nthreads / nblocks;
    GlobalThreadState *d_ts;
    cudaMalloc((void **) &d_ts, nthreads * sizeof(GlobalThreadState));
    GlobalThreadState *l_ts;
    l_ts = (GlobalThreadState *) calloc(nthreads, sizeof(GlobalThreadState));
    unsigned int *l_z;
    l_z = (unsigned int *) calloc(nthreads, WORDS * sizeof(unsigned int));

    //cudaPrintfInit((1<<20)*16);
    DPstreamInit(1<<16);
    struct timeval st, et;

    memset(l_ts, 0, nthreads * sizeof(GlobalThreadState));
    memset(l_z, 0, WORDS * sizeof(unsigned int));

//cout << "y = " << rep(zz_p_y) << "\n";
    ZZ_p zz_p_yr = target * zz_p_r;
//cout << "y = " << rep(zz_p_y) << "\n";
    ZZ_p zz_p_gr = base * zz_p_r;
    BytesFromZZ((unsigned char *) l_y, rep(zz_p_yr), WORDS * sizeof(unsigned int));
    BytesFromZZ((unsigned char *) l_g, rep(zz_p_gr), WORDS * sizeof(unsigned int));
    cudaMemcpyToSymbol(c_y, l_y, (WORDS) * sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
    cudaMemcpyToSymbol(c_g, l_g, (WORDS) * sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
    //dump("y", l_y, WORDS);

#ifdef CHECK_RESULTS
    int totmult = 0;
    bool fail = false;
    vec_ZZ aexp, bexp;  // The expected values of a and b at completion
    aexp.SetLength(nthreads);
    bexp.SetLength(nthreads);
#endif
    ZZ a = RandomBits_ZZ(orderbits-1);
    ZZ b = RandomBits_ZZ(orderbits-1);
    ZZ astep = RandomBits_ZZ(orderbits-1);
    ZZ bstep = RandomBits_ZZ(orderbits-1);
    ZZ_p zz_p_x = power(base, a) * power(target, b);
    ZZ_p zz_p_step = power(base, astep) * power(target, bstep);
    ZZ_p zz_p_xr = zz_p_x * zz_p_r;
    for (int t=0; t<nthreads; ++t)
    {
	zz_p_xr *= zz_p_step;
	a += astep;
	a %= order;
	b += bstep;
	b %= order;
	BytesFromZZ((unsigned char *) l_ts[t].x, rep(zz_p_xr), WORDS * sizeof(unsigned int));
	// dump("l_x", l_ts[t].x, WORDS);
	BytesFromZZ((unsigned char *) l_ts[t].a, a, 3 * sizeof(unsigned int));
	BytesFromZZ((unsigned char *) l_ts[t].b, b, 3 * sizeof(unsigned int));
	//dump("x", l_ts[t].x, WORDS);
//cout << "x = " << rep(zz_p_x) << "\n";
//cout << "xr = " << rep(zz_p_xr) << "\n";
	//dump("rho", rho, WORDS);
#ifdef CHECK_RESULTS
	for (int l = 0; l < nmult * nlaunch; ++l)
	{
	    unsigned int w = trunc_long(rep(zz_p_xr), 32);
	    if (w < TWO_32_DIV_3) {
		zz_p_xr *= base;
		a++;
	    } else if (w < TWO_32_DIV_3_X2) {
		zz_p_xr *= target;
		b++;
	    } else {
		zz_p_xr *= zz_p_xr * zz_p_r_inv;
		a += a;
		b += b;
	    }
	    if (a > order) a -= order;
	    if (b > order) b -= order;
	    ++totmult;
	}     
	BytesFromZZ((unsigned char *) (l_z + t * WORDS), rep(zz_p_xr), WORDS * sizeof(unsigned int));
	aexp[t] = a;
	bexp[t] = b;
#endif
//cout << "zr = " << rep(zz_p_xr) << "\n";
	//dump("z", l_z + t * WORDS, WORDS);
    }

    cudaMemcpy(d_ts, l_ts, nthreads * sizeof(GlobalThreadState), cudaMemcpyHostToDevice);
    dim3 tpb(threadDimx, threadsPerBlock/threadDimx);
    gettimeofday(&st, NULL);
    int launchcount = 0;
#ifdef CHECK_RESULTS
    for (int ln=0; ln<nlaunch; ++ln)
#else
    while(cbdata.found_collision == false)
#endif
    {
	cerr << getpid() << " Launch " << ++launchcount << "...\n";
	cudaMulmod<<< nblocks, nthreads/nblocks /*tpb*/, 0 >>>(d_ts, order_0, order_1, order_2);

	cudaThreadSynchronize();
	checkCUDAError("kernel launch");
	DPstreamParse(&cbdata);
	//cudaPrintfExtractDPE(dpcallback, &cbdata);
	//cudaPrintfDisplay(stdout, true);
	cerr << getpid() << "\n";
    }
    gettimeofday(&et, NULL);

    totmicros = (unsigned long long)(et.tv_sec - st.tv_sec) * 1000000 + (et.tv_usec - st.tv_usec);
#ifdef CHECK_RESULTS
    cudaMemcpy(l_ts, d_ts, nthreads * sizeof(GlobalThreadState), cudaMemcpyDeviceToHost);
    checkCUDAError("memcpy");

    int j = 0;
    for (int t = 0; t < nthreads; ++t)
    {
	unsigned int *l_Z = l_ts[t].x;
	for (int i = 0; i < WORDS; i++)
	{
	    if (l_Z[i] != l_z[i + t * WORDS])
	    {
		fail = true;
cout << i << ": " << l_Z[i] << " != " << l_z[i + t * WORDS] << "\n";
	    }
	}
	if (fail)
	{
	    dump("d_z", l_Z, WORDS);
	    dump("l_z", l_z + t * WORDS, WORDS);
	}
	ZZ ares = (to_ZZ(l_ts[t].a[1]) << 32) + to_ZZ(l_ts[t].a[0]);
	ZZ bres = (to_ZZ(l_ts[t].b[1]) << 32) + to_ZZ(l_ts[t].b[0]);
	if (ares != aexp[t] || bres != bexp[t]) {
	    cerr << "ares = " << ares << "\n";
	    cerr << "aexp = " << aexp[t] << "\n";
	    cerr << "bres = " << bres << "\n";
	    cerr << "bexp = " << bexp[t] << "\n";
	}
	++j;
    }
    if (!fail) {
	cerr << "Results correct.\n";
    }
#endif

    unsigned long long totnanos = totmicros * 1000;
    cout << totmicros << " us / " << nthreads << " = " << totmicros / nthreads << " us / " << (nmult*launchcount) << " = " << totnanos / ((unsigned long long)nthreads * nmult * launchcount) << " ns\n";
    //cudaPrintfEnd();
    DPstreamEnd();

    cudaFree(d_ts);
    free(rho);
    free(r_inv);
    free(r);
    free(l_ts);
    free(l_z);

    pbak.restore();

    return cbdata.expon;
}

#ifdef TEST_CUDA
int main(int argc, char** argv)
{
#ifdef DERANDOMIZE
    SetSeed(to_ZZ(1));
#else
    // Initialize the prng with some randomness from the kernel
    unsigned char randbuf[1024];
    ifstream urand("/dev/urandom");
    urand.read((char *) randbuf, sizeof(randbuf));
    urand.close();
    ZZ randzz = ZZFromBytes(randbuf, sizeof(randbuf));
    SetSeed(randzz);
#endif

    ZZ rho, p, q;
    vec_ZZ pfvec, qfvec;
    cin >> rho >> p >> pfvec >> q >> qfvec;

    if (argc > 1) nthreads = atoi(argv[1]);
    if (argc > 2) nblocks = atoi(argv[2]);

    ZZ_p::init(p);

    for (int iter = 0; iter < 1; ++iter) {
	// Create the subproblem for the f'th factor of p-1
	int f = iter;

	ZZ remorder = (p-1)/pfvec[f];
	ZZ_p g = power(to_ZZ_p(2), remorder);
	ZZ_p y = power(random_ZZ_p(), remorder);

	// g should now be of order pfvec[0].  Check that.
	if (g == 1 || power(g, pfvec[f]) != 1) {
	    cerr << "base has the wrong order!\n";
	    exit(1);
	}

	// Try to find the DL_g of y
	cerr << "DL_" << g << "(" << y << ") mod " << p << "\n";
	ZZ e = cuda_dl(g, y, pfvec[f], p);
	cerr << "e = " << e << "\n";
	cerr << ( (power(g,e) == y) ? "CORRECT!" : "INCORRECT!" ) << "\n";
    }

    return 0;
}
#endif