cudadl-0.9

Makefile

@@ -1,4 +1,4 @@
-#  cudadl version 0.8: Compute discrete logs in smooth group orders
+#  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

README

@@ -1,5 +1,5 @@
-cudadl-0.8
-23 Jan 2012
+cudadl-0.9
+21 Mar 2012
 Ryan Henry and Ian Goldberg
 {rhenry,iang}@cs.uwaterloo.ca
 http://crysp.uwaterloo.ca/software/
@@ -44,3 +44,13 @@ http://www.cacr.math.uwaterloo.ca/techreports/2012/cacr2012-02.pdf
 
 This program is covered under version 3 of the GNU General Public
 Licence; see the file COPYING for more information.
+
+Changelog:
+
+0.9 (21 Mar 2012)
+    Extend the code to handle smoothness levels (B) larger than 2^60.  Now
+    we can handle up to 2^92.  We have successfully run a test with
+    B = 2^80.
+
+0.8 (23 Jan 2012)
+    Initial public release

cudadl.h

@@ -1,5 +1,5 @@
 /*
- *  cudadl version 0.8: Compute discrete logs in smooth group orders
+ *  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

dlrho.cc

@@ -1,5 +1,5 @@
 /*
- *  cudadl version 0.8: Compute discrete logs in smooth group orders
+ *  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

dpstream.cu

@@ -1,5 +1,5 @@
 /*
- *  cudadl version 0.8: Compute discrete logs in smooth group orders
+ *  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
@@ -30,7 +30,7 @@ __constant__ static unsigned int DPbuffersize = 0;  // in words
 static unsigned int *DPbuffer_device = NULL;
 __device__ volatile static unsigned int DPbuffertail = 0;  // in words
 
-static const unsigned int DPrecordsize = (1+WORDS+4); // in words
+static const unsigned int DPrecordsize = (1+WORDS+6); // in words
 static unsigned int DPbufsize = 0;
 
 cudaError_t DPstreamInit(size_t numrecords)
@@ -61,14 +61,16 @@ void DPstreamEnd(void)
     DPbuffer_device = NULL;
 }
 
-__device__ inline void write_dp(unsigned int*, unsigned long, unsigned long);
+__device__ inline void write_dp(unsigned int*, unsigned int, unsigned int,
+    unsigned int, unsigned int, unsigned int, unsigned int);
 
 // Call this from the device to write the following into the buffer:
 //   block and thread indices (1 word total)
 //   x (WORDS words)
 //   a (2 words)
 //   b (2 words)
-__device__ inline void DPstreamWrite(unsigned long a, unsigned long b)
+__device__ inline void DPstreamWrite(unsigned int a_0, unsigned int a_1,
+    unsigned int a_2, unsigned int b_0, unsigned int b_1, unsigned int b_2)
 {
     if (!DPbuffer) return;
 
@@ -82,11 +84,12 @@ __device__ inline void DPstreamWrite(unsigned long a, unsigned long b)
     }
     unsigned int *ourbuffer = DPbuffer + ouroffset;
 
-    write_dp(ourbuffer, a, b);
+    write_dp(ourbuffer, a_0, a_1, a_2, b_0, b_1, b_2);
 }
 
 static void dpcallback(void *data, unsigned short threadId,
-    unsigned short blockId, string x, unsigned long a, unsigned long b);
+    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);
 
 void DPstreamParse(void *data)
 {
@@ -102,15 +105,15 @@ void DPstreamParse(void *data)
 	for (unsigned int *dp=dpbuf; dp < dpbuf+bufsize; dp += DPrecordsize) {
 	    unsigned short threadId = (unsigned short)(dp[0]);
 	    unsigned short blockId = dp[0]>>16;
-	    unsigned long a = dp[WORDS+2];
-	    a <<= 32;
-	    a += dp[WORDS+1];
-	    unsigned long b = dp[WORDS+4];
-	    b <<= 32;
-	    b += dp[WORDS+3];
+	    unsigned int a_0 = dp[WORDS+1];
+	    unsigned int a_1 = dp[WORDS+2];
+	    unsigned int a_2 = dp[WORDS+3];
+	    unsigned int b_0 = dp[WORDS+4];
+	    unsigned int b_1 = dp[WORDS+5];
+	    unsigned int b_2 = dp[WORDS+6];
 	    string x((const char *)(dp+1), WORDS*sizeof(unsigned int));
 
-	    dpcallback(data, threadId, blockId, x, a, b);
+	    dpcallback(data, threadId, blockId, x, a_0, a_1, a_2, b_0, b_1, b_2);
 	}
     }
     free(dpbuf);

gen_N.cc

@@ -1,5 +1,5 @@
 /*
- *  cudadl version 0.8: Compute discrete logs in smooth group orders
+ *  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

gencios_reg_20

@@ -1,6 +1,6 @@
 #!/usr/bin/perl
 
-#  cudadl version 0.8: Compute discrete logs in smooth group orders
+#  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
@@ -302,7 +302,9 @@ print "}\n";
 
 print <<'EOA';
 
-__device__ inline void write_dp(unsigned int *buf, unsigned long a, unsigned long b)
+__device__ inline void write_dp(unsigned int *buf, unsigned int a_0,
+    unsigned int a_1, unsigned int a_2, unsigned int b_0, unsigned int b_1,
+    unsigned int b_2)
 {
     unsigned short blockid = blockIdx.x + gridDim.x*blockIdx.y;
     unsigned short threadid = threadIdx.x + blockDim.x*threadIdx.y + blockDim.x*blockDim.y*threadIdx.z;
@@ -315,9 +317,11 @@ for($i=0;$i<$WORDS;++$i) {
 }
 print <<'EOA';
     : : "l" (buf));
-    buf[WORDS+1] = (unsigned int)a;
-    buf[WORDS+2] = a>>32;
-    buf[WORDS+3] = (unsigned int)b;
-    buf[WORDS+4] = b>>32;
+    buf[WORDS+1] = a_0;
+    buf[WORDS+2] = a_1;
+    buf[WORDS+3] = a_2;
+    buf[WORDS+4] = b_0;
+    buf[WORDS+5] = b_1;
+    buf[WORDS+6] = b_2;
 }
 EOA

parrhoasm.cu

@@ -1,5 +1,5 @@
 /*
- *  cudadl version 0.8: Compute discrete logs in smooth group orders
+ *  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
@@ -81,7 +81,7 @@ static void checkCUDAError(const char *msg)
 typedef struct
 {
     unsigned int x[SIZE_X];
-    unsigned int a[2], b[2];
+    unsigned int a[3], b[3];
 } GlobalThreadState;
 
 
@@ -187,14 +187,16 @@ __device__ inline bool _gt(const unsigned int * x)
 
 #define threadDimx 32
 
-__global__ void cudaMulmod(GlobalThreadState *global_ts, unsigned long long order)
+__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 long long a, b;
+    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;
@@ -204,12 +206,12 @@ __global__ void cudaMulmod(GlobalThreadState *global_ts, unsigned long long orde
 
     // cuPrintf("d_z = %08X%08X%08X\n", global_x_base[2], global_x_base[1], global_x_base[0]);
     loadx(global_x_base);
-    a = global_a_base[1];
-    a <<= 32;
-    a |= global_a_base[0];
-    b = global_b_base[1];
-    b <<= 32;
-    b |= global_b_base[0];
+    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]);
@@ -226,24 +228,45 @@ __global__ void cudaMulmod(GlobalThreadState *global_ts, unsigned long long orde
         if (xlow < TWO_32_DIV_3)
         {
             multtype = 0;
-	    a += 1;
+	    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;
-	    b += 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;
-	    a += a;
-	    b += b;
+	    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 > order) a -= order;
-	if (b > order) b -= order;
+	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));
 	/*
@@ -257,21 +280,23 @@ __global__ void cudaMulmod(GlobalThreadState *global_ts, unsigned long long orde
 
 	// Check for a distinguished point
 	asm("mov.u32 %0, $xr0;" : "=r" (xlow));
-	if ((xlow & 0x000003ff) == 0 || order < (1<<20)) {
-	    DPstreamWrite(a,b);
+	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[1] = (a>>32);
-    global_a_base[0] = (unsigned int)a;
-    global_b_base[1] = (b>>32);
-    global_b_base[0] = (unsigned int)b;
+    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) evilhack();
+    if (!(order_0|order_1|order_2)) evilhack();
 }
 
 int nthreads = 25600;
@@ -294,13 +319,22 @@ struct CBData {
 };
 
 static void dpcallback(void *cbdata, unsigned short threadId,
-    unsigned short blockId, string x, unsigned long a, unsigned long b)
+    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((unsigned long)a);
-    ZZ zz_b = to_ZZ((unsigned long)b);
+    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);
 
@@ -310,8 +344,8 @@ static void dpcallback(void *cbdata, unsigned short threadId,
     res = d->dtable.insert(DTable::value_type(x, ab));
     if (!res.second) {
 	// Collision!
-	ZZ adiff = res.first->second.first - a;
-	ZZ bdiff = b - res.first->second.second;
+	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) {
@@ -339,13 +373,14 @@ ZZ cuda_dl(const ZZ_p &base, const ZZ_p &target, const ZZ &order,
     }
 
     long orderbits = NumBits(order);
-    if (orderbits > 60) {
-	cerr << "order is larger than 60 bits.\n";
+    if (orderbits > 92) {
+	cerr << "order is larger than 92 bits.\n";
 	exit(1);
     }
 
-    unsigned long long orderll = trunc_long(order, 32) +
-	((unsigned long long)(trunc_long(order >> 32, 32)) << 32);
+	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";
@@ -442,8 +477,8 @@ ZZ cuda_dl(const ZZ_p &base, const ZZ_p &target, const ZZ &order,
 	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, 2 * sizeof(unsigned int));
-	BytesFromZZ((unsigned char *) l_ts[t].b, b, 2 * sizeof(unsigned int));
+	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";
@@ -486,7 +521,7 @@ ZZ cuda_dl(const ZZ_p &base, const ZZ_p &target, const ZZ &order,
 #endif
     {
 	cerr << getpid() << " Launch " << ++launchcount << "...\n";
-	cudaMulmod<<< nblocks, nthreads/nblocks /*tpb*/, 0 >>>(d_ts, orderll);
+	cudaMulmod<<< nblocks, nthreads/nblocks /*tpb*/, 0 >>>(d_ts, order_0, order_1, order_2);
 
 	cudaThreadSynchronize();
 	checkCUDAError("kernel launch");