From f936269eeeb3155076ffd47920a1f16989670b32 Mon Sep 17 00:00:00 2001
From: NIIBE Yutaka <gniibe@fsij.org>
Date: Mon, 10 Oct 2011 22:32:49 +0900
Subject: [PATCH 01/12] bignum 256-bit and 512-bit

---
 src/bn.c | 246 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 src/bn.h |  18 ++++
 2 files changed, 264 insertions(+)
 create mode 100644 src/bn.c
 create mode 100644 src/bn.h

diff --git a/src/bn.c b/src/bn.c
new file mode 100644
index 0000000..e0e7d95
--- /dev/null
+++ b/src/bn.c
@@ -0,0 +1,246 @@
+/*
+ * bn.c -- 256-bit and 512-bit bignum calculation
+ *
+ * Copyright (C) 2011 Free Software Initiative of Japan
+ * Author: NIIBE Yutaka <gniibe@fsij.org>
+ *
+ * This file is a part of Gnuk, a GnuPG USB Token implementation.
+ *
+ * Gnuk is free software: you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Gnuk 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 <stdint.h>
+#include "bn.h"
+
+int
+bn256_add (bn256 *X, const bn256 *A, const bn256 *B)
+{
+  int i;
+  uint32_t carry = 0;
+  uint32_t *px;
+  const uint32_t *pa, *pb;
+
+  px = X->words;
+  pa = A->words;
+  pb = B->words;
+
+  for (i = 0; i < BN256_WORDS; i++)
+    {
+      *px = *pa + carry;
+      carry = (*px < carry);
+      *px += *pb;
+      carry += (*px < *pb);
+      px++;
+      pa++;
+      pb++;
+    }
+
+  return carry;
+}
+
+int
+bn256_sub (bn256 *X, const bn256 *A, const bn256 *B)
+{
+  int i;
+  uint32_t borrow = 0;
+  uint32_t *px;
+  const uint32_t *pa, *pb;
+
+  px = X->words;
+  pa = A->words;
+  pb = B->words;
+
+  for (i = 0; i < BN256_WORDS; i++)
+    {
+      uint32_t borrow0 = (*pa < borrow);
+
+      *px = *pa - borrow;
+      borrow = (*px < *pb) + borrow0;
+      *px -= *pb;
+      px++;
+      pa++;
+      pb++;
+    }
+
+  return borrow;
+}
+
+void
+bn256_mul (bn512 *X, const bn256 *A, const bn256 *B)
+{
+  int i, j, k;
+  int i_beg, i_end;
+  uint32_t r0, r1, r2;
+
+  r0 = r1 = r2 = 0;
+  for (k = 0; k <= (BN256_WORDS - 1)*2; k++)
+    {
+      if (k < BN256_WORDS)
+	{
+	  i_beg = 0;
+	  i_end = k;
+	}
+      else
+	{
+	  i_beg = k - BN256_WORDS + 1;
+	  i_end = BN256_WORDS - 1;
+	}
+
+      for (i = i_beg; i <= i_end; i++)
+	{
+	  uint64_t uv;
+	  uint32_t u, v;
+	  uint32_t carry;
+
+	  j = k - i;
+
+	  uv = ((uint64_t )A->words[i])*((uint64_t )B->words[j]);
+	  v = uv;
+	  u = (uv >> 32);
+	  r0 += v;
+	  carry = (r0 < v);
+	  r1 += carry;
+	  carry = (r1 < carry);
+	  r1 += u;
+	  carry += (r1 < u);
+	  r2 += carry;
+	}
+
+      X->words[k] = r0;
+      r0 = r1;
+      r1 = r2;
+      r2 = 0;
+    }
+
+  X->words[k] = r0;
+}
+
+void
+bn256_sqr (bn512 *X, const bn256 *A)
+{
+  int i, j, k;
+  int i_beg, i_end;
+  uint32_t r0, r1, r2;
+
+  r0 = r1 = r2 = 0;
+  for (k = 0; k <= (BN256_WORDS - 1)*2; k++)
+    {
+      if (k < BN256_WORDS)
+	{
+	  i_beg = 0;
+	  i_end = k/2;
+	}
+      else
+	{
+	  i_beg = k - BN256_WORDS + 1;
+	  i_end = k/2;
+	}
+
+      for (i = i_beg; i <= i_end; i++)
+	{
+	  uint64_t uv;
+	  uint32_t u, v;
+	  uint32_t carry;
+
+	  j = k - i;
+
+	  uv = ((uint64_t )A->words[i])*((uint64_t )A->words[j]);
+	  if (i < j)
+	    {
+	      if ((uv >> 63) != 0)
+		r2++;
+	      uv <<= 1;
+	    }
+	  v = uv;
+	  u = (uv >> 32);
+	  r0 += v;
+	  carry = (r0 < v);
+	  r1 += carry;
+	  carry = (r1 < carry);
+	  r1 += u;
+	  carry += (r1 < u);
+	  r2 += carry;
+	}
+
+      X->words[k] = r0;
+      r0 = r1;
+      r1 = r2;
+      r2 = 0;
+    }
+
+  X->words[k] = r0;
+}
+
+int
+bn256_shift (bn256 *X, const bn256 *A, int shift)
+{
+  int i;
+  uint32_t carry = 0, next_carry;
+
+  if (shift > 0)
+    {
+      for (i = 0; i < BN256_WORDS; i++)
+	{
+	  next_carry = A->words[i] >> (32 - shift);
+	  X->words[i] = (A->words[i] << shift) | carry;
+	  carry = next_carry;
+	}
+    }
+  else
+    {
+      shift = -shift;
+
+      for (i = BN256_WORDS - 1; i >= 0; i--)
+	{
+	  next_carry = A->words[i] & ((1 << shift) - 1);
+	  X->words[i] = (A->words[i] >> shift) | (carry << (32 - shift));
+	  carry = next_carry;
+	}
+    }
+
+  return carry;
+}
+
+int
+bn256_is_zero (const bn256 *X)
+{
+  int i;
+
+  for (i = 0; i < BN256_WORDS; i++)
+    if (X->words[i] != 0)
+      return 0;
+
+  return 1;
+}
+
+int
+bn256_is_even (const bn256 *X)
+{
+  return !(X->words[0] & 1);
+}
+
+int
+bn256_is_ge (const bn256 *A, const bn256 *B)
+{
+  int i;
+
+  for (i = BN256_WORDS - 1; i >= 0; i--)
+    if (A->words[i] > B->words[i])
+      return 1;
+    else if (A->words[i] < B->words[i])
+      return 0;
+
+  return 1;
+}
diff --git a/src/bn.h b/src/bn.h
new file mode 100644
index 0000000..ac0adad
--- /dev/null
+++ b/src/bn.h
@@ -0,0 +1,18 @@
+#define BN256_WORDS 8
+typedef struct bn256 {
+  uint32_t words[ BN256_WORDS ]; /* Little endian */
+} bn256;
+
+#define BN512_WORDS 16
+typedef struct bn512 {
+  uint32_t words[ BN512_WORDS ]; /* Little endian */
+} bn512;
+
+int bn256_add (bn256 *X, const bn256 *A, const bn256 *B);
+int bn256_sub (bn256 *X, const bn256 *A, const bn256 *B);
+void bn256_mul (bn512 *X, const bn256 *A, const bn256 *B);
+void bn256_sqr (bn512 *X, const bn256 *A);
+int bn256_shift (bn256 *X, const bn256 *A, int shift);
+int bn256_is_zero (const bn256 *X);
+int bn256_is_even (const bn256 *X);
+int bn256_is_ge (const bn256 *A, const bn256 *B);

From 41272c2ab795ca4a3d734a0511388f3bdd133872 Mon Sep 17 00:00:00 2001
From: NIIBE Yutaka <gniibe@fsij.org>
Date: Mon, 10 Oct 2011 22:36:09 +0900
Subject: [PATCH 02/12] modulo P256 arithmetic

---
 src/modp256.c | 268 ++++++++++++++++++++++++++++++++++++++++++++++++++
 src/modp256.h |   9 ++
 2 files changed, 277 insertions(+)
 create mode 100644 src/modp256.c
 create mode 100644 src/modp256.h

diff --git a/src/modp256.c b/src/modp256.c
new file mode 100644
index 0000000..7ef0af6
--- /dev/null
+++ b/src/modp256.c
@@ -0,0 +1,268 @@
+/*
+ * modp256.c -- modulo P256 arithmetic
+ *
+ * Copyright (C) 2011 Free Software Initiative of Japan
+ * Author: NIIBE Yutaka <gniibe@fsij.org>
+ *
+ * This file is a part of Gnuk, a GnuPG USB Token implementation.
+ *
+ * Gnuk is free software: you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Gnuk 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/>.
+ *
+ */
+
+/*
+ * p256 =  2^256 - 2^224 + 2^192 + 2^96 - 1
+ */
+#include <stdint.h>
+#include <string.h>
+
+#include "bn.h"
+#include "modp256.h"
+
+/*
+256      224      192      160      128       96       64       32        0
+2^256
+  1 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000
+2^256 - 2^224
+  0 ffffffff 00000000 00000000 00000000 00000000 00000000 00000000 00000000
+2^256 - 2^224 + 2^192
+  0 ffffffff 00000001 00000000 00000000 00000000 00000000 00000000 00000000
+2^256 - 2^224 + 2^192 + 2^96
+  0 ffffffff 00000001 00000000 00000000 00000001 00000000 00000000 00000000
+2^256 - 2^224 + 2^192 + 2^96 - 1
+  0 ffffffff 00000001 00000000 00000000 00000000 ffffffff ffffffff ffffffff
+*/
+bn256 p256 = { {0xffffffff, 0xffffffff, 0xffffffff, 0x00000000,
+		0x00000000, 0x00000000, 0x00000001, 0xffffffff} };
+
+/**
+ * @brief  X = (A + B) mod p256
+ */
+void
+modp256_add (bn256 *X, const bn256 *A, const bn256 *B)
+{
+  int carry;
+
+  carry = bn256_add (X, A, B);
+  if (carry)
+    bn256_sub (X, X, P256);
+}
+
+/**
+ * @brief  X = (A - B) mod p256
+ */
+void
+modp256_sub (bn256 *X, bn256 *A, bn256 *B)
+{
+  int borrow;
+
+  borrow = bn256_sub (X, A, B);
+  if (borrow)
+    bn256_add (X, X, P256);
+}
+
+/**
+ * @brief  X = A mod p256
+ */
+void
+modp256_reduce (bn256 *X, bn512 *A)
+{
+  bn256 tmp[1];
+
+#define S1 X
+#define S2 tmp
+#define S3 tmp
+#define S4 tmp
+#define S5 tmp
+#define S6 tmp
+#define S7 tmp
+#define S8 tmp
+#define S9 tmp
+
+  S1->words[7] = A->words[7];
+  S1->words[6] = A->words[6];
+  S1->words[5] = A->words[5];
+  S1->words[4] = A->words[4];
+  S1->words[3] = A->words[3];
+  S1->words[2] = A->words[2];
+  S1->words[1] = A->words[1];
+  S1->words[0] = A->words[0];
+  /* X = S1 */
+
+  S2->words[7] = A->words[15];
+  S2->words[6] = A->words[14];
+  S2->words[5] = A->words[13];
+  S2->words[4] = A->words[12];
+  S2->words[3] = A->words[11];
+  S2->words[2] = S2->words[1] = S2->words[0] = 0;
+  /* X += 2 * S2 */
+  modp256_add (X, X, S2);
+  modp256_add (X, X, S2);
+
+  S3->words[7] = 0;
+  S3->words[6] = A->words[15];
+  S3->words[5] = A->words[14];
+  S3->words[4] = A->words[13];
+  S3->words[3] = A->words[12];
+  S3->words[2] = S3->words[1] = S3->words[0] = 0;
+  /* X += 2 * S3 */
+  modp256_add (X, X, S3);
+  modp256_add (X, X, S3);
+
+  S4->words[7] = A->words[15];
+  S4->words[6] = A->words[14];
+  S4->words[5] = S4->words[4] = S4->words[3] = 0;
+  S4->words[2] = A->words[10];
+  S4->words[1] = A->words[9];
+  S4->words[0] = A->words[8];
+  /* X += S4 */
+  modp256_add (X, X, S4);
+
+  S5->words[7] = A->words[8];
+  S5->words[6] = A->words[13];
+  S5->words[5] = A->words[15];
+  S5->words[4] = A->words[14];
+  S5->words[3] = A->words[13];
+  S5->words[2] = A->words[11];
+  S5->words[1] = A->words[10];
+  S5->words[0] = A->words[9];
+  /* X += S5 */
+  modp256_add (X, X, S5);
+
+  S6->words[7] = A->words[10];
+  S6->words[6] = A->words[8];
+  S6->words[5] = S6->words[4] = S6->words[3] = 0;
+  S6->words[2] = A->words[13];
+  S6->words[1] = A->words[12];
+  S6->words[0] = A->words[11];
+  /* X -= S6 */
+  modp256_sub (X, X, S6);
+
+  S7->words[7] = A->words[11];
+  S7->words[6] = A->words[9];
+  S7->words[5] = S7->words[4] = 0;
+  S7->words[3] = A->words[15];
+  S7->words[2] = A->words[14];
+  S7->words[1] = A->words[13];
+  S7->words[0] = A->words[12];
+  /* X -= S7 */
+  modp256_sub (X, X, S7);
+
+  S8->words[7] = A->words[12];
+  S8->words[6] = 0;
+  S8->words[5] = A->words[10];
+  S8->words[4] = A->words[9];
+  S8->words[3] = A->words[8];
+  S8->words[2] = A->words[15];
+  S8->words[1] = A->words[14];
+  S8->words[0] = A->words[13];
+  /* X -= S8 */
+  modp256_sub (X, X, S8);
+
+  S9->words[7] = A->words[13];
+  S9->words[6] = 0;
+  S9->words[5] = A->words[11];
+  S9->words[4] = A->words[10];
+  S9->words[3] = A->words[9];
+  S9->words[2] = 0;
+  S9->words[1] = A->words[15];
+  S9->words[0] = A->words[14];
+  /* X -= S9 */
+  modp256_sub (X, X, S9);
+
+  if (bn256_is_ge (X, P256))
+    bn256_sub (X, X, P256);
+}
+
+/**
+ * @brief  X = (A * B) mod p256
+ */
+void
+modp256_mul (bn256 *X, bn256 *A, bn256 *B)
+{
+  bn512 AB[1];
+
+  bn256_mul (AB, A, B);
+  modp256_reduce (X, AB);
+}
+
+/**
+ * @brief  X = A * A mod p256
+ */
+void
+modp256_sqr (bn256 *X, bn256 *A)
+{
+  bn512 AA[1];
+
+  bn256_sqr (AA, A);
+  modp256_reduce (X, AA);
+}
+
+/**
+ * @brief  C = (1 / a)  mod p256
+ */
+void
+modp256_inv (bn256 *C, const bn256 *a)
+{
+  bn256 u[1], v[1];
+  bn256 A[1] = { { { 1, 0, 0, 0, 0, 0, 0, 0 } } };
+
+  memset (C, 0, sizeof (bn256));
+  memcpy (u, a, sizeof (bn256));
+  memcpy (v, P256, sizeof (bn256));
+
+  while (!bn256_is_zero (u))
+    {
+      while (bn256_is_even (u))
+	{
+	  bn256_shift (u, u, -1);
+	  if (bn256_is_even (A))
+	    bn256_shift (A, A, -1);
+	  else
+	    {
+	      int carry = bn256_add (A, A, P256);
+
+	      bn256_shift (A, A, -1);
+	      if (carry)
+		A->words[7] |= 0x80000000;
+	    }
+	}
+
+      while (bn256_is_even (v))
+	{
+	  bn256_shift (v, v, -1);
+	  if (bn256_is_even (C))
+	    bn256_shift (C, C, -1);
+	  else
+	    {
+	      int carry = bn256_add (C, C, P256);
+
+	      bn256_shift (C, C, -1);
+	      if (carry)
+		C->words[7] |= 0x80000000;
+	    }
+	}
+
+      if (bn256_is_ge (u, v))
+	{
+	  bn256_sub (u, u, v);
+	  modp256_sub (A, A, C);
+	}
+      else
+	{
+	  bn256_sub (v, v, u);
+	  modp256_sub (C, C, A);
+	}
+    }
+}
diff --git a/src/modp256.h b/src/modp256.h
new file mode 100644
index 0000000..ed74b69
--- /dev/null
+++ b/src/modp256.h
@@ -0,0 +1,9 @@
+extern bn256 p256;
+#define P256 (&p256)
+
+void modp256_add (bn256 *X, const bn256 *A, const bn256 *B);
+void modp256_sub (bn256 *X, bn256 *A, bn256 *B);
+void modp256_reduce (bn256 *X, bn512 *A);
+void modp256_mul (bn256 *X, bn256 *A, bn256 *B);
+void modp256_sqr (bn256 *X, bn256 *A);
+void modp256_inv (bn256 *C, const bn256 *a);

From 89bc2ed0a3cc32920705bbe6f7c37588fe291844 Mon Sep 17 00:00:00 2001
From: NIIBE Yutaka <gniibe@fsij.org>
Date: Tue, 11 Oct 2011 09:39:21 +0900
Subject: [PATCH 03/12] fix API for const qualifier and add modp256_shift

---
 src/modp256.c | 33 +++++++++++++++++++++++++++++----
 src/modp256.h |  9 +++++----
 2 files changed, 34 insertions(+), 8 deletions(-)

diff --git a/src/modp256.c b/src/modp256.c
index 7ef0af6..be392af 100644
--- a/src/modp256.c
+++ b/src/modp256.c
@@ -63,7 +63,7 @@ modp256_add (bn256 *X, const bn256 *A, const bn256 *B)
  * @brief  X = (A - B) mod p256
  */
 void
-modp256_sub (bn256 *X, bn256 *A, bn256 *B)
+modp256_sub (bn256 *X, const bn256 *A, const bn256 *B)
 {
   int borrow;
 
@@ -76,7 +76,7 @@ modp256_sub (bn256 *X, bn256 *A, bn256 *B)
  * @brief  X = A mod p256
  */
 void
-modp256_reduce (bn256 *X, bn512 *A)
+modp256_reduce (bn256 *X, const bn512 *A)
 {
   bn256 tmp[1];
 
@@ -189,7 +189,7 @@ modp256_reduce (bn256 *X, bn512 *A)
  * @brief  X = (A * B) mod p256
  */
 void
-modp256_mul (bn256 *X, bn256 *A, bn256 *B)
+modp256_mul (bn256 *X, const bn256 *A, const bn256 *B)
 {
   bn512 AB[1];
 
@@ -201,7 +201,7 @@ modp256_mul (bn256 *X, bn256 *A, bn256 *B)
  * @brief  X = A * A mod p256
  */
 void
-modp256_sqr (bn256 *X, bn256 *A)
+modp256_sqr (bn256 *X, const bn256 *A)
 {
   bn512 AA[1];
 
@@ -266,3 +266,28 @@ modp256_inv (bn256 *C, const bn256 *a)
 	}
     }
 }
+
+/**
+ * @brief  X = (A << shift) mod p256
+ * @note   shift <= 32
+ */
+void
+modp256_shift (bn256 *X, const bn256 *A, int shift)
+{
+  int carry;
+
+  carry = bn256_shift (X, A, shift);
+  if (shift < 0)
+    return;
+
+  while (carry)
+    {
+      int borrow;
+
+      borrow = bn256_sub (X, X, P256);
+      carry -= borrow;
+    }
+
+  if (bn256_is_ge (X, P256))
+    bn256_sub (X, X, P256);
+}
diff --git a/src/modp256.h b/src/modp256.h
index ed74b69..3c63d80 100644
--- a/src/modp256.h
+++ b/src/modp256.h
@@ -2,8 +2,9 @@ extern bn256 p256;
 #define P256 (&p256)
 
 void modp256_add (bn256 *X, const bn256 *A, const bn256 *B);
-void modp256_sub (bn256 *X, bn256 *A, bn256 *B);
-void modp256_reduce (bn256 *X, bn512 *A);
-void modp256_mul (bn256 *X, bn256 *A, bn256 *B);
-void modp256_sqr (bn256 *X, bn256 *A);
+void modp256_sub (bn256 *X, const bn256 *A, const bn256 *B);
+void modp256_reduce (bn256 *X, const bn512 *A);
+void modp256_mul (bn256 *X, const bn256 *A, const bn256 *B);
+void modp256_sqr (bn256 *X, const bn256 *A);
 void modp256_inv (bn256 *C, const bn256 *a);
+void modp256_shift (bn256 *X, const bn256 *A, int shift);

From 2f67ef0142a1a23e90c0eab4175be935adc68a19 Mon Sep 17 00:00:00 2001
From: NIIBE Yutaka <gniibe@fsij.org>
Date: Tue, 11 Oct 2011 17:12:30 +0900
Subject: [PATCH 04/12] add arithmetic on Jacobian projective coordinates and
 computation of kG

---
 src/ec_p256.c | 247 ++++++++++++++++++++++++++++++++++++++++++++++++++
 src/jpc-ac.h  |  22 +++++
 src/jpc.c     | 135 +++++++++++++++++++++++++++
 3 files changed, 404 insertions(+)
 create mode 100644 src/ec_p256.c
 create mode 100644 src/jpc-ac.h
 create mode 100644 src/jpc.c

diff --git a/src/ec_p256.c b/src/ec_p256.c
new file mode 100644
index 0000000..2b5ea45
--- /dev/null
+++ b/src/ec_p256.c
@@ -0,0 +1,247 @@
+/*
+ * ec_p256.c - Elliptic curve over GF(p256)
+ */
+
+#include <stdint.h>
+#include <string.h>
+#include "bn.h"
+#include "modp256.h"
+#include "jpc-ac.h"
+
+#if 0
+/*
+ * Generator of Elliptic curve over GF(p256)
+ */
+bn256 Gx[1] = {
+  {{  0xd898c296, 0xf4a13945, 0x2deb33a0, 0x77037d81,
+      0x63a440f2, 0xf8bce6e5, 0xe12c4247, 0x6b17d1f2  }}
+};
+
+bn256 Gy[1] = {
+  {{  0x37bf51f5, 0xcbb64068, 0x6b315ece, 0x2bce3357,
+      0x7c0f9e16, 0x8ee7eb4a, 0xfe1a7f9b, 0x4fe342e2  }}
+};
+#endif
+
+/*
+ * w = 4
+ * m = 256
+ * d = 64
+ * e = 32
+ */
+
+static const ac precomputed_KG[15] = {
+  {
+    {{{ 0xd898c296, 0xf4a13945, 0x2deb33a0, 0x77037d81,
+	0x63a440f2, 0xf8bce6e5, 0xe12c4247, 0x6b17d1f2 }}},
+    {{{ 0x37bf51f5, 0xcbb64068, 0x6b315ece, 0x2bce3357,
+	0x7c0f9e16, 0x8ee7eb4a, 0xfe1a7f9b, 0x4fe342e2 }}}
+  }, {
+    {{{ 0x8e14db63, 0x90e75cb4, 0xad651f7e, 0x29493baa,
+	0x326e25de, 0x8492592e, 0x2811aaa5, 0x0fa822bc }}},
+    {{{ 0x5f462ee7, 0xe4112454, 0x50fe82f5, 0x34b1a650,
+	0xb3df188b, 0x6f4ad4bc, 0xf5dba80d, 0xbff44ae8 }}}
+  }, {
+    {{{ 0x097992af, 0x93391ce2, 0x0d35f1fa, 0xe96c98fd,
+	0x95e02789, 0xb257c0de, 0x89d6726f, 0x300a4bbc }}},
+    {{{ 0xc08127a0, 0xaa54a291, 0xa9d806a5, 0x5bb1eead,
+	0xff1e3c6f, 0x7f1ddb25, 0xd09b4644, 0x72aac7e0 }}}
+  }, {
+    {{{ 0xd789bd85, 0x57c84fc9, 0xc297eac3, 0xfc35ff7d,
+	0x88c6766e, 0xfb982fd5, 0xeedb5e67, 0x447d739b }}},
+    {{{ 0x72e25b32, 0x0c7e33c9, 0xa7fae500, 0x3d349b95,
+	0x3a4aaff7, 0xe12e9d95, 0x834131ee, 0x2d4825ab }}}
+  }, {
+    {{{ 0x2a1d367f, 0x13949c93, 0x1a0a11b7, 0xef7fbd2b,
+	0xb91dfc60, 0xddc6068b, 0x8a9c72ff, 0xef951932 }}},
+    {{{ 0x7376d8a8, 0x196035a7, 0x95ca1740, 0x23183b08,
+	0x022c219c, 0xc1ee9807, 0x7dbb2c9b, 0x611e9fc3 }}}
+  }, {
+    {{{ 0x0b57f4bc, 0xcae2b192, 0xc6c9bc36, 0x2936df5e,
+	0xe11238bf, 0x7dea6482, 0x7b51f5d8, 0x55066379 }}},
+    {{{ 0x348a964c, 0x44ffe216, 0xdbdefbe1, 0x9fb3d576,
+	0x8d9d50e5, 0x0afa4001, 0x8aecb851, 0x15716484 }}}
+  }, {
+    {{{ 0xfc5cde01, 0xe48ecaff, 0x0d715f26, 0x7ccd84e7,
+	0xf43e4391, 0xa2e8f483, 0xb21141ea, 0xeb5d7745 }}},
+    {{{ 0x731a3479, 0xcac917e2, 0x2844b645, 0x85f22cfe,
+	0x58006cee, 0x0990e6a1, 0xdbecc17b, 0xeafd72eb }}}
+  }, {
+    {{{ 0x313728be, 0x6cf20ffb, 0xa3c6b94a, 0x96439591,
+	0x44315fc5, 0x2736ff83, 0xa7849276, 0xa6d39677 }}},
+    {{{ 0xc357f5f4, 0xf2bab833, 0x2284059b, 0x824a920c,
+	0x2d27ecdf, 0x66b8babd, 0x9b0b8816, 0x674f8474 }}}
+  }, {
+    {{{ 0x677c8a3e, 0x2df48c04, 0x0203a56b, 0x74e02f08,
+	0xb8c7fedb, 0x31855f7d, 0x72c9ddad, 0x4e769e76 }}},
+    {{{ 0xb824bbb0, 0xa4c36165, 0x3b9122a5, 0xfb9ae16f,
+	0x06947281, 0x1ec00572, 0xde830663, 0x42b99082 }}}
+  }, {
+    {{{ 0xdda868b9, 0x6ef95150, 0x9c0ce131, 0xd1f89e79,
+	0x08a1c478, 0x7fdc1ca0, 0x1c6ce04d, 0x78878ef6 }}},
+    {{{ 0x1fe0d976, 0x9c62b912, 0xbde08d4f, 0x6ace570e,
+	0x12309def, 0xde53142c, 0x7b72c321, 0xb6cb3f5d }}}
+  }, {
+    {{{ 0xc31a3573, 0x7f991ed2, 0xd54fb496, 0x5b82dd5b,
+	0x812ffcae, 0x595c5220, 0x716b1287, 0x0c88bc4d }}},
+    {{{ 0x5f48aca8, 0x3a57bf63, 0xdf2564f3, 0x7c8181f4,
+	0x9c04e6aa, 0x18d1b5b3, 0xf3901dc6, 0xdd5ddea3 }}}
+  }, {
+    {{{ 0x3e72ad0c, 0xe96a79fb, 0x42ba792f, 0x43a0a28c,
+	0x083e49f3, 0xefe0a423, 0x6b317466, 0x68f344af }}},
+    {{{ 0x3fb24d4a, 0xcdfe17db, 0x71f5c626, 0x668bfc22,
+	0x24d67ff3, 0x604ed93c, 0xf8540a20, 0x31b9c405 }}}
+  }, {
+    {{{ 0xa2582e7f, 0xd36b4789, 0x4ec39c28, 0xd1a1014,
+	0xedbad7a0, 0x663c62c3, 0x6f461db9, 0x4052bf4b }}},
+    {{{ 0x188d25eb, 0x235a27c3, 0x99bfcc5b, 0xe724f339,
+	0x71d70cc8, 0x862be6bd, 0x90b0fc61, 0xfecf4d51 }}}
+  }, {
+    {{{ 0xa1d4cfac, 0x74346c10, 0x8526a7a4, 0xafdf5cc0,
+	0xf62bff7a, 0x123202a8, 0xc802e41a, 0x1eddbae2 }}},
+    {{{ 0xd603f844, 0x8fa0af2d, 0x4c701917, 0x36e06b7e,
+	0x73db33a0, 0x0c45f452, 0x560ebcfc, 0x43104d86 }}}
+  }, {
+    {{{ 0x0d1d78e5, 0x9615b511, 0x25c4744b, 0x66b0de32,
+	0x6aaf363a, 0x0a4a46fb, 0x84f7a21c, 0xb48e26b4 }}},
+    {{{ 0x21a01b2d, 0x06ebb0f6, 0x8b7b0f98, 0xc004e404,
+	0xfed6f668, 0x64131bcd, 0x4d4d3dab, 0xfac01540 }}}
+  }
+};
+
+static const ac precomputed_2E_KG[15] = {
+  {
+    {{{ 0x185a5943, 0x3a5a9e22, 0x5c65dfb6, 0x1ab91936,
+	0x262c71da, 0x21656b32, 0xaf22af89, 0x7fe36b40 }}},
+    {{{ 0x699ca101, 0xd50d152c, 0x7b8af212, 0x74b3d586,
+	0x07dca6f1, 0x9f09f404, 0x25b63624, 0xe697d458 }}}
+  }, {
+    {{{ 0x7512218e, 0xa84aa939, 0x74ca0141, 0xe9a521b0,
+	0x18a2e902, 0x57880b3a, 0x12a677a6, 0x4a5b5066 }}},
+    {{{ 0x4c4f3840, 0x0beada7a, 0x19e26d9d, 0x626db154,
+	0xe1627d40, 0xc42604fb, 0xeac089f1, 0xeb13461c }}}
+  }, {
+    {{{ 0x27a43281, 0xf9faed09, 0x4103ecbc, 0x5e52c414,
+	0xa815c857, 0xc342967a, 0x1c6a220a, 0x0781b829 }}},
+    {{{ 0xeac55f80, 0x5a8343ce, 0xe54a05e3, 0x88f80eee,
+	0x12916434, 0x97b2a14f, 0xf0151593, 0x690cde8d }}}
+  }, {
+    {{{ 0xf7f82f2a, 0xaee9c75d, 0x4afdf43a, 0x9e4c3587,
+	0x37371326, 0xf5622df4, 0x6ec73617, 0x8a535f56 }}},
+    {{{ 0x223094b7, 0xc5f9a0ac, 0x4c8c7669, 0xcde53386,
+	0x085a92bf, 0x37e02819, 0x68b08bd7, 0x0455c084 }}}
+  }, {
+    {{{ 0x9477b5d9, 0x0c0a6e2c, 0x876dc444, 0xf9a4bf62,
+	0xb6cdc279, 0x5050a949, 0xb77f8276, 0x06bada7a }}},
+    {{{ 0xea48dac9, 0xc8b4aed1, 0x7ea1070f, 0xdebd8a4b,
+	0x1366eb70, 0x427d4910, 0x0e6cb18a, 0x5b476dfd }}}
+  }, {
+    {{{ 0x278c340a, 0x7c5c3e44, 0x12d66f3b, 0x4d546068,
+	0xae23c5d8, 0x29a751b1, 0x8a2ec908, 0x3e29864e }}},
+    {{{ 0x26dbb850, 0x142d2a66, 0x765bd780, 0xad1744c4,
+	0xe322d1ed, 0x1f150e68, 0x3dc31e7e, 0x239b90ea }}}
+  }, {
+    {{{ 0x7a53322a, 0x78c41652, 0x09776f8e, 0x305dde67,
+	0xf8862ed4, 0xdbcab759, 0x49f72ff7, 0x820f4dd9 }}},
+    {{{ 0x2b5debd4, 0x6cc544a6, 0x7b4e8cc4, 0x75be5d93,
+	0x215c14d3, 0x1b481b1b, 0x783a05ec, 0x140406ec }}}
+  }, {
+    {{{ 0xe895df07, 0x6a703f10, 0x01876bd8, 0xfd75f3fa,
+	0x0ce08ffe, 0xeb5b06e7, 0x2783dfee, 0x68f6b854 }}},
+    {{{ 0x78712655, 0x90c76f8a, 0xf310bf7f, 0xcf5293d2,
+	0xfda45028, 0xfbc8044d, 0x92e40ce6, 0xcbe1feba }}}
+  }, {
+    {{{ 0x4396e4c1, 0xe998ceea, 0x6acea274, 0xfc82ef0b,
+	0x2250e927, 0x230f729f, 0x2f420109, 0xd0b2f94d }}},
+    {{{ 0xb38d4966, 0x4305addd, 0x624c3b45, 0x10b838f8,
+	0x58954e7a, 0x7db26366, 0x8b0719e5, 0x97145982 }}}
+  }, {
+    {{{ 0x23369fc9, 0x4bd6b726, 0x53d0b876, 0x57f2929e,
+	0xf2340687, 0xc2d5cba4, 0x4a866aba, 0x96161000 }}},
+    {{{ 0x2e407a5e, 0x49997bcd, 0x92ddcb24, 0x69ab197d,
+	0x8fe5131c, 0x2cf1f243, 0xcee75e44, 0x7acb9fad }}}
+  }, {
+    {{{ 0x23d2d4c0, 0x254e8394, 0x7aea685b, 0xf57f0c91,
+	0x6f75aaea, 0xa60d880f, 0xa333bf5b, 0x24eb9acc }}},
+    {{{ 0x1cda5dea, 0xe3de4ccb, 0xc51a6b4f, 0xfeef9341,
+	0x8bac4c4d, 0x743125f8, 0xacd079cc, 0x69f891c5 }}}
+  }, {
+    {{{ 0x702476b5, 0xeee44b35, 0xe45c2258, 0x7ed031a0,
+	0xbd6f8514, 0xb422d1e7, 0x5972a107, 0xe51f547c }}},
+    {{{ 0xc9cf343d, 0xa25bcd6f, 0x097c184e, 0x8ca922ee,
+	0xa9fe9a06, 0xa62f98b3, 0x25bb1387, 0x1c309a2b }}}
+  }, {
+    {{{ 0x1967c459, 0x9295dbeb, 0x3472c98e, 0xb0014883,
+	0x08011828, 0xc5049777, 0xa2c4e503, 0x20b87b8a }}},
+    {{{ 0xe057c277, 0x3063175d, 0x8fe582dd, 0x1bd53933,
+	0x5f69a044, 0x0d11adef, 0x919776be, 0xf5c6fa49 }}}
+  }, {
+    {{{ 0x0fd59e11, 0x8c944e76, 0x102fad5f, 0x3876cba1,
+	0xd83faa56, 0xa454c3fa, 0x332010b9, 0x1ed7d1b9 }}},
+    {{{ 0x0024b889, 0xa1011a27, 0xac0cd344, 0x05e4d0dc,
+	0xeb6a2a24, 0x52b520f0, 0x3217257a, 0x3a2b03f0 }}}
+  }, {
+    {{{ 0xdf1d043d, 0xf20fc2af, 0xb58d5a62, 0xf330240d,
+	0xa0058c3b, 0xfc7d229c, 0xc78dd9f6, 0x15fee545 }}},
+    {{{ 0x5bc98cda, 0x501e8288, 0xd046ac04, 0x41ef80e5,
+	0x461210fb, 0x557d9f49, 0xb8753f81, 0x4ab5b6b2 }}}
+  }
+};
+
+void
+calculate_kG (ac *X, const bn256 *K)
+{
+  int i;
+  int q_infinite = 1;
+  jpc Q[1];
+
+  for (i = 31; i >= 0; i--)
+    {
+      int k_i, k_i_e;
+
+      if (!q_infinite)
+	jpc_double (Q, Q);
+
+      k_i = (((K->words[6] >> i) & 1) << 3)
+	| (((K->words[4] >> i) & 1) << 2)
+	| (((K->words[2] >> i) & 1) << 1)
+	| ((K->words[0] >> i) & 1);
+      k_i_e = (((K->words[7] >> i) & 1) << 3)
+	| (((K->words[5] >> i) & 1) << 2)
+	| (((K->words[3] >> i) & 1) << 1)
+	| ((K->words[1] >> i) & 1);
+
+      if (k_i)
+	{
+	  if (q_infinite)
+	    {
+	      memcpy (Q->x, (&precomputed_KG[k_i - 1])->x, sizeof (bn256));
+	      memcpy (Q->y, (&precomputed_KG[k_i - 1])->y, sizeof (bn256));
+	      Q->z->words[0] = 1;
+	      Q->z->words[1] = Q->z->words[2] = Q->z->words[3]
+		= Q->z->words[4] = Q->z->words[5] = Q->z->words[6]
+		= Q->z->words[7] = 0;
+	      q_infinite = 0;
+	    }
+	  else
+	    jpc_add_ac (Q, Q, &precomputed_KG[k_i - 1]);
+	}
+      if (k_i_e)
+	{
+	  if (q_infinite)
+	    {
+	      memcpy (Q->x, (&precomputed_KG[k_i - 1])->x, sizeof (bn256));
+	      memcpy (Q->y, (&precomputed_KG[k_i - 1])->y, sizeof (bn256));
+	      Q->z->words[0] = 1;
+	      Q->z->words[1] = Q->z->words[2] = Q->z->words[3]
+		= Q->z->words[4] = Q->z->words[5] = Q->z->words[6]
+		= Q->z->words[7] = 0;
+	      q_infinite = 0;
+	    }
+	  else
+	    jpc_add_ac (Q, Q, &precomputed_2E_KG[k_i_e - 1]);
+	}
+    }
+
+  jpc_to_ac (X, Q);
+}
diff --git a/src/jpc-ac.h b/src/jpc-ac.h
new file mode 100644
index 0000000..8aee7f1
--- /dev/null
+++ b/src/jpc-ac.h
@@ -0,0 +1,22 @@
+/**
+ * @brief	Jacobian projective coordinates
+ */
+typedef struct
+{
+  bn256 x[1];
+  bn256 y[1];
+  bn256 z[1];
+} jpc;
+
+/**
+ * @brief	Affin coordinates
+ */
+typedef struct
+{
+  bn256 x[1];
+  bn256 y[1];
+} ac;
+
+void jpc_double (jpc *X, const jpc *A);
+void jpc_add_ac (jpc *X, const jpc *A, const ac *B);
+void jpc_to_ac (ac *X, const jpc *A);
diff --git a/src/jpc.c b/src/jpc.c
new file mode 100644
index 0000000..6ccfe35
--- /dev/null
+++ b/src/jpc.c
@@ -0,0 +1,135 @@
+/*
+ * jpc.c -- arithmetic on Jacobian projective coordinates and Affin coordinates
+ *
+ * Copyright (C) 2011 Free Software Initiative of Japan
+ * Author: NIIBE Yutaka <gniibe@fsij.org>
+ *
+ * This file is a part of Gnuk, a GnuPG USB Token implementation.
+ *
+ * Gnuk is free software: you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Gnuk 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 <stdint.h>
+#include <string.h>
+#include "bn.h"
+#include "modp256.h"
+#include "jpc-ac.h"
+
+/**
+ * @brief	X = 2 * A
+ *
+ * @param X	Destination JPC
+ * @param A	JPC
+ */
+void
+jpc_double (jpc *X, const jpc *A)
+{
+  bn256 a[1], b[1], c[1], tmp0[1];
+  bn256 *d = X->x;
+
+  modp256_sqr (a, A->y);
+  memcpy (b, a, sizeof (bn256));
+  modp256_mul (a, a, A->x);
+  modp256_shift (a, a, 2);
+
+  modp256_sqr (b, b);
+  modp256_shift (b, b, 3);
+
+  modp256_sqr (tmp0, A->z);
+  modp256_sub (c, A->x, tmp0);
+  modp256_add (tmp0, tmp0, A->x);
+  modp256_mul (tmp0, tmp0, c);
+  modp256_shift (c, tmp0, 1);
+  modp256_add (c, c, tmp0);
+
+  modp256_sqr (d, c);
+  modp256_shift (tmp0, a, 1);
+  modp256_sub (d, d, tmp0);
+
+  modp256_mul (X->z, A->y, A->z);
+  modp256_shift (X->z, X->z, 1);
+
+  modp256_sub (tmp0, a, d);
+  modp256_mul (tmp0, c, tmp0);
+  modp256_sub (X->y, tmp0, b);
+}
+
+/**
+ * @brief	X = A + B
+ *
+ * @param X	Destination JPC
+ * @param A	JPC
+ * @param B	AC
+ */
+void
+jpc_add_ac (jpc *X, const jpc *A, const ac *B)
+{
+  bn256 a[1], b[1], c[1], d[1];
+#define c_sqr a
+#define c_cube b
+#define x1_c_sqr c
+#define x1_c_sqr_2 c
+#define c_cube_plus_x1_c_sqr_2 c
+#define x1_c_sqr_copy a
+#define y3_tmp c
+#define y1_c_cube a
+
+  modp256_sqr (a, A->z);
+  memcpy (b, a, sizeof (bn256));
+  modp256_mul (a, a, B->x);
+
+  modp256_mul (b, b, A->z);
+  modp256_mul (b, b, B->y);
+
+  modp256_sub (c, a, A->x);
+  modp256_sub (d, b, A->y);
+
+  modp256_mul (X->z, A->z, c);
+
+  modp256_sqr (c_sqr, c);
+  modp256_mul (c_cube, c_sqr, c);
+
+  modp256_mul (x1_c_sqr, A->x, c_sqr);
+
+  modp256_sqr (X->x, d);
+  memcpy (x1_c_sqr_copy, x1_c_sqr, sizeof (bn256));
+  modp256_shift (x1_c_sqr_2, x1_c_sqr, 1);
+  modp256_add (c_cube_plus_x1_c_sqr_2, x1_c_sqr_2, c_cube);
+  modp256_sub (X->x, X->x, c_cube_plus_x1_c_sqr_2);
+
+  modp256_sub (y3_tmp, x1_c_sqr_copy, X->x);
+  modp256_mul (y3_tmp, y3_tmp, d);
+  modp256_mul (y1_c_cube, A->y, c_cube);
+  modp256_sub (X->y, y3_tmp, y1_c_cube);
+}
+
+/**
+ * @brief	X = convert A
+ *
+ * @param X	Destination AC
+ * @param A	JPC
+ */
+void
+jpc_to_ac (ac *X, const jpc *A)
+{
+  bn256 z_inv[1], z_inv_sqr[1];
+
+  modp256_inv (z_inv, A->z);
+  modp256_sqr (z_inv_sqr, z_inv);
+  modp256_mul (z_inv, z_inv, z_inv_sqr);
+
+  modp256_mul (X->x, A->x, z_inv_sqr);
+  modp256_mul (X->y, A->y, z_inv);
+}

From 331c1a15635e70c6f84e38e8c8f80727588cda03 Mon Sep 17 00:00:00 2001
From: NIIBE Yutaka <gniibe@fsij.org>
Date: Tue, 11 Oct 2011 22:44:16 +0900
Subject: [PATCH 05/12] change return value of carry

---
 src/bn.c | 6 +++---
 src/bn.h | 6 +++---
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/bn.c b/src/bn.c
index e0e7d95..9d97d90 100644
--- a/src/bn.c
+++ b/src/bn.c
@@ -24,7 +24,7 @@
 #include <stdint.h>
 #include "bn.h"
 
-int
+uint32_t
 bn256_add (bn256 *X, const bn256 *A, const bn256 *B)
 {
   int i;
@@ -50,7 +50,7 @@ bn256_add (bn256 *X, const bn256 *A, const bn256 *B)
   return carry;
 }
 
-int
+uint32_t
 bn256_sub (bn256 *X, const bn256 *A, const bn256 *B)
 {
   int i;
@@ -183,7 +183,7 @@ bn256_sqr (bn512 *X, const bn256 *A)
   X->words[k] = r0;
 }
 
-int
+uint32_t
 bn256_shift (bn256 *X, const bn256 *A, int shift)
 {
   int i;
diff --git a/src/bn.h b/src/bn.h
index ac0adad..a4e2f63 100644
--- a/src/bn.h
+++ b/src/bn.h
@@ -8,11 +8,11 @@ typedef struct bn512 {
   uint32_t words[ BN512_WORDS ]; /* Little endian */
 } bn512;
 
-int bn256_add (bn256 *X, const bn256 *A, const bn256 *B);
-int bn256_sub (bn256 *X, const bn256 *A, const bn256 *B);
+uint32_t bn256_add (bn256 *X, const bn256 *A, const bn256 *B);
+uint32_t bn256_sub (bn256 *X, const bn256 *A, const bn256 *B);
 void bn256_mul (bn512 *X, const bn256 *A, const bn256 *B);
 void bn256_sqr (bn512 *X, const bn256 *A);
-int bn256_shift (bn256 *X, const bn256 *A, int shift);
+uint32_t bn256_shift (bn256 *X, const bn256 *A, int shift);
 int bn256_is_zero (const bn256 *X);
 int bn256_is_even (const bn256 *X);
 int bn256_is_ge (const bn256 *A, const bn256 *B);

From 434cf67470bacfdbe6187486a71a13b64c139615 Mon Sep 17 00:00:00 2001
From: NIIBE Yutaka <gniibe@fsij.org>
Date: Tue, 11 Oct 2011 22:44:50 +0900
Subject: [PATCH 06/12] fix bug of ec_p256

---
 src/ec_p256.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/ec_p256.c b/src/ec_p256.c
index 2b5ea45..9cd7b98 100644
--- a/src/ec_p256.c
+++ b/src/ec_p256.c
@@ -230,8 +230,8 @@ calculate_kG (ac *X, const bn256 *K)
 	{
 	  if (q_infinite)
 	    {
-	      memcpy (Q->x, (&precomputed_KG[k_i - 1])->x, sizeof (bn256));
-	      memcpy (Q->y, (&precomputed_KG[k_i - 1])->y, sizeof (bn256));
+	      memcpy (Q->x, (&precomputed_2E_KG[k_i_e - 1])->x, sizeof (bn256));
+	      memcpy (Q->y, (&precomputed_2E_KG[k_i_e - 1])->y, sizeof (bn256));
 	      Q->z->words[0] = 1;
 	      Q->z->words[1] = Q->z->words[2] = Q->z->words[3]
 		= Q->z->words[4] = Q->z->words[5] = Q->z->words[6]

From 23c411188cc092ae78c1c7ddbf038cf91efbdb48 Mon Sep 17 00:00:00 2001
From: NIIBE Yutaka <gniibe@fsij.org>
Date: Tue, 11 Oct 2011 22:46:39 +0900
Subject: [PATCH 07/12] fast reduction for modp256_shift

---
 src/modp256.c | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/src/modp256.c b/src/modp256.c
index be392af..69ced35 100644
--- a/src/modp256.c
+++ b/src/modp256.c
@@ -274,19 +274,23 @@ modp256_inv (bn256 *C, const bn256 *a)
 void
 modp256_shift (bn256 *X, const bn256 *A, int shift)
 {
-  int carry;
+  uint32_t carry;
+  bn256 tmp[1];
 
   carry = bn256_shift (X, A, shift);
   if (shift < 0)
     return;
 
-  while (carry)
-    {
-      int borrow;
+  memset (tmp, 0, sizeof (bn256));
+  tmp->words[7] = carry;
+  tmp->words[0] = carry;
+  modp256_add (X, X, tmp);
 
-      borrow = bn256_sub (X, X, P256);
-      carry -= borrow;
-    }
+  tmp->words[7] = 0;
+  tmp->words[0] = 0;
+  tmp->words[6] = carry;
+  tmp->words[3] = carry;
+  modp256_sub (X, X, tmp);
 
   if (bn256_is_ge (X, P256))
     bn256_sub (X, X, P256);

From d6dc9645367e1a5cc6845e6f20115917e0a1a4dd Mon Sep 17 00:00:00 2001
From: NIIBE Yutaka <gniibe@fsij.org>
Date: Wed, 12 Oct 2011 14:40:31 +0900
Subject: [PATCH 08/12] modulo arithmetic

---
 src/mod.c | 195 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 src/mod.h |   2 +
 2 files changed, 197 insertions(+)
 create mode 100644 src/mod.c
 create mode 100644 src/mod.h

diff --git a/src/mod.c b/src/mod.c
new file mode 100644
index 0000000..018f6a3
--- /dev/null
+++ b/src/mod.c
@@ -0,0 +1,195 @@
+/*
+ * mod.c -- modulo arithmetic
+ *
+ * Copyright (C) 2011 Free Software Initiative of Japan
+ * Author: NIIBE Yutaka <gniibe@fsij.org>
+ *
+ * This file is a part of Gnuk, a GnuPG USB Token implementation.
+ *
+ * Gnuk is free software: you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Gnuk 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 <stdint.h>
+#include <string.h>
+#include "bn.h"
+
+/**
+ * @brief X = A mod B (using MU=(1<<(256)+MU_lower)) (Barret reduction)
+ * 
+ */
+void
+mod_reduce (bn256 *X, const bn512 *A, const bn256 *B, const bn256 *MU_lower)
+{
+  bn256 q[1];
+  bn512 q_big[1], tmp[1];
+  uint32_t carry;
+#define borrow carry
+  uint32_t borrow_next;
+
+  memset (q, 0, sizeof (bn256));
+  q->words[0] = A->words[15];
+  bn256_mul (tmp, q, MU_lower);
+  tmp->words[8] += A->words[15];
+  carry = (tmp->words[8] < A->words[15]);
+  tmp->words[9] += carry;
+
+  q->words[7] = A->words[14];
+  q->words[6] = A->words[13];
+  q->words[5] = A->words[12];
+  q->words[4] = A->words[11];
+  q->words[3] = A->words[10];
+  q->words[2] = A->words[9];
+  q->words[1] = A->words[8];
+  q->words[0] = A->words[7];
+  bn256_mul (q_big, q, MU_lower);
+  bn256_add ((bn256 *)&q_big->words[8], (bn256 *)&q_big->words[8], q);
+
+  q->words[0] = q_big->words[9] + tmp->words[1];
+  carry = (q->words[0] < tmp->words[1]);
+
+  q->words[1] = q_big->words[10] + carry;
+  carry = (q->words[1] < carry);
+  q->words[1] += tmp->words[2];
+  carry += (q->words[1] < tmp->words[2]);
+
+  q->words[2] = q_big->words[11] + carry;
+  carry = (q->words[2] < carry);
+  q->words[2] += tmp->words[3];
+  carry += (q->words[2] < tmp->words[3]);
+
+  q->words[3] = q_big->words[12] + carry;
+  carry = (q->words[3] < carry);
+  q->words[3] += tmp->words[4];
+  carry += (q->words[3] < tmp->words[4]);
+
+  q->words[4] = q_big->words[13] + carry;
+  carry = (q->words[4] < carry);
+  q->words[4] += tmp->words[5];
+  carry += (q->words[4] < tmp->words[5]);
+
+  q->words[5] = q_big->words[14] + carry;
+  carry = (q->words[5] < carry);
+  q->words[5] += tmp->words[6];
+  carry += (q->words[5] < tmp->words[6]);
+
+  q->words[6] = q_big->words[15] + carry;
+  carry = (q->words[6] < carry);
+  q->words[6] += tmp->words[7];
+  carry += (q->words[6] < tmp->words[7]);
+
+  q->words[7] = carry;
+  q->words[7] += tmp->words[8];
+  carry = (q->words[7] < tmp->words[8]);
+
+  memset (q_big, 0, sizeof (bn512));
+  q_big->words[8] = A->words[8];
+  q_big->words[7] = A->words[7];
+  q_big->words[6] = A->words[6];
+  q_big->words[5] = A->words[5];
+  q_big->words[4] = A->words[4];
+  q_big->words[3] = A->words[3];
+  q_big->words[2] = A->words[2];
+  q_big->words[1] = A->words[1];
+  q_big->words[0] = A->words[0];
+
+  bn256_mul (tmp, q, B);
+  if (carry)
+    tmp->words[8] += B->words[0];
+  tmp->words[15] = tmp->words[14] = tmp->words[13] = tmp->words[12]
+    = tmp->words[11] = tmp->words[10] = tmp->words[9] = 0;
+
+  borrow = bn256_sub (X, (bn256 *)&q_big->words[0], (bn256 *)&tmp->words[0]);
+  borrow_next = (q_big->words[8] < borrow);
+  q_big->words[8] -= borrow;
+  borrow_next += (q_big->words[8] < tmp->words[8]);
+  q_big->words[8] -= tmp->words[8];
+
+  carry = q_big->words[8];
+  while (carry)
+    {
+      borrow_next = bn256_sub (X, X, B);
+      carry -= borrow_next;
+    }
+
+  if (bn256_is_ge (X, B))
+    bn256_sub (X, X, B);
+}
+
+/**
+ * @brief C = X^(-1) mod N
+ * 
+ */
+void
+mod_inv (bn256 *C, const bn256 *X, const bn256 *N)
+{
+  bn256 u[1], v[1];
+  bn256 A[1] = { { { 1, 0, 0, 0, 0, 0, 0, 0 } } };
+
+  memset (C, 0, sizeof (bn256));
+  memcpy (u, X, sizeof (bn256));
+  memcpy (v, N, sizeof (bn256));
+
+  while (!bn256_is_zero (u))
+    {
+      while (bn256_is_even (u))
+	{
+	  bn256_shift (u, u, -1);
+	  if (bn256_is_even (A))
+	    bn256_shift (A, A, -1);
+	  else
+	    {
+	      int carry = bn256_add (A, A, N);
+
+	      bn256_shift (A, A, -1);
+	      if (carry)
+		A->words[7] |= 0x80000000;
+	    }
+	}
+
+      while (bn256_is_even (v))
+	{
+	  bn256_shift (v, v, -1);
+	  if (bn256_is_even (C))
+	    bn256_shift (C, C, -1);
+	  else
+	    {
+	      int carry = bn256_add (C, C, N);
+
+	      bn256_shift (C, C, -1);
+	      if (carry)
+		C->words[7] |= 0x80000000;
+	    }
+	}
+
+      if (bn256_is_ge (u, v))
+	{
+	  int borrow;
+
+	  bn256_sub (u, u, v);
+	  borrow = bn256_sub (A, A, C);
+	  if (borrow)
+	    bn256_add (A, A, N);
+	}
+      else
+	{
+	  int borrow;
+
+	  bn256_sub (v, v, u);
+	  borrow = bn256_sub (C, C, A);
+	  if (borrow)
+	    bn256_add (C, C, N);
+	}
+    }
+}
diff --git a/src/mod.h b/src/mod.h
new file mode 100644
index 0000000..8f4b9c1
--- /dev/null
+++ b/src/mod.h
@@ -0,0 +1,2 @@
+void mod_reduce (bn256 *X, const bn512 *A, const bn256 *B, const bn256 *MU_lower);
+void mod_inv (bn256 *X, const bn256 *A, const bn256 *N);

From fb97562fdc13c8dc1d9cccc53a9666281081e9c6 Mon Sep 17 00:00:00 2001
From: NIIBE Yutaka <gniibe@fsij.org>
Date: Wed, 12 Oct 2011 18:07:52 +0900
Subject: [PATCH 09/12] ecdsa

---
 src/bn.h      |  1 +
 src/ec_p256.c | 57 +++++++++++++++++++++++++++++++++++++++++++++++++++
 src/ec_p256.h |  3 +++
 3 files changed, 61 insertions(+)
 create mode 100644 src/ec_p256.h

diff --git a/src/bn.h b/src/bn.h
index a4e2f63..5def1c4 100644
--- a/src/bn.h
+++ b/src/bn.h
@@ -16,3 +16,4 @@ uint32_t bn256_shift (bn256 *X, const bn256 *A, int shift);
 int bn256_is_zero (const bn256 *X);
 int bn256_is_even (const bn256 *X);
 int bn256_is_ge (const bn256 *A, const bn256 *B);
+void bn256_random (bn256 *X);
diff --git a/src/ec_p256.c b/src/ec_p256.c
index 9cd7b98..b5b2895 100644
--- a/src/ec_p256.c
+++ b/src/ec_p256.c
@@ -7,6 +7,7 @@
 #include "bn.h"
 #include "modp256.h"
 #include "jpc-ac.h"
+#include "mod.h"
 
 #if 0
 /*
@@ -245,3 +246,59 @@ calculate_kG (ac *X, const bn256 *K)
 
   jpc_to_ac (X, Q);
 }
+
+
+/*
+ * N: order of G
+ */
+static const bn256 N[1] = {
+  {{ 0xfc632551, 0xf3b9cac2, 0xa7179e84, 0xbce6faad,
+     0xffffffff, 0xffffffff, 0x00000000, 0xffffffff }}
+};
+
+/*
+ * MU = 2^512 / N
+ * MU = ( (1 << 256) | MU_lower )
+ */
+static const bn256 MU_lower[1] = {
+  {{ 0xeedf9bfe, 0x012ffd85, 0xdf1a6c21, 0x43190552,
+     0xffffffff, 0xfffffffe, 0xffffffff, 0x00000000 }}
+};
+
+
+/**
+ * @brief Compute signature (r,s) of hash string z with secret key dA
+ */
+void
+ecdsa (bn256 *r, bn256 *s, const bn256 *z, const bn256 *d)
+{
+  bn256 k[1];
+  ac KG[1];
+  bn512 tmp[1];
+  bn256 k_inv[1];
+  uint32_t carry;
+
+  do
+    {
+      do
+	{
+	  bn256_random (k);
+	  calculate_kG (KG, k);
+	  if (bn256_is_ge (KG->x, N))
+	    bn256_sub (r, KG->x, N);
+	  else
+	    memcpy (r, KG->x, sizeof (bn256));
+	}
+      while (bn256_is_zero (r));
+
+      mod_inv (k_inv, k, N);
+      bn256_mul (tmp, r, d);
+      mod_reduce (s, tmp, N, MU_lower);
+      carry = bn256_add (s, s, z);
+      if (carry)
+	bn256_sub (s, s, N);
+      bn256_mul (tmp, s, k_inv);
+      mod_reduce (s, tmp, N, MU_lower);
+    }
+  while (bn256_is_zero (s));
+}
diff --git a/src/ec_p256.h b/src/ec_p256.h
new file mode 100644
index 0000000..ac9df77
--- /dev/null
+++ b/src/ec_p256.h
@@ -0,0 +1,3 @@
+void calculate_kG (ac *X, const bn256 *K);
+void ecdsa (bn256 *r, bn256 *s, const bn256 *z, const bn256 *d);
+

From 8b6713bc523d12bd6d35a9cef59dd4f99fbb6011 Mon Sep 17 00:00:00 2001
From: NIIBE Yutaka <gniibe@fsij.org>
Date: Thu, 13 Oct 2011 15:43:58 +0900
Subject: [PATCH 10/12] compute_KP (for ECDH)

---
 src/Makefile.in |   2 +-
 src/bn.c        |  74 +++++++++++++++++++++++
 src/bn.h        |   3 +
 src/ec_p256.c   | 155 ++++++++++++++++++++++++++++++++++++++++++++++--
 src/ec_p256.h   |  10 +++-
 src/jpc-ac.h    |   1 +
 src/jpc.c       |  25 +++++++-
 7 files changed, 261 insertions(+), 9 deletions(-)

diff --git a/src/Makefile.in b/src/Makefile.in
index cc77479..c8b0827 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -91,7 +91,7 @@ CSRC = $(PORTSRC) \
        main.c usb_lld.c \
 	usb_desc.c usb_prop.c \
 	usb-icc.c openpgp.c ac.c openpgp-do.c flash.c hardclock.c \
-	random.c neug.c
+	random.c neug.c bn.c modp256.c jpc.c mod.c ec_p256.c
 
 ifneq ($(ENABLE_DEBUG),)
 CSRC += debug.c
diff --git a/src/bn.c b/src/bn.c
index 9d97d90..9987d3d 100644
--- a/src/bn.c
+++ b/src/bn.c
@@ -77,6 +77,65 @@ bn256_sub (bn256 *X, const bn256 *A, const bn256 *B)
   return borrow;
 }
 
+uint32_t
+bn256_add_uint (bn256 *X, const bn256 *A, uint32_t w)
+{
+  int i;
+  uint32_t carry = 0;
+  uint32_t *px;
+  const uint32_t *pa;
+
+  px = X->words;
+  pa = A->words;
+
+  for (i = 0; i < BN256_WORDS; i++)
+    {
+      *px = *pa + carry;
+      carry = (*px < carry);
+      if (i == 0)
+	{
+	  *px += w;
+	  carry += (*px < w);
+	}
+
+      px++;
+      pa++;
+    }
+
+  return carry;
+}
+
+uint32_t
+bn256_sub_uint (bn256 *X, const bn256 *A, uint32_t w)
+{
+  int i;
+  uint32_t borrow = 0;
+  uint32_t *px;
+  const uint32_t *pa;
+
+  px = X->words;
+  pa = A->words;
+
+  for (i = 0; i < BN256_WORDS; i++)
+    {
+      uint32_t borrow0 = (*pa < borrow);
+
+      *px = *pa - borrow;
+      if (i == 0)
+	{
+	  borrow = (*px < w) + borrow0;
+	  *px -= w;
+	}
+      else
+	borrow = borrow0;
+
+      px++;
+      pa++;
+    }
+
+  return borrow;
+}
+
 void
 bn256_mul (bn512 *X, const bn256 *A, const bn256 *B)
 {
@@ -244,3 +303,18 @@ bn256_is_ge (const bn256 *A, const bn256 *B)
 
   return 1;
 }
+
+void
+bn256_random (bn256 *X)
+{
+#if 1
+  X->words[7] = 0x01234567;
+  X->words[6] = 0x89abcdef;
+  X->words[5] = 0xff00ff00;
+  X->words[4] = 0x00ff00ff;
+  X->words[3] = 0xee55ee55;
+  X->words[2] = 0x55ee55ee;
+  X->words[1] = 0x01234567;
+  X->words[0] = 0x89abcdef;
+#endif
+}
diff --git a/src/bn.h b/src/bn.h
index 5def1c4..32d9d62 100644
--- a/src/bn.h
+++ b/src/bn.h
@@ -10,6 +10,9 @@ typedef struct bn512 {
 
 uint32_t bn256_add (bn256 *X, const bn256 *A, const bn256 *B);
 uint32_t bn256_sub (bn256 *X, const bn256 *A, const bn256 *B);
+uint32_t bn256_add_uint (bn256 *X, const bn256 *A, uint32_t w);
+uint32_t bn256_sub_uint (bn256 *X, const bn256 *A, uint32_t w);
+
 void bn256_mul (bn512 *X, const bn256 *A, const bn256 *B);
 void bn256_sqr (bn512 *X, const bn256 *A);
 uint32_t bn256_shift (bn256 *X, const bn256 *A, int shift);
diff --git a/src/ec_p256.c b/src/ec_p256.c
index b5b2895..2b8a031 100644
--- a/src/ec_p256.c
+++ b/src/ec_p256.c
@@ -8,6 +8,7 @@
 #include "modp256.h"
 #include "jpc-ac.h"
 #include "mod.h"
+#include "ec_p256.h"
 
 #if 0
 /*
@@ -190,7 +191,7 @@ static const ac precomputed_2E_KG[15] = {
 };
 
 void
-calculate_kG (ac *X, const bn256 *K)
+compute_kG (ac *X, const bn256 *K)
 {
   int i;
   int q_infinite = 1;
@@ -233,10 +234,8 @@ calculate_kG (ac *X, const bn256 *K)
 	    {
 	      memcpy (Q->x, (&precomputed_2E_KG[k_i_e - 1])->x, sizeof (bn256));
 	      memcpy (Q->y, (&precomputed_2E_KG[k_i_e - 1])->y, sizeof (bn256));
+	      memset (Q->z, 0, sizeof (bn256));
 	      Q->z->words[0] = 1;
-	      Q->z->words[1] = Q->z->words[2] = Q->z->words[3]
-		= Q->z->words[4] = Q->z->words[5] = Q->z->words[6]
-		= Q->z->words[7] = 0;
 	      q_infinite = 0;
 	    }
 	  else
@@ -248,6 +247,152 @@ calculate_kG (ac *X, const bn256 *K)
 }
 
 
+#define NAF_K_SIGN(k)	(k&8)
+#define NAF_K_INDEX(k)	((k&7)-1)
+
+void
+naf4_257_set (naf4_257 *NAF_K, int i, int ki)
+{
+  if (ki != 0)
+    {
+      if (ki > 0)
+	ki = (ki+1)/2;
+      else
+	ki = (1-ki)/2 | 8;
+    }
+
+  if (i == 256)
+    NAF_K->last_nibble = ki;
+  else
+    {
+      NAF_K->words[i/8] &= ~(0x0f << ((i & 0x07)*4));
+      NAF_K->words[i/8] |= (ki << ((i & 0x07)*4));
+    }
+}
+
+int
+naf4_257_get (const naf4_257 *NAF_K, int i)
+{
+  int ki;
+
+  if (i == 256)
+    ki = NAF_K->last_nibble;
+  else
+    {
+      ki = NAF_K->words[i/8] >> ((i & 0x07)*4);
+      ki &= 0x0f;
+    }
+
+  return ki;
+}
+
+void
+compute_naf4_257 (naf4_257 *NAF_K, const bn256 *K)
+{
+  int i = 0;
+  bn256 K_tmp[1];
+  uint32_t carry = 0;
+
+  memcpy (K_tmp, K, sizeof (bn256));
+  memset (NAF_K, 0, sizeof (naf4_257));
+
+  while (!bn256_is_zero (K_tmp))
+    {
+      if (bn256_is_even (K_tmp))
+	naf4_257_set (NAF_K, i, 0);
+      else
+	{
+	  int ki = (K_tmp->words[0]) & 0x0f;
+
+	  if ((ki & 0x08))
+	    {
+	      carry = bn256_add_uint (K_tmp, K_tmp, 16 - ki);
+	      ki = ki - 16;
+	    }
+	  else
+	    K_tmp->words[0] &= 0xfffffff0;
+
+	  naf4_257_set (NAF_K, i, ki);
+	}
+
+      bn256_shift (K_tmp, K_tmp, -1);
+      if (carry)
+	{
+	  K_tmp->words[7] |= 0x80000000;
+	  carry = 0;
+	}
+      i++;
+    }
+}
+
+void
+compute_kP (ac *X, const naf4_257 *NAF_K, const ac *P)
+{
+  int i;
+  int q_infinite = 1;
+  jpc Q[1];
+  ac P3[1], P5[1], P7[1];
+  const ac *p_Pi[4];
+
+  p_Pi[0] = P;
+  p_Pi[1] = P3;
+  p_Pi[2] = P5;
+  p_Pi[3] = P7;
+
+  {
+    jpc Q1[1];
+
+    memcpy (Q->x, P->x, sizeof (bn256));
+    memcpy (Q->y, P->y, sizeof (bn256));
+    memset (Q->z, 0, sizeof (bn256));
+    Q->z->words[0] = 1;
+
+    jpc_double (Q, Q);
+    jpc_add_ac (Q1, Q, P);
+    jpc_to_ac (P3, Q1);
+    jpc_double (Q, Q);
+    jpc_add_ac (Q1, Q, P);
+    jpc_to_ac (P5, Q1);
+
+    memcpy (Q->x, P3->x, sizeof (bn256));
+    memcpy (Q->y, P3->y, sizeof (bn256));
+    memset (Q->z, 0, sizeof (bn256));
+    Q->z->words[0] = 1;
+    jpc_double (Q, Q);
+    jpc_add_ac (Q1, Q, P);
+    jpc_to_ac (P7, Q1);
+  }
+
+  for (i = 256; i >= 0; i--)
+    {
+      int k_i;
+
+      if (!q_infinite)
+	jpc_double (Q, Q);
+
+      k_i = naf4_257_get (NAF_K, i);
+      if (k_i)
+	{
+	  if (q_infinite)
+	    {
+	      memcpy (Q->x, p_Pi[NAF_K_INDEX(k_i)]->x, sizeof (bn256));
+	      if (NAF_K_SIGN (k_i))
+		bn256_sub (Q->y, P256, p_Pi[NAF_K_INDEX(k_i)]->y);
+	      else
+		memcpy (Q->y, p_Pi[NAF_K_INDEX(k_i)]->y, sizeof (bn256));
+	      memset (Q->z, 0, sizeof (bn256));
+	      Q->z->words[0] = 1;
+	      q_infinite = 0;
+	    }
+	  else
+	    jpc_add_ac_signed (Q, Q, p_Pi[NAF_K_INDEX(k_i)], NAF_K_SIGN (k_i));
+	}
+    }
+
+  jpc_to_ac (X, Q);
+}
+
+
 /*
  * N: order of G
  */
@@ -283,7 +428,7 @@ ecdsa (bn256 *r, bn256 *s, const bn256 *z, const bn256 *d)
       do
 	{
 	  bn256_random (k);
-	  calculate_kG (KG, k);
+	  compute_kG (KG, k);
 	  if (bn256_is_ge (KG->x, N))
 	    bn256_sub (r, KG->x, N);
 	  else
diff --git a/src/ec_p256.h b/src/ec_p256.h
index ac9df77..3e4c199 100644
--- a/src/ec_p256.h
+++ b/src/ec_p256.h
@@ -1,3 +1,11 @@
-void calculate_kG (ac *X, const bn256 *K);
+typedef struct naf4_257 {
+  uint32_t words[ BN256_WORDS*4 ]; /* Little endian */
+  uint8_t last_nibble;		   /* most significant nibble */
+} naf4_257;
+
+void compute_naf4_257 (naf4_257 *NAF_K, const bn256 *K);
+void compute_kP (ac *X, const naf4_257 *NAF_K, const ac *P);
+
+void compute_kG (ac *X, const bn256 *K);
 void ecdsa (bn256 *r, bn256 *s, const bn256 *z, const bn256 *d);
 
diff --git a/src/jpc-ac.h b/src/jpc-ac.h
index 8aee7f1..6917b2d 100644
--- a/src/jpc-ac.h
+++ b/src/jpc-ac.h
@@ -19,4 +19,5 @@ typedef struct
 
 void jpc_double (jpc *X, const jpc *A);
 void jpc_add_ac (jpc *X, const jpc *A, const ac *B);
+void jpc_add_ac_signed (jpc *X, const jpc *A, const ac *B, int minus);
 void jpc_to_ac (ac *X, const jpc *A);
diff --git a/src/jpc.c b/src/jpc.c
index 6ccfe35..f2a3d13 100644
--- a/src/jpc.c
+++ b/src/jpc.c
@@ -72,11 +72,13 @@ jpc_double (jpc *X, const jpc *A)
  * @param X	Destination JPC
  * @param A	JPC
  * @param B	AC
+ * @param MINUS if 1 subtraction, addition otherwise.
  */
 void
-jpc_add_ac (jpc *X, const jpc *A, const ac *B)
+jpc_add_ac_signed (jpc *X, const jpc *A, const ac *B, int minus)
 {
   bn256 a[1], b[1], c[1], d[1];
+#define minus_B_y c
 #define c_sqr a
 #define c_cube b
 #define x1_c_sqr c
@@ -91,7 +93,13 @@ jpc_add_ac (jpc *X, const jpc *A, const ac *B)
   modp256_mul (a, a, B->x);
 
   modp256_mul (b, b, A->z);
-  modp256_mul (b, b, B->y);
+  if (minus)
+    {
+      bn256_sub (minus_B_y, P256, B->y);
+      modp256_mul (b, b, minus_B_y);
+    }
+  else
+    modp256_mul (b, b, B->y);
 
   modp256_sub (c, a, A->x);
   modp256_sub (d, b, A->y);
@@ -115,6 +123,19 @@ jpc_add_ac (jpc *X, const jpc *A, const ac *B)
   modp256_sub (X->y, y3_tmp, y1_c_cube);
 }
 
+/**
+ * @brief	X = A + B
+ *
+ * @param X	Destination JPC
+ * @param A	JPC
+ * @param B	AC
+ */
+void
+jpc_add_ac (jpc *X, const jpc *A, const ac *B)
+{
+  jpc_add_ac_signed (X, A, B, 0);
+}
+
 /**
  * @brief	X = convert A
  *

From e929963d30f3db450a3dc8320dd2bd8d0910dd39 Mon Sep 17 00:00:00 2001
From: NIIBE Yutaka <gniibe@fsij.org>
Date: Thu, 13 Oct 2011 16:54:08 +0900
Subject: [PATCH 11/12] use NeuG for bn256_random

---
 src/bn.c     | 22 ++++++++++++----------
 src/random.c |  6 +++---
 2 files changed, 15 insertions(+), 13 deletions(-)

diff --git a/src/bn.c b/src/bn.c
index 9987d3d..1e08362 100644
--- a/src/bn.c
+++ b/src/bn.c
@@ -307,14 +307,16 @@ bn256_is_ge (const bn256 *A, const bn256 *B)
 void
 bn256_random (bn256 *X)
 {
-#if 1
-  X->words[7] = 0x01234567;
-  X->words[6] = 0x89abcdef;
-  X->words[5] = 0xff00ff00;
-  X->words[4] = 0x00ff00ff;
-  X->words[3] = 0xee55ee55;
-  X->words[2] = 0x55ee55ee;
-  X->words[1] = 0x01234567;
-  X->words[0] = 0x89abcdef;
-#endif
+  const uint8_t *rand = random_bytes_get ();
+
+  X->words[7] = ((uint32_t *)rand)[7];
+  X->words[6] = ((uint32_t *)rand)[6];
+  X->words[5] = ((uint32_t *)rand)[5];
+  X->words[4] = ((uint32_t *)rand)[4];
+  X->words[3] = ((uint32_t *)rand)[3];
+  X->words[2] = ((uint32_t *)rand)[2];
+  X->words[1] = ((uint32_t *)rand)[1];
+  X->words[0] = ((uint32_t *)rand)[0];
+
+  random_bytes_free (rand);
 }
diff --git a/src/random.c b/src/random.c
index 95299f6..d099cf1 100644
--- a/src/random.c
+++ b/src/random.c
@@ -26,7 +26,7 @@
 #include "gnuk.h"
 #include "neug.h"
 
-#define RANDOM_BYTES_LENGTH 16
+#define RANDOM_BYTES_LENGTH 32
 static uint32_t random_word[RANDOM_BYTES_LENGTH/sizeof (uint32_t)];
 
 void
@@ -43,7 +43,7 @@ random_init (void)
 }
 
 /*
- * Return pointer to random 16-byte
+ * Return pointer to random 32-byte
  */
 const uint8_t *
 random_bytes_get (void)
@@ -53,7 +53,7 @@ random_bytes_get (void)
 }
 
 /*
- * Free pointer to random 16-byte
+ * Free pointer to random 32-byte
  */
 void
 random_bytes_free (const uint8_t *p)

From d70951e941edded4ab9863012a4953c5edf92e3e Mon Sep 17 00:00:00 2001
From: NIIBE Yutaka <gniibe@fsij.org>
Date: Sun, 17 Feb 2013 22:20:30 +0900
Subject: [PATCH 12/12] update and add ecc-cdh.c

---
 src/bn.c      |  2 +-
 src/ec_p256.c | 83 ++++++++++++++++++++++++++++++++++++++++++---------
 src/ec_p256.h |  1 +
 src/ecc-cdh.c | 68 +++++++++++++++++++++++++++++++++++++++++
 src/mod.h     |  3 +-
 5 files changed, 141 insertions(+), 16 deletions(-)
 create mode 100644 src/ecc-cdh.c

diff --git a/src/bn.c b/src/bn.c
index 1e08362..9127f61 100644
--- a/src/bn.c
+++ b/src/bn.c
@@ -1,5 +1,5 @@
 /*
- * bn.c -- 256-bit and 512-bit bignum calculation
+ * bn.c -- 256-bit (and 512-bit) bignum calculation
  *
  * Copyright (C) 2011 Free Software Initiative of Japan
  * Author: NIIBE Yutaka <gniibe@fsij.org>
diff --git a/src/ec_p256.c b/src/ec_p256.c
index 2b8a031..38cf738 100644
--- a/src/ec_p256.c
+++ b/src/ec_p256.c
@@ -1,5 +1,37 @@
 /*
  * ec_p256.c - Elliptic curve over GF(p256)
+ *
+ * Copyright (C) 2011 Free Software Initiative of Japan
+ * Author: NIIBE Yutaka <gniibe@fsij.org>
+ *
+ * This file is a part of Gnuk, a GnuPG USB Token implementation.
+ *
+ * Gnuk is free software: you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Gnuk 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/>.
+ *
+ */
+
+/*
+ * References:
+ *
+ * [1] Suite B Implementer's Guide to FIPS 186-3 (ECDSA), February 3, 2010.
+ *
+ * [2] Michael Brown, Darrel Hankerson, Julio López, and Alfred Menezes,
+ *     Software Implementation of the NIST Elliptic Curves Over Prime Fields,
+ *     Proceedings of the 2001 Conference on Topics in Cryptology: The
+ *     Cryptographer's Track at RSA
+ *     Pages 250-265, Springer-Verlag London, UK, 2001
+ *     ISBN:3-540-41898-9
  */
 
 #include <stdint.h>
@@ -10,7 +42,7 @@
 #include "mod.h"
 #include "ec_p256.h"
 
-#if 0
+#if TEST
 /*
  * Generator of Elliptic curve over GF(p256)
  */
@@ -25,6 +57,7 @@ bn256 Gy[1] = {
 };
 #endif
 
+
 /*
  * w = 4
  * m = 256
@@ -190,18 +223,23 @@ static const ac precomputed_2E_KG[15] = {
   }
 };
 
+/**
+ * @brief	X  = k * G
+ *
+ * @param K	scalar k
+ */
 void
 compute_kG (ac *X, const bn256 *K)
 {
   int i;
-  int q_infinite = 1;
+  int q_is_infinite = 1;
   jpc Q[1];
 
   for (i = 31; i >= 0; i--)
     {
       int k_i, k_i_e;
 
-      if (!q_infinite)
+      if (!q_is_infinite)
 	jpc_double (Q, Q);
 
       k_i = (((K->words[6] >> i) & 1) << 3)
@@ -215,7 +253,7 @@ compute_kG (ac *X, const bn256 *K)
 
       if (k_i)
 	{
-	  if (q_infinite)
+	  if (q_is_infinite)
 	    {
 	      memcpy (Q->x, (&precomputed_KG[k_i - 1])->x, sizeof (bn256));
 	      memcpy (Q->y, (&precomputed_KG[k_i - 1])->y, sizeof (bn256));
@@ -223,20 +261,20 @@ compute_kG (ac *X, const bn256 *K)
 	      Q->z->words[1] = Q->z->words[2] = Q->z->words[3]
 		= Q->z->words[4] = Q->z->words[5] = Q->z->words[6]
 		= Q->z->words[7] = 0;
-	      q_infinite = 0;
+	      q_is_infinite = 0;
 	    }
 	  else
 	    jpc_add_ac (Q, Q, &precomputed_KG[k_i - 1]);
 	}
       if (k_i_e)
 	{
-	  if (q_infinite)
+	  if (q_is_infinite)
 	    {
 	      memcpy (Q->x, (&precomputed_2E_KG[k_i_e - 1])->x, sizeof (bn256));
 	      memcpy (Q->y, (&precomputed_2E_KG[k_i_e - 1])->y, sizeof (bn256));
 	      memset (Q->z, 0, sizeof (bn256));
 	      Q->z->words[0] = 1;
-	      q_infinite = 0;
+	      q_is_infinite = 0;
 	    }
 	  else
 	    jpc_add_ac (Q, Q, &precomputed_2E_KG[k_i_e - 1]);
@@ -250,7 +288,7 @@ compute_kG (ac *X, const bn256 *K)
 #define NAF_K_SIGN(k)	(k&8)
 #define NAF_K_INDEX(k)	((k&7)-1)
 
-void
+static void
 naf4_257_set (naf4_257 *NAF_K, int i, int ki)
 {
   if (ki != 0)
@@ -270,7 +308,7 @@ naf4_257_set (naf4_257 *NAF_K, int i, int ki)
     }
 }
 
-int
+static int
 naf4_257_get (const naf4_257 *NAF_K, int i)
 {
   int ki;
@@ -286,6 +324,10 @@ naf4_257_get (const naf4_257 *NAF_K, int i)
   return ki;
 }
 
+
+/*
+ * convert 256-bit bignum into non-adjacent form (NAF)
+ */
 void
 compute_naf4_257 (naf4_257 *NAF_K, const bn256 *K)
 {
@@ -325,11 +367,17 @@ compute_naf4_257 (naf4_257 *NAF_K, const bn256 *K)
     }
 }
 
+/**
+ * @brief	X  = k * P
+ *
+ * @param NAF_K	NAF representation of k
+ * @param P	P in affin coordiate
+ */
 void
 compute_kP (ac *X, const naf4_257 *NAF_K, const ac *P)
 {
   int i;
-  int q_infinite = 1;
+  int q_is_infinite = 1;
   jpc Q[1];
   ac P3[1], P5[1], P7[1];
   const ac *p_Pi[4];
@@ -367,13 +415,13 @@ compute_kP (ac *X, const naf4_257 *NAF_K, const ac *P)
     {
       int k_i;
 
-      if (!q_infinite)
+      if (!q_is_infinite)
 	jpc_double (Q, Q);
 
       k_i = naf4_257_get (NAF_K, i);
       if (k_i)
 	{
-	  if (q_infinite)
+	  if (q_is_infinite)
 	    {
 	      memcpy (Q->x, p_Pi[NAF_K_INDEX(k_i)]->x, sizeof (bn256));
 	      if (NAF_K_SIGN (k_i))
@@ -382,7 +430,7 @@ compute_kP (ac *X, const naf4_257 *NAF_K, const ac *P)
 		memcpy (Q->y, p_Pi[NAF_K_INDEX(k_i)]->y, sizeof (bn256));
 	      memset (Q->z, 0, sizeof (bn256));
 	      Q->z->words[0] = 1;
-	      q_infinite = 0;
+	      q_is_infinite = 0;
 	    }
 	  else
 	    jpc_add_ac_signed (Q, Q, p_Pi[NAF_K_INDEX(k_i)], NAF_K_SIGN (k_i));
@@ -412,7 +460,7 @@ static const bn256 MU_lower[1] = {
 
 
 /**
- * @brief Compute signature (r,s) of hash string z with secret key dA
+ * @brief Compute signature (r,s) of hash string z with secret key d
  */
 void
 ecdsa (bn256 *r, bn256 *s, const bn256 *z, const bn256 *d)
@@ -423,11 +471,18 @@ ecdsa (bn256 *r, bn256 *s, const bn256 *z, const bn256 *d)
   bn256 k_inv[1];
   uint32_t carry;
 
+#define tmp_k k_inv
+
   do
     {
       do
 	{
 	  bn256_random (k);
+	  if (bn256_sub (tmp_k, k, N) == 0)	/* > N, it's too big.  */
+	    continue;
+	  if (bn256_add_uint (tmp_k, tmp_k, 2)) /* > N - 2, still big.  */
+	    continue;
+	  bn256_add_uint (k, 1);
 	  compute_kG (KG, k);
 	  if (bn256_is_ge (KG->x, N))
 	    bn256_sub (r, KG->x, N);
diff --git a/src/ec_p256.h b/src/ec_p256.h
index 3e4c199..0c465fd 100644
--- a/src/ec_p256.h
+++ b/src/ec_p256.h
@@ -1,3 +1,4 @@
+/* Non-adjacent form */
 typedef struct naf4_257 {
   uint32_t words[ BN256_WORDS*4 ]; /* Little endian */
   uint8_t last_nibble;		   /* most significant nibble */
diff --git a/src/ecc-cdh.c b/src/ecc-cdh.c
new file mode 100644
index 0000000..c792094
--- /dev/null
+++ b/src/ecc-cdh.c
@@ -0,0 +1,68 @@
+/*
+ * ecc-cdh.c - One-Pass Diffie-Hellman method implementation
+ *             C(1, 1, ECC CDH) for EC DH of OpenPGP ECC
+ *
+ * Copyright (C) 2013 Free Software Initiative of Japan
+ * Author: NIIBE Yutaka <gniibe@fsij.org>
+ *
+ * This file is a part of Gnuk, a GnuPG USB Token implementation.
+ *
+ * Gnuk is free software: you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Gnuk 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/>.
+ *
+ */
+
+/*
+ * References:
+ *
+ * [1] A. Jivsov, Elliptic Curve Cryptography (ECC) in OpenPGP, RFC 6637,
+ *     June 2012.
+ *
+ * [2] Suite B Implementer's Guide to NIST SP 800-56A, July 28, 2009.
+ *
+ */
+
+static const char param[] = {
+  /**/
+  curve_OID_len,
+  curve_OID,
+  public-key_alg_ID, /*ecdh*/
+  0x03,
+  0x01,
+  KDF_hash_ID, /*sha256*/
+  KEK_alg_ID, /*aes128*/
+  "Anonymous Sender    ",
+  my_finger_print /*20-byte*/
+};
+
+/*
+ *
+ */
+int
+ecdh (unsigned char *key,
+      const unsigned char *key_encrypted, const ac *P,
+      const naf4_257 *naf_d, const unsigned char *fp)
+{
+  ac S[1];
+  sha256_context ctx;
+  unsigned char kek[32];
+
+  compute_kP (S, naf_d, P);	/* Get shared key.  */
+
+  /* kdf (kek, S, parameter) */
+  sha256_start (&ctx);
+  sha256_update (&ctx, "\x00\x00\x00\x01", 4);
+  sha256_update (&ctx, (const char *)S, size of S); /* XXX 04, X, Y bigendian!! */
+  sha256_update (&ctx, (const char *)param, size of param);
+  sha256_finish (&ctx, kek);
+}
diff --git a/src/mod.h b/src/mod.h
index 8f4b9c1..117bcc7 100644
--- a/src/mod.h
+++ b/src/mod.h
@@ -1,2 +1,3 @@
-void mod_reduce (bn256 *X, const bn512 *A, const bn256 *B, const bn256 *MU_lower);
+void mod_reduce (bn256 *X, const bn512 *A, const bn256 *B,
+		 const bn256 *MU_lower);
 void mod_inv (bn256 *X, const bn256 *A, const bn256 *N);