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 <NTL/vec_ZZ_p.h>
#include <string.h>
#include <limits.h>
#include <sys/time.h>
#include <time.h>
#include <unistd.h>
#include <string>
#include <sstream>
#include <utility>
#include <map>

#include "atomic_iostream.h"
#include "worker.h"

// #define CHECK_RESULTS

 #define VERBOSE

// #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)

__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];

#include "cios.asm"

#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__ unsigned int getgy(int i)
{
    return c_rho[i] + c_y[i] + c_g[i];
}
*/

#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,
	unsigned int dpfreq)
// x <- x * y mod rho
{
    CIOS_REG_DECLARE;

    // 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]);

    CIOS_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]);

    #pragma unroll 1 // to prevent namespace conflicts caused by loop unrolling when 'nmult' is small
    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 volatile("mov.u32 %0, $xr0;" : "=r" (xlow));
	// volatile qualifier needed to prevent the compiler from moving it outside of the loop
	// change can be seen by building with 'make ptx' and reading the 'parrhoasm.ptx' file

        if (xlow < TWO_32_DIV_3)
        {
            multtype = 0;
	    asm volatile("add.cc.u32 %0, %1, %2;" : "=r"(a_0) : "r"(a_0), "r"((unsigned int)1U));
	    asm volatile("addc.cc.u32 %0, %1, %2;" : "=r"(a_1) : "r"(a_1), "r"((unsigned int)0U));
	    asm volatile("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 volatile("add.cc.u32 %0, %1, %2;" : "=r"(b_0) : "r"(b_0), "r"((unsigned int)1U));
	    asm volatile("addc.cc.u32 %0, %1, %2;" : "=r"(b_1) : "r"(b_1), "r"((unsigned int)0U));
	    asm volatile("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 volatile("add.cc.u32 %0, %1, %2;" : "=r"(a_0) : "r"(a_0), "r"(a_0));
	    asm volatile("addc.cc.u32 %0, %1, %2;" : "=r"(a_1) : "r"(a_1), "r"(a_1));
	    asm volatile("addc.u32 %0, %1, %2;" : "=r"(a_2) : "r"(a_2), "r"(a_2));
	    asm volatile("add.cc.u32 %0, %1, %2;" : "=r"(b_0) : "r"(b_0), "r"(b_0));
	    asm volatile("addc.cc.u32 %0, %1, %2;" : "=r"(b_1) : "r"(b_1), "r"(b_1));
	    asm volatile("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 volatile("sub.cc.u32 %0, %1, %2;" : "=r"(a_0) : "r"(a_0), "r"(order_0));
	    asm volatile("subc.cc.u32 %0, %1, %2;" : "=r"(a_1) : "r"(a_1), "r"(order_1));
	    asm volatile("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 volatile("sub.cc.u32 %0, %1, %2;" : "=r"(b_0) : "r"(b_0), "r"(order_0));
	    asm volatile("subc.cc.u32 %0, %1, %2;" : "=r"(b_1) : "r"(b_1), "r"(order_1));
	    asm volatile("subc.u32 %0, %1, %2;" : "=r"(b_2) : "r"(b_2), "r"(order_2));
	}
		
        CIOS_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 volatile("mov.u32 %0, $xr0;" : "=r" (xlow));
	if (xlow <= dpfreq) {
	    unsigned int *ourbuffer = DPstreamAlloc();
	    if (ourbuffer) {
		CIOS_WRITE_DP(ourbuffer, a_0, a_1, a_2, b_0, b_1, b_2);
	    }
	}
    }

    CIOS_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 defined(NTHREADS) || defined(NBLOCKS)
// must define both together
int nthreads = NTHREADS;
int nblocks = NBLOCKS;
#else
int nthreads = 25600;
int nblocks = 50;
#endif

void cuda_dl(const ZZ_p &base, const ZZ_p &target, const ZZ &order,
		const ZZ &modulus, unsigned int dpfreq, void *cbdata,
		unsigned int* launch_count, bool* filled_dp_buffer)
{
    unsigned long long totmicros = 0;
    ZZ_pBak pbak;
    pbak.save();

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

    long orderbits = NumBits(order);
    if (orderbits > 92) {
	AtomicWriter(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;

//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) {
	AtomicWriter(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
    vec_ZZ astart, bstart;  // the starting values of a and b
    vec_ZZ_p zz_p_xr_start;
    astart.SetLength(nthreads);
    bstart.SetLength(nthreads);
    zz_p_xr_start.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));
#ifdef CHECK_RESULTS
	astart[t] = a;
	bstart[t] = b;
	zz_p_xr_start[t] = zz_p_xr;
#endif
	//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);
    }

    cudaMemcpy(d_ts, l_ts, nthreads * sizeof(GlobalThreadState), cudaMemcpyHostToDevice);
    dim3 tpb(threadDimx, threadsPerBlock/threadDimx);
    gettimeofday(&st, NULL);
    int local_launchcount = 0;
    bool stop_computing = false;
    unsigned int num_dps = 0;

#ifdef VERBOSE
    AtomicWriter(cerr) << getpid() << " Using (nblocks, nthreads): " << nblocks << ", " << nthreads << "\n";
#endif

#ifdef CHECK_RESULTS
    for (int ln=0; ln<nlaunch && stop_computing == false; ++ln)
#else
    while(stop_computing == false)
#endif
    {
	local_launchcount++;
#ifdef VERBOSE
	AtomicWriter(cerr) << getpid() << " Launch " << local_launchcount << "...\n";
#endif
	cudaMulmod<<< nblocks, nthreads/nblocks /*tpb*/, 0 >>>(d_ts,
	    order_0, order_1, order_2, dpfreq);
	checkCUDAError("kernel launch");

	cudaThreadSynchronize();
	checkCUDAError("thread sync");
	bool single_launch_filled_dp_buffer = false;
	stop_computing = DPstreamParse(cbdata, &num_dps, &single_launch_filled_dp_buffer);
	*filled_dp_buffer = (*filled_dp_buffer || single_launch_filled_dp_buffer);
#ifdef VERBOSE
	AtomicWriter(cerr) << getpid() << " Finished " << local_launchcount << " (" << num_dps << " DPs)\n";
#endif
	//cudaPrintfExtractDPE(dpcallback, cbdata);
	//cudaPrintfDisplay(stdout, true);
    }
    gettimeofday(&et, NULL);

    totmicros = (unsigned long long)(et.tv_sec - st.tv_sec) * 1000000 + (et.tv_usec - st.tv_usec);

#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);
    ZZ atemp;
    ZZ btemp;
    ZZ_p zz_p_xr_temp;

    for (int t=0; t<nthreads; ++t)
    {
	atemp = astart[t];
	btemp = bstart[t];
	zz_p_xr_temp = zz_p_xr_start[t];

	// local_launchcount should be <= nlaunch
	for (int l = 0; l < nmult * local_launchcount; ++l)
	{
	    unsigned int w = trunc_long(rep(zz_p_xr_temp), 32);
	    if (w < TWO_32_DIV_3) {
		zz_p_xr_temp *= base;
		atemp++;
	    } else if (w < TWO_32_DIV_3_X2) {
		zz_p_xr_temp *= target;
		btemp++;
	    } else {
		zz_p_xr_temp *= zz_p_xr_temp * zz_p_r_inv;
		atemp += atemp;
		btemp += btemp;
	    }
	    if (atemp > order) atemp -= order;
	    if (btemp > order) btemp -= order;
	    ++totmult;
	}
	BytesFromZZ((unsigned char *) (l_z + t * WORDS), rep(zz_p_xr_temp), WORDS * sizeof(unsigned int));
	aexp[t] = atemp;
	bexp[t] = btemp;
	//cout << "zr = " << rep(zz_p_xr_temp) << "\n";
	//dump("z", l_z + t * WORDS, WORDS);
    }
    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[2]) << 64) + (to_ZZ(l_ts[t].a[1]) << 32) + to_ZZ(l_ts[t].a[0]);
	ZZ bres = (to_ZZ(l_ts[t].b[2]) << 64) + (to_ZZ(l_ts[t].b[1]) << 32) + to_ZZ(l_ts[t].b[0]);
	if (ares != aexp[t] || bres != bexp[t]) {
	    AtomicWriter atomic_cerr(cerr);
	    atomic_cerr << "ares = " << ares << "\n";
	    atomic_cerr << "aexp = " << aexp[t] << "\n";
	    atomic_cerr << "bres = " << bres << "\n";
	    atomic_cerr << "bexp = " << bexp[t] << "\n";
	}
	++j;
    }
    if (!fail) {
	AtomicWriter(cerr) << "Results correct.\n";
    }
#endif

    char timestamp[20];
    sprintf(timestamp, "%d.%06d", et.tv_sec, et.tv_usec);

    unsigned long long totnanos = totmicros * 1000;
    AtomicWriter(cout) << timestamp << ":" << output_prefix << ": " << totmicros << " us / " << nthreads << " = " << totmicros / nthreads << " us / " << (nmult*local_launchcount) << " = " << totnanos / ((unsigned long long)nthreads * nmult * local_launchcount) << " ns\n";
    cout.flush();
    //cudaPrintfEnd();
    DPstreamEnd();

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

    pbak.restore();

    *launch_count = local_launchcount;
}

#ifdef TEST_CUDA

string output_prefix("test_cuda");

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) {}
};

// This function is called from inside cuda_dl for each DP it encounters.
// It calls the function named "dpcallback" directly.  It would be
// cleaner if this were passed as a function pointer to cuda_dl, but
// that makes nvcc 3.1 segfault. :-p
// dp points to an array of WORDS+7 unsigned ints:
// - 1 word of threadID/blockID
// - WORDS words of the dp value
// - 3 words of a
// - 3 words of b
bool dpcallback(void *cbdata, unsigned int *dpwords)
{
    CBData *d = (CBData*)cbdata;
    // WARNING: this assumes
    //    sizeof(unsigned long) == sizeof(unsigned long long) !
    ZZ zz_a = to_ZZ(dpwords[WORDS+3]);
    zz_a <<= 32;
    zz_a += dpwords[WORDS+2];
    zz_a <<= 32;
    zz_a += dpwords[WORDS+1];
    ZZ zz_b = to_ZZ(dpwords[WORDS+6]);
    zz_b <<= 32;
    zz_b += dpwords[WORDS+5];
    zz_b <<= 32;
    zz_b += dpwords[WORDS+4];

    // ZZ_p dp = power(d->base, zz_a) * power(d->target, zz_b);
    // cerr << zz_a << " " << zz_b << " " << dp << "\n";

    // cerr << "DP " << ++(d->numdp) << "\r";
    pair<ZZ,ZZ> ab(zz_a,zz_b);
    pair<DTable::iterator, bool> res;
    string x((const char *)(dpwords+1), WORDS*sizeof(unsigned int));

    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";
	    }
	}
    }

    return d->found_collision;
}

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";
	CBData cbdata(g, y, pfvec[f]);
	unsigned int* launch_count;
	cuda_dl(g, y, pfvec[f], p, 4294967, &cbdata, &launch_count);
	cerr << "e = " << cbdata.expon << "\n";
	cerr << ( (power(g,cbdata.expon) == y) ? "CORRECT!" : "INCORRECT!" ) << "\n";
    }

    return 0;
}
#endif