diff --git a/jni/libzrtp/sources/bnlib/ec/curve25519-donna.c b/jni/libzrtp/sources/bnlib/ec/curve25519-donna.c
new file mode 100644
index 0000000..de11280
--- /dev/null
+++ b/jni/libzrtp/sources/bnlib/ec/curve25519-donna.c
@@ -0,0 +1,731 @@
+/* Copyright 2008, Google Inc.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are
+ * met:
+ *
+ *     * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *     * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following disclaimer
+ * in the documentation and/or other materials provided with the
+ * distribution.
+ *     * Neither the name of Google Inc. nor the names of its
+ * contributors may be used to endorse or promote products derived from
+ * this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+ * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+ * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+ * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * curve25519-donna: Curve25519 elliptic curve, public key function
+ *
+ * http://code.google.com/p/curve25519-donna/
+ *
+ * Adam Langley <agl@imperialviolet.org>
+ *
+ * Derived from public domain C code by Daniel J. Bernstein <djb@cr.yp.to>
+ *
+ * More information about curve25519 can be found here
+ *   http://cr.yp.to/ecdh.html
+ *
+ * djb's sample implementation of curve25519 is written in a special assembly
+ * language called qhasm and uses the floating point registers.
+ *
+ * This is, almost, a clean room reimplementation from the curve25519 paper. It
+ * uses many of the tricks described therein. Only the crecip function is taken
+ * from the sample implementation.
+ */
+
+#include <string.h>
+#include <stdint.h>
+
+#ifdef _MSC_VER
+#define inline __inline
+#endif
+
+typedef uint8_t u8;
+typedef int32_t s32;
+typedef int64_t limb;
+
+/* Field element representation:
+ *
+ * Field elements are written as an array of signed, 64-bit limbs, least
+ * significant first. The value of the field element is:
+ *   x[0] + 2^26·x[1] + x^51·x[2] + 2^102·x[3] + ...
+ *
+ * i.e. the limbs are 26, 25, 26, 25, ... bits wide.
+ */
+
+/* Sum two numbers: output += in */
+static void fsum(limb *output, const limb *in) {
+  unsigned i;
+  for (i = 0; i < 10; i += 2) {
+    output[0+i] = (output[0+i] + in[0+i]);
+    output[1+i] = (output[1+i] + in[1+i]);
+  }
+}
+
+/* Find the difference of two numbers: output = in - output
+ * (note the order of the arguments!)
+ */
+static void fdifference(limb *output, const limb *in) {
+  unsigned i;
+  for (i = 0; i < 10; ++i) {
+    output[i] = (in[i] - output[i]);
+  }
+}
+
+/* Multiply a number by a scalar: output = in * scalar */
+static void fscalar_product(limb *output, const limb *in, const limb scalar) {
+  unsigned i;
+  for (i = 0; i < 10; ++i) {
+    output[i] = in[i] * scalar;
+  }
+}
+
+/* Multiply two numbers: output = in2 * in
+ *
+ * output must be distinct to both inputs. The inputs are reduced coefficient
+ * form, the output is not.
+ */
+static void fproduct(limb *output, const limb *in2, const limb *in) {
+  output[0] =       ((limb) ((s32) in2[0])) * ((s32) in[0]);
+  output[1] =       ((limb) ((s32) in2[0])) * ((s32) in[1]) +
+                    ((limb) ((s32) in2[1])) * ((s32) in[0]);
+  output[2] =  2 *  ((limb) ((s32) in2[1])) * ((s32) in[1]) +
+                    ((limb) ((s32) in2[0])) * ((s32) in[2]) +
+                    ((limb) ((s32) in2[2])) * ((s32) in[0]);
+  output[3] =       ((limb) ((s32) in2[1])) * ((s32) in[2]) +
+                    ((limb) ((s32) in2[2])) * ((s32) in[1]) +
+                    ((limb) ((s32) in2[0])) * ((s32) in[3]) +
+                    ((limb) ((s32) in2[3])) * ((s32) in[0]);
+  output[4] =       ((limb) ((s32) in2[2])) * ((s32) in[2]) +
+               2 * (((limb) ((s32) in2[1])) * ((s32) in[3]) +
+                    ((limb) ((s32) in2[3])) * ((s32) in[1])) +
+                    ((limb) ((s32) in2[0])) * ((s32) in[4]) +
+                    ((limb) ((s32) in2[4])) * ((s32) in[0]);
+  output[5] =       ((limb) ((s32) in2[2])) * ((s32) in[3]) +
+                    ((limb) ((s32) in2[3])) * ((s32) in[2]) +
+                    ((limb) ((s32) in2[1])) * ((s32) in[4]) +
+                    ((limb) ((s32) in2[4])) * ((s32) in[1]) +
+                    ((limb) ((s32) in2[0])) * ((s32) in[5]) +
+                    ((limb) ((s32) in2[5])) * ((s32) in[0]);
+  output[6] =  2 * (((limb) ((s32) in2[3])) * ((s32) in[3]) +
+                    ((limb) ((s32) in2[1])) * ((s32) in[5]) +
+                    ((limb) ((s32) in2[5])) * ((s32) in[1])) +
+                    ((limb) ((s32) in2[2])) * ((s32) in[4]) +
+                    ((limb) ((s32) in2[4])) * ((s32) in[2]) +
+                    ((limb) ((s32) in2[0])) * ((s32) in[6]) +
+                    ((limb) ((s32) in2[6])) * ((s32) in[0]);
+  output[7] =       ((limb) ((s32) in2[3])) * ((s32) in[4]) +
+                    ((limb) ((s32) in2[4])) * ((s32) in[3]) +
+                    ((limb) ((s32) in2[2])) * ((s32) in[5]) +
+                    ((limb) ((s32) in2[5])) * ((s32) in[2]) +
+                    ((limb) ((s32) in2[1])) * ((s32) in[6]) +
+                    ((limb) ((s32) in2[6])) * ((s32) in[1]) +
+                    ((limb) ((s32) in2[0])) * ((s32) in[7]) +
+                    ((limb) ((s32) in2[7])) * ((s32) in[0]);
+  output[8] =       ((limb) ((s32) in2[4])) * ((s32) in[4]) +
+               2 * (((limb) ((s32) in2[3])) * ((s32) in[5]) +
+                    ((limb) ((s32) in2[5])) * ((s32) in[3]) +
+                    ((limb) ((s32) in2[1])) * ((s32) in[7]) +
+                    ((limb) ((s32) in2[7])) * ((s32) in[1])) +
+                    ((limb) ((s32) in2[2])) * ((s32) in[6]) +
+                    ((limb) ((s32) in2[6])) * ((s32) in[2]) +
+                    ((limb) ((s32) in2[0])) * ((s32) in[8]) +
+                    ((limb) ((s32) in2[8])) * ((s32) in[0]);
+  output[9] =       ((limb) ((s32) in2[4])) * ((s32) in[5]) +
+                    ((limb) ((s32) in2[5])) * ((s32) in[4]) +
+                    ((limb) ((s32) in2[3])) * ((s32) in[6]) +
+                    ((limb) ((s32) in2[6])) * ((s32) in[3]) +
+                    ((limb) ((s32) in2[2])) * ((s32) in[7]) +
+                    ((limb) ((s32) in2[7])) * ((s32) in[2]) +
+                    ((limb) ((s32) in2[1])) * ((s32) in[8]) +
+                    ((limb) ((s32) in2[8])) * ((s32) in[1]) +
+                    ((limb) ((s32) in2[0])) * ((s32) in[9]) +
+                    ((limb) ((s32) in2[9])) * ((s32) in[0]);
+  output[10] = 2 * (((limb) ((s32) in2[5])) * ((s32) in[5]) +
+                    ((limb) ((s32) in2[3])) * ((s32) in[7]) +
+                    ((limb) ((s32) in2[7])) * ((s32) in[3]) +
+                    ((limb) ((s32) in2[1])) * ((s32) in[9]) +
+                    ((limb) ((s32) in2[9])) * ((s32) in[1])) +
+                    ((limb) ((s32) in2[4])) * ((s32) in[6]) +
+                    ((limb) ((s32) in2[6])) * ((s32) in[4]) +
+                    ((limb) ((s32) in2[2])) * ((s32) in[8]) +
+                    ((limb) ((s32) in2[8])) * ((s32) in[2]);
+  output[11] =      ((limb) ((s32) in2[5])) * ((s32) in[6]) +
+                    ((limb) ((s32) in2[6])) * ((s32) in[5]) +
+                    ((limb) ((s32) in2[4])) * ((s32) in[7]) +
+                    ((limb) ((s32) in2[7])) * ((s32) in[4]) +
+                    ((limb) ((s32) in2[3])) * ((s32) in[8]) +
+                    ((limb) ((s32) in2[8])) * ((s32) in[3]) +
+                    ((limb) ((s32) in2[2])) * ((s32) in[9]) +
+                    ((limb) ((s32) in2[9])) * ((s32) in[2]);
+  output[12] =      ((limb) ((s32) in2[6])) * ((s32) in[6]) +
+               2 * (((limb) ((s32) in2[5])) * ((s32) in[7]) +
+                    ((limb) ((s32) in2[7])) * ((s32) in[5]) +
+                    ((limb) ((s32) in2[3])) * ((s32) in[9]) +
+                    ((limb) ((s32) in2[9])) * ((s32) in[3])) +
+                    ((limb) ((s32) in2[4])) * ((s32) in[8]) +
+                    ((limb) ((s32) in2[8])) * ((s32) in[4]);
+  output[13] =      ((limb) ((s32) in2[6])) * ((s32) in[7]) +
+                    ((limb) ((s32) in2[7])) * ((s32) in[6]) +
+                    ((limb) ((s32) in2[5])) * ((s32) in[8]) +
+                    ((limb) ((s32) in2[8])) * ((s32) in[5]) +
+                    ((limb) ((s32) in2[4])) * ((s32) in[9]) +
+                    ((limb) ((s32) in2[9])) * ((s32) in[4]);
+  output[14] = 2 * (((limb) ((s32) in2[7])) * ((s32) in[7]) +
+                    ((limb) ((s32) in2[5])) * ((s32) in[9]) +
+                    ((limb) ((s32) in2[9])) * ((s32) in[5])) +
+                    ((limb) ((s32) in2[6])) * ((s32) in[8]) +
+                    ((limb) ((s32) in2[8])) * ((s32) in[6]);
+  output[15] =      ((limb) ((s32) in2[7])) * ((s32) in[8]) +
+                    ((limb) ((s32) in2[8])) * ((s32) in[7]) +
+                    ((limb) ((s32) in2[6])) * ((s32) in[9]) +
+                    ((limb) ((s32) in2[9])) * ((s32) in[6]);
+  output[16] =      ((limb) ((s32) in2[8])) * ((s32) in[8]) +
+               2 * (((limb) ((s32) in2[7])) * ((s32) in[9]) +
+                    ((limb) ((s32) in2[9])) * ((s32) in[7]));
+  output[17] =      ((limb) ((s32) in2[8])) * ((s32) in[9]) +
+                    ((limb) ((s32) in2[9])) * ((s32) in[8]);
+  output[18] = 2 *  ((limb) ((s32) in2[9])) * ((s32) in[9]);
+}
+
+/* Reduce a long form to a short form by taking the input mod 2^255 - 19. */
+static void freduce_degree(limb *output) {
+  /* Each of these shifts and adds ends up multiplying the value by 19. */
+  output[8] += output[18] << 4;
+  output[8] += output[18] << 1;
+  output[8] += output[18];
+  output[7] += output[17] << 4;
+  output[7] += output[17] << 1;
+  output[7] += output[17];
+  output[6] += output[16] << 4;
+  output[6] += output[16] << 1;
+  output[6] += output[16];
+  output[5] += output[15] << 4;
+  output[5] += output[15] << 1;
+  output[5] += output[15];
+  output[4] += output[14] << 4;
+  output[4] += output[14] << 1;
+  output[4] += output[14];
+  output[3] += output[13] << 4;
+  output[3] += output[13] << 1;
+  output[3] += output[13];
+  output[2] += output[12] << 4;
+  output[2] += output[12] << 1;
+  output[2] += output[12];
+  output[1] += output[11] << 4;
+  output[1] += output[11] << 1;
+  output[1] += output[11];
+  output[0] += output[10] << 4;
+  output[0] += output[10] << 1;
+  output[0] += output[10];
+}
+
+#if (-1 & 3) != 3
+#error "This code only works on a two's complement system"
+#endif
+
+/* return v / 2^26, using only shifts and adds. */
+static limb div_by_2_26(const limb v)
+{
+  /* High word of v; no shift needed*/
+  const uint32_t highword = (uint32_t) (((uint64_t) v) >> 32);
+  /* Set to all 1s if v was negative; else set to 0s. */
+  const int32_t sign = ((int32_t) highword) >> 31;
+  /* Set to 0x3ffffff if v was negative; else set to 0. */
+  const int32_t roundoff = ((uint32_t) sign) >> 6;
+  /* Should return v / (1<<26) */
+  return (v + roundoff) >> 26;
+}
+
+/* return v / (2^25), using only shifts and adds. */
+static limb div_by_2_25(const limb v)
+{
+  /* High word of v; no shift needed*/
+  const uint32_t highword = (uint32_t) (((uint64_t) v) >> 32);
+  /* Set to all 1s if v was negative; else set to 0s. */
+  const int32_t sign = ((int32_t) highword) >> 31;
+  /* Set to 0x1ffffff if v was negative; else set to 0. */
+  const int32_t roundoff = ((uint32_t) sign) >> 7;
+  /* Should return v / (1<<25) */
+  return (v + roundoff) >> 25;
+}
+
+static s32 div_s32_by_2_25(const s32 v)
+{
+   const s32 roundoff = ((uint32_t)(v >> 31)) >> 7;
+   return (v + roundoff) >> 25;
+}
+
+/* Reduce all coefficients of the short form input so that |x| < 2^26.
+ *
+ * On entry: |output[i]| < 2^62
+ */
+static void freduce_coefficients(limb *output) {
+  unsigned i;
+
+  output[10] = 0;
+
+  for (i = 0; i < 10; i += 2) {
+    limb over = div_by_2_26(output[i]);
+    output[i] -= over << 26;
+    output[i+1] += over;
+
+    over = div_by_2_25(output[i+1]);
+    output[i+1] -= over << 25;
+    output[i+2] += over;
+  }
+  /* Now |output[10]| < 2 ^ 38 and all other coefficients are reduced. */
+  output[0] += output[10] << 4;
+  output[0] += output[10] << 1;
+  output[0] += output[10];
+
+  output[10] = 0;
+
+  /* Now output[1..9] are reduced, and |output[0]| < 2^26 + 19 * 2^38
+   * So |over| will be no more than 77825  */
+  {
+    limb over = div_by_2_26(output[0]);
+    output[0] -= over << 26;
+    output[1] += over;
+  }
+
+  /* Now output[0,2..9] are reduced, and |output[1]| < 2^25 + 77825
+   * So |over| will be no more than 1. */
+  {
+    /* output[1] fits in 32 bits, so we can use div_s32_by_2_25 here. */
+    s32 over32 = div_s32_by_2_25((s32) output[1]);
+    output[1] -= over32 << 25;
+    output[2] += over32;
+  }
+
+  /* Finally, output[0,1,3..9] are reduced, and output[2] is "nearly reduced":
+   * we have |output[2]| <= 2^26.  This is good enough for all of our math,
+   * but it will require an extra freduce_coefficients before fcontract. */
+}
+
+/* A helpful wrapper around fproduct: output = in * in2.
+ *
+ * output must be distinct to both inputs. The output is reduced degree and
+ * reduced coefficient.
+ */
+static void
+fmul(limb *output, const limb *in, const limb *in2) {
+  limb t[19];
+  fproduct(t, in, in2);
+  freduce_degree(t);
+  freduce_coefficients(t);
+  memcpy(output, t, sizeof(limb) * 10);
+}
+
+static void fsquare_inner(limb *output, const limb *in) {
+  output[0] =       ((limb) ((s32) in[0])) * ((s32) in[0]);
+  output[1] =  2 *  ((limb) ((s32) in[0])) * ((s32) in[1]);
+  output[2] =  2 * (((limb) ((s32) in[1])) * ((s32) in[1]) +
+                    ((limb) ((s32) in[0])) * ((s32) in[2]));
+  output[3] =  2 * (((limb) ((s32) in[1])) * ((s32) in[2]) +
+                    ((limb) ((s32) in[0])) * ((s32) in[3]));
+  output[4] =       ((limb) ((s32) in[2])) * ((s32) in[2]) +
+               4 *  ((limb) ((s32) in[1])) * ((s32) in[3]) +
+               2 *  ((limb) ((s32) in[0])) * ((s32) in[4]);
+  output[5] =  2 * (((limb) ((s32) in[2])) * ((s32) in[3]) +
+                    ((limb) ((s32) in[1])) * ((s32) in[4]) +
+                    ((limb) ((s32) in[0])) * ((s32) in[5]));
+  output[6] =  2 * (((limb) ((s32) in[3])) * ((s32) in[3]) +
+                    ((limb) ((s32) in[2])) * ((s32) in[4]) +
+                    ((limb) ((s32) in[0])) * ((s32) in[6]) +
+               2 *  ((limb) ((s32) in[1])) * ((s32) in[5]));
+  output[7] =  2 * (((limb) ((s32) in[3])) * ((s32) in[4]) +
+                    ((limb) ((s32) in[2])) * ((s32) in[5]) +
+                    ((limb) ((s32) in[1])) * ((s32) in[6]) +
+                    ((limb) ((s32) in[0])) * ((s32) in[7]));
+  output[8] =       ((limb) ((s32) in[4])) * ((s32) in[4]) +
+               2 * (((limb) ((s32) in[2])) * ((s32) in[6]) +
+                    ((limb) ((s32) in[0])) * ((s32) in[8]) +
+               2 * (((limb) ((s32) in[1])) * ((s32) in[7]) +
+                    ((limb) ((s32) in[3])) * ((s32) in[5])));
+  output[9] =  2 * (((limb) ((s32) in[4])) * ((s32) in[5]) +
+                    ((limb) ((s32) in[3])) * ((s32) in[6]) +
+                    ((limb) ((s32) in[2])) * ((s32) in[7]) +
+                    ((limb) ((s32) in[1])) * ((s32) in[8]) +
+                    ((limb) ((s32) in[0])) * ((s32) in[9]));
+  output[10] = 2 * (((limb) ((s32) in[5])) * ((s32) in[5]) +
+                    ((limb) ((s32) in[4])) * ((s32) in[6]) +
+                    ((limb) ((s32) in[2])) * ((s32) in[8]) +
+               2 * (((limb) ((s32) in[3])) * ((s32) in[7]) +
+                    ((limb) ((s32) in[1])) * ((s32) in[9])));
+  output[11] = 2 * (((limb) ((s32) in[5])) * ((s32) in[6]) +
+                    ((limb) ((s32) in[4])) * ((s32) in[7]) +
+                    ((limb) ((s32) in[3])) * ((s32) in[8]) +
+                    ((limb) ((s32) in[2])) * ((s32) in[9]));
+  output[12] =      ((limb) ((s32) in[6])) * ((s32) in[6]) +
+               2 * (((limb) ((s32) in[4])) * ((s32) in[8]) +
+               2 * (((limb) ((s32) in[5])) * ((s32) in[7]) +
+                    ((limb) ((s32) in[3])) * ((s32) in[9])));
+  output[13] = 2 * (((limb) ((s32) in[6])) * ((s32) in[7]) +
+                    ((limb) ((s32) in[5])) * ((s32) in[8]) +
+                    ((limb) ((s32) in[4])) * ((s32) in[9]));
+  output[14] = 2 * (((limb) ((s32) in[7])) * ((s32) in[7]) +
+                    ((limb) ((s32) in[6])) * ((s32) in[8]) +
+               2 *  ((limb) ((s32) in[5])) * ((s32) in[9]));
+  output[15] = 2 * (((limb) ((s32) in[7])) * ((s32) in[8]) +
+                    ((limb) ((s32) in[6])) * ((s32) in[9]));
+  output[16] =      ((limb) ((s32) in[8])) * ((s32) in[8]) +
+               4 *  ((limb) ((s32) in[7])) * ((s32) in[9]);
+  output[17] = 2 *  ((limb) ((s32) in[8])) * ((s32) in[9]);
+  output[18] = 2 *  ((limb) ((s32) in[9])) * ((s32) in[9]);
+}
+
+static void
+fsquare(limb *output, const limb *in) {
+  limb t[19];
+  fsquare_inner(t, in);
+  freduce_degree(t);
+  freduce_coefficients(t);
+  memcpy(output, t, sizeof(limb) * 10);
+}
+
+/* Take a little-endian, 32-byte number and expand it into polynomial form */
+static void
+fexpand(limb *output, const u8 *input) {
+#define F(n,start,shift,mask) \
+  output[n] = ((((limb) input[start + 0]) | \
+                ((limb) input[start + 1]) << 8 | \
+                ((limb) input[start + 2]) << 16 | \
+                ((limb) input[start + 3]) << 24) >> shift) & mask;
+  F(0, 0, 0, 0x3ffffff);
+  F(1, 3, 2, 0x1ffffff);
+  F(2, 6, 3, 0x3ffffff);
+  F(3, 9, 5, 0x1ffffff);
+  F(4, 12, 6, 0x3ffffff);
+  F(5, 16, 0, 0x1ffffff);
+  F(6, 19, 1, 0x3ffffff);
+  F(7, 22, 3, 0x1ffffff);
+  F(8, 25, 4, 0x3ffffff);
+  F(9, 28, 6, 0x1ffffff);
+#undef F
+}
+
+#if (-32 >> 1) != -16
+#error "This code only works when >> does sign-extension on negative numbers"
+#endif
+
+/* Take a fully reduced polynomial form number and contract it into a
+ * little-endian, 32-byte array
+ */
+static void
+fcontract(u8 *output, limb *input) {
+  int i;
+  int j;
+
+  for (j = 0; j < 2; ++j) {
+    for (i = 0; i < 9; ++i) {
+      if ((i & 1) == 1) {
+        /* This calculation is a time-invariant way to make input[i] positive
+           by borrowing from the next-larger limb.
+        */
+        const s32 mask = (s32)(input[i]) >> 31;
+        const s32 carry = -(((s32)(input[i]) & mask) >> 25);
+        input[i] = (s32)(input[i]) + (carry << 25);
+        input[i+1] = (s32)(input[i+1]) - carry;
+      } else {
+        const s32 mask = (s32)(input[i]) >> 31;
+        const s32 carry = -(((s32)(input[i]) & mask) >> 26);
+        input[i] = (s32)(input[i]) + (carry << 26);
+        input[i+1] = (s32)(input[i+1]) - carry;
+      }
+    }
+    {
+      const s32 mask = (s32)(input[9]) >> 31;
+      const s32 carry = -(((s32)(input[9]) & mask) >> 25);
+      input[9] = (s32)(input[9]) + (carry << 25);
+      input[0] = (s32)(input[0]) - (carry * 19);
+    }
+  }
+
+  /* The first borrow-propagation pass above ended with every limb
+     except (possibly) input[0] non-negative.
+
+     Since each input limb except input[0] is decreased by at most 1
+     by a borrow-propagation pass, the second borrow-propagation pass
+     could only have wrapped around to decrease input[0] again if the
+     first pass left input[0] negative *and* input[1] through input[9]
+     were all zero.  In that case, input[1] is now 2^25 - 1, and this
+     last borrow-propagation step will leave input[1] non-negative.
+  */
+  {
+    const s32 mask = (s32)(input[0]) >> 31;
+    const s32 carry = -(((s32)(input[0]) & mask) >> 26);
+    input[0] = (s32)(input[0]) + (carry << 26);
+    input[1] = (s32)(input[1]) - carry;
+  }
+
+  /* Both passes through the above loop, plus the last 0-to-1 step, are
+     necessary: if input[9] is -1 and input[0] through input[8] are 0,
+     negative values will remain in the array until the end.
+   */
+
+  input[1] <<= 2;
+  input[2] <<= 3;
+  input[3] <<= 5;
+  input[4] <<= 6;
+  input[6] <<= 1;
+  input[7] <<= 3;
+  input[8] <<= 4;
+  input[9] <<= 6;
+#define F(i, s) \
+  output[s+0] |=  input[i] & 0xff; \
+  output[s+1]  = (input[i] >> 8) & 0xff; \
+  output[s+2]  = (input[i] >> 16) & 0xff; \
+  output[s+3]  = (input[i] >> 24) & 0xff;
+  output[0] = 0;
+  output[16] = 0;
+  F(0,0);
+  F(1,3);
+  F(2,6);
+  F(3,9);
+  F(4,12);
+  F(5,16);
+  F(6,19);
+  F(7,22);
+  F(8,25);
+  F(9,28);
+#undef F
+}
+
+/* Input: Q, Q', Q-Q'
+ * Output: 2Q, Q+Q'
+ *
+ *   x2 z3: long form
+ *   x3 z3: long form
+ *   x z: short form, destroyed
+ *   xprime zprime: short form, destroyed
+ *   qmqp: short form, preserved
+ */
+static void fmonty(limb *x2, limb *z2,  /* output 2Q */
+                   limb *x3, limb *z3,  /* output Q + Q' */
+                   limb *x, limb *z,    /* input Q */
+                   limb *xprime, limb *zprime,  /* input Q' */
+                   const limb *qmqp /* input Q - Q' */) {
+  limb origx[10], origxprime[10], zzz[19], xx[19], zz[19], xxprime[19],
+        zzprime[19], zzzprime[19], xxxprime[19];
+
+  memcpy(origx, x, 10 * sizeof(limb));
+  fsum(x, z);
+  fdifference(z, origx);  /* does x - z */
+
+  memcpy(origxprime, xprime, sizeof(limb) * 10);
+  fsum(xprime, zprime);
+  fdifference(zprime, origxprime);
+  fproduct(xxprime, xprime, z);
+  fproduct(zzprime, x, zprime);
+  freduce_degree(xxprime);
+  freduce_coefficients(xxprime);
+  freduce_degree(zzprime);
+  freduce_coefficients(zzprime);
+  memcpy(origxprime, xxprime, sizeof(limb) * 10);
+  fsum(xxprime, zzprime);
+  fdifference(zzprime, origxprime);
+  fsquare(xxxprime, xxprime);
+  fsquare(zzzprime, zzprime);
+  fproduct(zzprime, zzzprime, qmqp);
+  freduce_degree(zzprime);
+  freduce_coefficients(zzprime);
+  memcpy(x3, xxxprime, sizeof(limb) * 10);
+  memcpy(z3, zzprime, sizeof(limb) * 10);
+
+  fsquare(xx, x);
+  fsquare(zz, z);
+  fproduct(x2, xx, zz);
+  freduce_degree(x2);
+  freduce_coefficients(x2);
+  fdifference(zz, xx);         /* does zz = xx - zz */
+  memset(zzz + 10, 0, sizeof(limb) * 9);
+  fscalar_product(zzz, zz, 121665);
+  /* No need to call freduce_degree here:
+     fscalar_product doesn't increase the degree of its input.
+   */
+  freduce_coefficients(zzz);
+  fsum(zzz, xx);
+  fproduct(z2, zz, zzz);
+  freduce_degree(z2);
+  freduce_coefficients(z2);
+}
+
+/* Conditionally swap two reduced-form limb arrays if 'iswap' is 1, but leave
+ * them unchanged if 'iswap' is 0.  Runs in data-invariant time to avoid
+ * side-channel attacks.
+ *
+ * NOTE that this function requires that 'iswap' be 1 or 0; other values give
+ * wrong results.  Also, the two limb arrays must be in reduced-coefficient,
+ * reduced-degree form: the values in a[10..19] or b[10..19] aren't swapped,
+ * and all all values in a[0..9],b[0..9] must have magnitude less than
+ * INT32_MAX.
+ */
+static void
+swap_conditional(limb a[19], limb b[19], limb iswap) {
+  unsigned i;
+  const s32 swap = (s32) -iswap;
+
+  for (i = 0; i < 10; ++i) {
+    const s32 x = swap & ( ((s32)a[i]) ^ ((s32)b[i]) );
+    a[i] = ((s32)a[i]) ^ x;
+    b[i] = ((s32)b[i]) ^ x;
+  }
+}
+
+/* Calculates nQ where Q is the x-coordinate of a point on the curve
+ *
+ *   resultx/resultz: the x coordinate of the resulting curve point (short form)
+ *   n: a little endian, 32-byte number
+ *   q: a point of the curve (short form)
+ */
+static void
+cmult(limb *resultx, limb *resultz, const u8 *n, const limb *q) {
+  limb a[19] = {0}, b[19] = {1}, c[19] = {1}, d[19] = {0};
+  limb *nqpqx = a, *nqpqz = b, *nqx = c, *nqz = d, *t;
+  limb e[19] = {0}, f[19] = {1}, g[19] = {0}, h[19] = {1};
+  limb *nqpqx2 = e, *nqpqz2 = f, *nqx2 = g, *nqz2 = h;
+
+  unsigned i, j;
+
+  memcpy(nqpqx, q, sizeof(limb) * 10);
+
+  for (i = 0; i < 32; ++i) {
+    u8 byte = n[31 - i];
+    for (j = 0; j < 8; ++j) {
+      const limb bit = byte >> 7;
+
+      swap_conditional(nqx, nqpqx, bit);
+      swap_conditional(nqz, nqpqz, bit);
+      fmonty(nqx2, nqz2,
+             nqpqx2, nqpqz2,
+             nqx, nqz,
+             nqpqx, nqpqz,
+             q);
+      swap_conditional(nqx2, nqpqx2, bit);
+      swap_conditional(nqz2, nqpqz2, bit);
+
+      t = nqx;
+      nqx = nqx2;
+      nqx2 = t;
+      t = nqz;
+      nqz = nqz2;
+      nqz2 = t;
+      t = nqpqx;
+      nqpqx = nqpqx2;
+      nqpqx2 = t;
+      t = nqpqz;
+      nqpqz = nqpqz2;
+      nqpqz2 = t;
+
+      byte <<= 1;
+    }
+  }
+
+  memcpy(resultx, nqx, sizeof(limb) * 10);
+  memcpy(resultz, nqz, sizeof(limb) * 10);
+}
+
+/* -----------------------------------------------------------------------------
+ * Shamelessly copied from djb's code
+ * ----------------------------------------------------------------------------- */
+static void
+crecip(limb *out, const limb *z) {
+  limb z2[10];
+  limb z9[10];
+  limb z11[10];
+  limb z2_5_0[10];
+  limb z2_10_0[10];
+  limb z2_20_0[10];
+  limb z2_50_0[10];
+  limb z2_100_0[10];
+  limb t0[10];
+  limb t1[10];
+  int i;
+
+  /* 2 */ fsquare(z2,z);
+  /* 4 */ fsquare(t1,z2);
+  /* 8 */ fsquare(t0,t1);
+  /* 9 */ fmul(z9,t0,z);
+  /* 11 */ fmul(z11,z9,z2);
+  /* 22 */ fsquare(t0,z11);
+  /* 2^5 - 2^0 = 31 */ fmul(z2_5_0,t0,z9);
+
+  /* 2^6 - 2^1 */ fsquare(t0,z2_5_0);
+  /* 2^7 - 2^2 */ fsquare(t1,t0);
+  /* 2^8 - 2^3 */ fsquare(t0,t1);
+  /* 2^9 - 2^4 */ fsquare(t1,t0);
+  /* 2^10 - 2^5 */ fsquare(t0,t1);
+  /* 2^10 - 2^0 */ fmul(z2_10_0,t0,z2_5_0);
+
+  /* 2^11 - 2^1 */ fsquare(t0,z2_10_0);
+  /* 2^12 - 2^2 */ fsquare(t1,t0);
+  /* 2^20 - 2^10 */ for (i = 2;i < 10;i += 2) { fsquare(t0,t1); fsquare(t1,t0); }
+  /* 2^20 - 2^0 */ fmul(z2_20_0,t1,z2_10_0);
+
+  /* 2^21 - 2^1 */ fsquare(t0,z2_20_0);
+  /* 2^22 - 2^2 */ fsquare(t1,t0);
+  /* 2^40 - 2^20 */ for (i = 2;i < 20;i += 2) { fsquare(t0,t1); fsquare(t1,t0); }
+  /* 2^40 - 2^0 */ fmul(t0,t1,z2_20_0);
+
+  /* 2^41 - 2^1 */ fsquare(t1,t0);
+  /* 2^42 - 2^2 */ fsquare(t0,t1);
+  /* 2^50 - 2^10 */ for (i = 2;i < 10;i += 2) { fsquare(t1,t0); fsquare(t0,t1); }
+  /* 2^50 - 2^0 */ fmul(z2_50_0,t0,z2_10_0);
+
+  /* 2^51 - 2^1 */ fsquare(t0,z2_50_0);
+  /* 2^52 - 2^2 */ fsquare(t1,t0);
+  /* 2^100 - 2^50 */ for (i = 2;i < 50;i += 2) { fsquare(t0,t1); fsquare(t1,t0); }
+  /* 2^100 - 2^0 */ fmul(z2_100_0,t1,z2_50_0);
+
+  /* 2^101 - 2^1 */ fsquare(t1,z2_100_0);
+  /* 2^102 - 2^2 */ fsquare(t0,t1);
+  /* 2^200 - 2^100 */ for (i = 2;i < 100;i += 2) { fsquare(t1,t0); fsquare(t0,t1); }
+  /* 2^200 - 2^0 */ fmul(t1,t0,z2_100_0);
+
+  /* 2^201 - 2^1 */ fsquare(t0,t1);
+  /* 2^202 - 2^2 */ fsquare(t1,t0);
+  /* 2^250 - 2^50 */ for (i = 2;i < 50;i += 2) { fsquare(t0,t1); fsquare(t1,t0); }
+  /* 2^250 - 2^0 */ fmul(t0,t1,z2_50_0);
+
+  /* 2^251 - 2^1 */ fsquare(t1,t0);
+  /* 2^252 - 2^2 */ fsquare(t0,t1);
+  /* 2^253 - 2^3 */ fsquare(t1,t0);
+  /* 2^254 - 2^4 */ fsquare(t0,t1);
+  /* 2^255 - 2^5 */ fsquare(t1,t0);
+  /* 2^255 - 21 */ fmul(out,t1,z11);
+}
+
+int curve25519_donna(u8 *, const u8 *, const u8 *);
+
+int curve25519_donna(u8 *mypublic, const u8 *secret, const u8 *basepoint) {
+  limb bp[10], x[10], z[11], zmone[10];
+  uint8_t e[32];
+  int i;
+
+  for (i = 0; i < 32; ++i) e[i] = secret[i];
+  e[0] &= 248;
+  e[31] &= 127;
+  e[31] |= 64;
+
+  fexpand(bp, basepoint);
+  cmult(x, z, e, bp);
+  crecip(zmone, z);
+  fmul(z, x, zmone);
+  freduce_coefficients(z);
+  fcontract(mypublic, z);
+  return 0;
+}
diff --git a/jni/libzrtp/sources/bnlib/ec/ec.c b/jni/libzrtp/sources/bnlib/ec/ec.c
index 1f4123a..18e612f 100644
--- a/jni/libzrtp/sources/bnlib/ec/ec.c
+++ b/jni/libzrtp/sources/bnlib/ec/ec.c
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2012 Werner Dittmann
+ * Copyright (C) 2012-2013 Werner Dittmann
  * All rights reserved. For licensing and other legal details, see the file legal.c.
  *
  * @author Werner Dittmann <Werner.Dittmann@t-online.de>
@@ -100,6 +100,46 @@
         "11839296a789a3bc0045c8a5fb42c7d1bd998f54449579b446817afbd17273e662c97ee72995ef42640c550b9013fad0761353c7086a272c24088be94769fd16650",
 };
 
+
+/*
+ * The data for curve3617 copied from:
+ * http://safecurves.cr.yp.to/field.html
+ * http://safecurves.cr.yp.to/base.html
+ */
+static curveData curve3617 = {
+    "3fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffef",  /* Prime */
+    "7ffffffffffffffffffffffffffffffffffffffffffffffffffeb3cc92414cf706022b36f1c0338ad63cf181b0e71a5e106af79",   /* order */
+    "",                                                                                                          /* SEED */
+    "",                                                                                                          /* c */
+    "",                                                                                                          /* b */
+    "1a334905141443300218c0631c326e5fcd46369f44c03ec7f57ff35498a4ab4d6d6ba111301a73faa8537c64c4fd3812f3cbc595",  /* Gx*/
+    "22",                                                                                                        /* Gy (radix 16) */
+};
+
+/*
+ * The data for curve25519 copied from:
+ * http://safecurves.cr.yp.to/field.html
+ * http://safecurves.cr.yp.to/base.html
+ * 
+ * Note: 
+ * The data for Curve25519 is here for the sake of completeness and to have the same
+ * set of initialization. One exception if the base point X coordinate (Gx) that we use to
+ * compute the DH public value, refer to function ecdhGeneratePublic(...) in ecdh.c.
+ * 
+ * Otherwise the functions use EcCurve structure only to get the pointers to the Curve25519
+ * wrapper functions.
+ * 
+ */
+static curveData curve25519 = {
+    "7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffed",   /* Prime */
+    "1000000000000000000000000000000014def9dea2f79cd65812631a5cf5d3ed",   /* order */
+    "",                                                                   /* SEED */
+    "",                                                                   /* c */
+    "",                                                                   /* b */
+    "9",                                                                  /* Gx */
+    "20ae19a1b8a086b4e01edd2c7748d14c923d4d7e6d7c61b229e9c5a27eced3d9",   /* Gy */
+};
+
 /*============================================================================*/
 /*    Bignum Shorthand Functions                                              */
 /*============================================================================*/
@@ -140,53 +180,90 @@
     return 0;
 }
 
-int bnMulMod_ (struct BigNum *rslt, struct BigNum *n1, struct BigNum *n2, struct BigNum *mod)
+int bnMulMod_ (struct BigNum *rslt, struct BigNum *n1, struct BigNum *n2, struct BigNum *mod, const EcCurve *curve)
 {
     bnMul (rslt, n1, n2);
-    bnMod (rslt, rslt, mod);
+    if (curve)
+        curve->modOp(rslt, rslt, mod);
+    else
+        bnMod(rslt, rslt, mod);
     return 0;
 }
 
-int bnMulQMod_ (struct BigNum *rslt, struct BigNum *n1, unsigned n2, struct BigNum *mod)
+int bnMulQMod_ (struct BigNum *rslt, struct BigNum *n1, unsigned n2, struct BigNum *mod, const EcCurve *curve)
 {
     bnMulQ (rslt, n1, n2);
-    bnMod (rslt, rslt, mod);
-    return 0;
+    if (curve)
+        curve->modOp(rslt, rslt, mod);
+    else
+        bnMod(rslt, rslt, mod);
+   return 0;
 }
 
-int bnSquareMod_ (struct BigNum *rslt, struct BigNum *n1, struct BigNum *mod)
+int bnSquareMod_ (struct BigNum *rslt, struct BigNum *n1, struct BigNum *mod, const EcCurve *curve)
 {
     bnSquare (rslt, n1);
-    bnMod (rslt, rslt, mod);
+    if (curve)
+        curve->modOp(rslt, rslt, mod);
+    else
+        bnMod(rslt, rslt, mod);
     return 0;
 }
 
-int ecGetCurveNistECp(NistCurves curveId, NistECpCurve *curve)
+/*
+ * Note on the Curve25519 functions and usage of BigNumber:
+ * In most cases the functions to compute Curve25519 data are small wrapper functions
+ * that implement the same API as for the other curve functions. The wrapper functions
+ * then call the very specific, high-efficient function in curve25519-donna.c .
+ * 
+ * For Curve25519 we don't have a real implementation for point add, point doubling, modulo
+ * and check public key. Please refer to the actual implementations below.
+ */
+
+static int ecGetAffineNist(const EcCurve *curve, EcPoint *R, const EcPoint *P);
+static int ecGetAffineEd(const EcCurve *curve, EcPoint *R, const EcPoint *P);
+static int ecGetAffine25519(const EcCurve *curve, EcPoint *R, const EcPoint *P);
+
+static int ecDoublePointNist(const EcCurve *curve, EcPoint *R, const EcPoint *P);
+static int ecDoublePointEd(const EcCurve *curve, EcPoint *R, const EcPoint *P);
+static int ecDoublePoint25519(const EcCurve *curve, EcPoint *R, const EcPoint *P);
+
+static int ecAddPointNist(const EcCurve *curve, EcPoint *R, const EcPoint *P, const EcPoint *Q);
+static int ecAddPointEd(const EcCurve *curve, EcPoint *R, const EcPoint *P, const EcPoint *Q);
+static int ecAddPoint25519(const EcCurve *curve, EcPoint *R, const EcPoint *P, const EcPoint *Q);
+
+static int ecCheckPubKeyNist(const EcCurve *curve, const EcPoint *pub);
+static int ecCheckPubKey3617(const EcCurve *curve, const EcPoint *pub);
+static int ecCheckPubKey25519(const EcCurve *curve, const EcPoint *pub);
+
+static int ecGenerateRandomNumberNist(const EcCurve *curve, BigNum *d);
+static int ecGenerateRandomNumber3617(const EcCurve *curve, BigNum *d);
+static int ecGenerateRandomNumber25519(const EcCurve *curve, BigNum *d);
+
+static int ecMulPointScalarNormal(const EcCurve *curve, EcPoint *R, const EcPoint *P, const BigNum *scalar);
+static int ecMulPointScalar25519(const EcCurve *curve, EcPoint *R, const EcPoint *P, const BigNum *scalar);
+
+/* Forward declaration of new modulo functions for the EC curves */
+static int newMod192(BigNum *r, const BigNum *a, const BigNum *modulo);
+static int newMod256(BigNum *r, const BigNum *a, const BigNum *modulo);
+static int newMod384(BigNum *r, const BigNum *a, const BigNum *modulo);
+static int newMod521(BigNum *r, const BigNum *a, const BigNum *modulo);
+
+static int mod3617(BigNum *r, const BigNum *a, const BigNum *modulo);
+static int mod25519(BigNum *r, const BigNum *a, const BigNum *modulo);
+
+static void commonInit()
 {
-    size_t maxBits;
-    curveData *cd;
+    bnBegin(mpiZero); bnSetQ(mpiZero, 0);
+    bnBegin(mpiOne); bnSetQ(mpiOne, 1);
+    bnBegin(mpiTwo); bnSetQ(mpiTwo, 2);
+    bnBegin(mpiThree); bnSetQ(mpiThree, 3);
+    bnBegin(mpiFour); bnSetQ(mpiFour, 4);
+    bnBegin(mpiEight); bnSetQ(mpiEight, 8);
+}
 
-    if (!initialized) {
-        bnBegin(mpiZero); bnSetQ(mpiZero, 0);
-        bnBegin(mpiOne); bnSetQ(mpiOne, 1);
-        bnBegin(mpiTwo); bnSetQ(mpiTwo, 2);
-        bnBegin(mpiThree); bnSetQ(mpiThree, 3);
-        bnBegin(mpiFour); bnSetQ(mpiFour, 4);
-        bnBegin(mpiEight); bnSetQ(mpiEight, 8);
-        initialized = 1;
-    }
-    if (curve == NULL)
-        return -2;
-
-    bnBegin(&curve->_p);    curve->p = &curve->_p;
-    bnBegin(&curve->_n);    curve->n = &curve->_n;
-    bnBegin(&curve->_SEED); curve->SEED = &curve->_SEED;
-    bnBegin(&curve->_c);    curve->c = &curve->_c;
-    bnBegin(&curve->_a);    curve->a = &curve->_a;
-    bnBegin(&curve->_b);    curve->b = &curve->_b;
-    bnBegin(&curve->_Gx);   curve->Gx = &curve->_Gx;
-    bnBegin(&curve->_Gy);   curve->Gy = &curve->_Gy;
-
+static void curveCommonInit(EcCurve *curve)
+{
     /* Initialize scratchpad variables and their pointers */
     bnBegin(&curve->_S1); curve->S1 = &curve->_S1;
     bnBegin(&curve->_U1); curve->U1 = &curve->_U1;
@@ -196,41 +273,11 @@
     bnBegin(&curve->_t1); curve->t1 = &curve->_t1;
     bnBegin(&curve->_t2); curve->t2 = &curve->_t2;
     bnBegin(&curve->_t3); curve->t3 = &curve->_t3;
+}
 
-    switch (curveId) {
-    case NIST192P:
-        cd = &nist192;
-        break;
-
-    case NIST224P:
-        cd = &nist224;
-        break;
-
-    case NIST256P:
-        cd = &nist256;
-        break;
-
-    case NIST384P:
-        cd = &nist384;
-        break;
-
-    case NIST521P:
-        cd = &nist521;
-        break;
-
-    default:
-        return -2;
-    }
-
-    bnReadAscii(curve->p, cd->p, 10);
-    bnReadAscii(curve->n, cd->n, 10);
-    bnReadAscii(curve->SEED, cd->SEED, 16);
-    bnReadAscii(curve->c, cd->c, 16);
-    bnCopy(curve->a, curve->p);
-    bnSub(curve->a, mpiThree);
-    bnReadAscii(curve->b, cd->b, 16);
-    bnReadAscii(curve->Gx, cd->Gx, 16);
-    bnReadAscii(curve->Gy, cd->Gy, 16);
+static void curveCommonPrealloc(EcCurve *curve)
+{
+    size_t maxBits;
 
     /* variables must be able to hold p^2, plus one nimb (min. 15 bits) for overflow */
     maxBits = bnBits(curve->p) * 2 + 15;
@@ -245,16 +292,151 @@
     bnPrealloc(curve->t1, maxBits);
     bnPrealloc(curve->t2, maxBits);
     bnPrealloc(curve->t3, maxBits);
-
-    return 0;
-
-/*     ecFreeCurveNistECp(curve);
-     return ret;
-*/
 }
 
+int ecGetCurveNistECp(Curves curveId, EcCurve *curve)
+{
+    curveData *cd;
 
-void ecFreeCurveNistECp(NistECpCurve *curve) 
+    if (curveId >= Curve25519 && curveId <= Curve3617)
+        return ecGetCurvesCurve(curveId, curve);
+
+    if (!initialized) {
+        commonInit();
+        initialized = 1;
+    }
+    if (curve == NULL)
+        return -2;
+
+    bnBegin(&curve->_p);    curve->p = &curve->_p;
+    bnBegin(&curve->_n);    curve->n = &curve->_n;
+    bnBegin(&curve->_SEED); curve->SEED = &curve->_SEED;
+    bnBegin(&curve->_c);    curve->c = &curve->_c;
+    bnBegin(&curve->_a);    curve->a = &curve->_a;
+    bnBegin(&curve->_b);    curve->b = &curve->_b;
+    bnBegin(&curve->_Gx);   curve->Gx = &curve->_Gx;
+    bnBegin(&curve->_Gy);   curve->Gy = &curve->_Gy;
+
+    curveCommonInit(curve);
+
+    switch (curveId) {
+    case NIST192P:
+        cd = &nist192;
+        curve->modOp = newMod192;
+        break;
+
+    case NIST224P:
+        cd = &nist224;
+        curve->modOp = bnMod;
+        break;
+
+    case NIST256P:
+        cd = &nist256;
+        curve->modOp = bnMod;
+        break;
+
+    case NIST384P:
+        cd = &nist384;
+        curve->modOp = newMod384;
+        break;
+
+    case NIST521P:
+        cd = &nist521;
+        curve->modOp = newMod521;
+        break;
+
+    default:
+        return -2;
+    }
+
+    curve->affineOp = ecGetAffineNist;
+    curve->doubleOp = ecDoublePointNist;
+    curve->addOp = ecAddPointNist;
+    curve->checkPubOp = ecCheckPubKeyNist;
+    curve->randomOp = ecGenerateRandomNumberNist;
+    curve->mulScalar = ecMulPointScalarNormal;
+
+    bnReadAscii(curve->p, cd->p, 10);
+    bnReadAscii(curve->n, cd->n, 10);
+    bnReadAscii(curve->SEED, cd->SEED, 16);
+    bnReadAscii(curve->c, cd->c, 16);
+    bnCopy(curve->a, curve->p);
+    bnSub(curve->a, mpiThree);
+    bnReadAscii(curve->b, cd->b, 16);
+    bnReadAscii(curve->Gx, cd->Gx, 16);
+    bnReadAscii(curve->Gy, cd->Gy, 16);
+
+    curveCommonPrealloc(curve);
+    curve->id = curveId;
+
+    return 0;
+}
+
+int ecGetCurvesCurve(Curves curveId, EcCurve *curve)
+{
+    curveData *cd;
+
+    if (!initialized) {
+        commonInit();
+        initialized = 1;
+    }
+    if (curve == NULL)
+        return -2;
+
+    /* set-up all bignum structures, simplifies "free" handling */
+    bnBegin(&curve->_p);    curve->p = &curve->_p;
+    bnBegin(&curve->_n);    curve->n = &curve->_n;
+    bnBegin(&curve->_SEED); curve->SEED = &curve->_SEED;
+    bnBegin(&curve->_c);    curve->c = &curve->_c;
+    bnBegin(&curve->_a);    curve->a = &curve->_a;
+    bnBegin(&curve->_b);    curve->b = &curve->_b;
+    bnBegin(&curve->_Gx);   curve->Gx = &curve->_Gx;
+    bnBegin(&curve->_Gy);   curve->Gy = &curve->_Gy;
+
+    curveCommonInit(curve);
+
+    switch (curveId) {
+    case Curve3617:
+        cd = &curve3617;
+        curve->modOp = mod3617;
+        curve->affineOp = ecGetAffineEd;
+        curve->doubleOp = ecDoublePointEd;
+        curve->addOp = ecAddPointEd;
+        curve->checkPubOp = ecCheckPubKey3617;
+        curve->randomOp = ecGenerateRandomNumber3617;
+        curve->mulScalar = ecMulPointScalarNormal;
+
+        bnReadAscii(curve->a, "3617", 10);
+        break;
+
+    case Curve25519:
+        cd = &curve25519;
+        curve->modOp = mod25519;
+        curve->affineOp = ecGetAffine25519;
+        curve->doubleOp = ecDoublePoint25519;
+        curve->addOp = ecAddPoint25519;
+        curve->checkPubOp = ecCheckPubKey25519;
+        curve->randomOp = ecGenerateRandomNumber25519;
+        curve->mulScalar = ecMulPointScalar25519;
+
+        bnReadAscii(curve->a, "486662", 10);
+        break;
+
+    default:
+        return -2;
+    }
+    bnReadAscii(curve->p, cd->p, 16);
+    bnReadAscii(curve->n, cd->n, 16);
+
+    bnReadAscii(curve->Gx, cd->Gx, 16);
+    bnReadAscii(curve->Gy, cd->Gy, 16);
+
+    curveCommonPrealloc(curve);
+    curve->id = curveId;
+    return 0;
+}
+
+void ecFreeCurveNistECp(EcCurve *curve) 
 {
     if (curve == NULL)
         return;
@@ -277,12 +459,40 @@
     bnEnd(curve->t3);
 }
 
+/*
+ * EC point helper functions
+ */
+
+void ecInitPoint(EcPoint *P)
+{
+    INIT_EC_POINT(P);
+}
+
+void ecFreePoint(EcPoint *P)
+{
+    FREE_EC_POINT(P);
+}
+
+void ecSetBasePoint(EcCurve *C, EcPoint *P)
+{
+    SET_EC_BASE_POINT(C, P);
+}
+
+void ecFreeCurvesCurve(EcCurve *curve)
+{
+    ecFreeCurveNistECp(curve);
+}
 
 /*============================================================================*/
 /*    Elliptic Curve arithmetic                                               */
 /*============================================================================*/
 
-int ecGetAffine(const NistECpCurve *curve, EcPoint *R, const EcPoint *P)
+int ecGetAffine(const EcCurve *curve, EcPoint *R, const EcPoint *P)
+{
+    return curve->affineOp(curve, R, P);
+}
+
+static int ecGetAffineNist(const EcCurve *curve, EcPoint *R, const EcPoint *P)
 {
     int ret = 0;
 
@@ -292,13 +502,13 @@
     bnBegin(&z_2);
 
     /* affine x = X / Z^2 */
-    bnInv (&z_1, P->z, curve->p);          /* z_1 = Z^(-1) */
-    bnMulMod_(&z_2, &z_1, &z_1, curve->p); /* z_2 = Z^(-2) */
-    bnMulMod_(R->x, P->x, &z_2, curve->p);
-    
+    bnInv (&z_1, P->z, curve->p);                 /* z_1 = Z^(-1) */
+    bnMulMod_(&z_2, &z_1, &z_1, curve->p, curve); /* z_2 = Z^(-2) */
+    bnMulMod_(R->x, P->x, &z_2, curve->p, curve);
+
     /* affine y = Y / Z^3 */
-    bnMulMod_(&z_2, &z_2, &z_1, curve->p); /* z_2 = Z^(-3) */
-    bnMulMod_(R->y, P->y, &z_2, curve->p);
+    bnMulMod_(&z_2, &z_2, &z_1, curve->p, curve); /* z_2 = Z^(-3) */
+    bnMulMod_(R->y, P->y, &z_2, curve->p, curve);
 
     bnSetQ(R->z, 1);
 
@@ -307,7 +517,48 @@
     return ret;
 }
 
-int ecDoublePoint(const NistECpCurve *curve, EcPoint *R, const EcPoint *P)
+static int ecGetAffineEd(const EcCurve *curve, EcPoint *R, const EcPoint *P)
+{
+    int ret = 0;
+
+    struct BigNum z_1;
+
+    bnBegin(&z_1);
+
+    /* affine x = X / Z */
+    bnInv (&z_1, P->z, curve->p);                 /* z_1 = Z^(-1) */
+    bnMulMod_(R->x, P->x, &z_1, curve->p, curve);
+
+    /* affine y = Y / Z */
+    bnMulMod_(R->y, P->y, &z_1, curve->p, curve);
+
+    bnSetQ(R->z, 1);
+
+    bnEnd(&z_1);
+    return ret;
+
+}
+
+/* 
+ * If the arguments do not point to the same EcPoint then copy P to result.
+ * Curve25519 has no specific GetAffine function, it's all inside curve25519-donna
+ */
+static int ecGetAffine25519(const EcCurve *curve, EcPoint *R, const EcPoint *P)
+{
+    if (R != P) {
+        bnCopy(R->x, P->x);
+        bnCopy(R->y, P->y);
+        bnCopy(R->z, P->z);
+    }
+    return 0;
+}
+
+int ecDoublePoint(const EcCurve *curve, EcPoint *R, const EcPoint *P)
+{
+    return curve->doubleOp(curve, R, P);
+}
+
+static int ecDoublePointNist(const EcCurve *curve, EcPoint *R, const EcPoint *P)
 {
     int ret = 0;
 
@@ -333,37 +584,37 @@
         ptP = P;
 
     /* S = 4*X*Y^2, save Y^2 in t1 for later use */
-    bnMulMod_(curve->t1, ptP->y, ptP->y, curve->p);       /* t1 = Y^2 */
-    bnMulMod_(curve->t0, ptP->x, mpiFour, curve->p);      /* t0 = 4 * X */
-    bnMulMod_(curve->S1, curve->t0, curve->t1, curve->p); /* S1 = t0 * t1 */
+    bnMulMod_(curve->t1, ptP->y, ptP->y, curve->p, curve);       /* t1 = Y^2 */
+    bnMulMod_(curve->t0, ptP->x, mpiFour, curve->p, curve);      /* t0 = 4 * X */
+    bnMulMod_(curve->S1, curve->t0, curve->t1, curve->p, curve); /* S1 = t0 * t1 */
 
     /* M = 3*(X + Z^2)*(X - Z^2), use scratch variable U1 to store M value */
-    bnMulMod_(curve->t2, ptP->z, ptP->z, curve->p);       /* t2 = Z^2 */
+    bnMulMod_(curve->t2, ptP->z, ptP->z, curve->p, curve);       /* t2 = Z^2 */
     bnCopy(curve->t0, ptP->x);
-    bnAddMod_(curve->t0, curve->t2, curve->p);            /* t0 = X + t2  */
-    bnMulMod_(curve->t3, curve->t0, mpiThree, curve->p);  /* t3 = 3 * t0 */
+    bnAddMod_(curve->t0, curve->t2, curve->p);                   /* t0 = X + t2  */
+    bnMulMod_(curve->t3, curve->t0, mpiThree, curve->p, curve);  /* t3 = 3 * t0 */
     bnCopy(curve->t0, ptP->x);
-    bnSubMod_(curve->t0, curve->t2, curve->p);            /* t0 = X - t2 */
-    bnMulMod_(curve->U1, curve->t3, curve->t0, curve->p); /* M = t3 * t0 */
+    bnSubMod_(curve->t0, curve->t2, curve->p);                   /* t0 = X - t2 */
+    bnMulMod_(curve->U1, curve->t3, curve->t0, curve->p, curve); /* M = t3 * t0 */
     
     /* X' = M^2 - 2*S */
-    bnMulMod_(curve->t2, curve->U1, curve->U1, curve->p); /* t2 = M^2 */
-    bnMulMod_(curve->t0, curve->S1, mpiTwo, curve->p);    /* t0 = S * 2 */
+    bnMulMod_(curve->t2, curve->U1, curve->U1, curve->p, curve); /* t2 = M^2 */
+    bnMulMod_(curve->t0, curve->S1, mpiTwo, curve->p, curve);    /* t0 = S * 2 */
     bnCopy(R->x, curve->t2);
-    bnSubMod_(R->x, curve->t0, curve->p);                 /* X' = t2 - t0 */
+    bnSubMod_(R->x, curve->t0, curve->p);                        /* X' = t2 - t0 */
 
     /* Y' = M*(S - X') - 8*Y^4 */
-    bnMulMod_(curve->t3, curve->t1, curve->t1, curve->p); /* t3 = Y^4 (t1 saved above) */
-    bnMulMod_(curve->t2, curve->t3, mpiEight, curve->p); /* t2 = t3 * 8 */
+    bnMulMod_(curve->t3, curve->t1, curve->t1, curve->p, curve); /* t3 = Y^4 (t1 saved above) */
+    bnMulMod_(curve->t2, curve->t3, mpiEight, curve->p, curve);  /* t2 = t3 * 8 */
     bnCopy(curve->t3, curve->S1);
-    bnSubMod_(curve->t3, R->x, curve->p);                 /* t3 = S - X' */
-    bnMulMod_(curve->t0, curve->U1, curve->t3, curve->p); /* t0 = M * t3 */
+    bnSubMod_(curve->t3, R->x, curve->p);                        /* t3 = S - X' */
+    bnMulMod_(curve->t0, curve->U1, curve->t3, curve->p, curve); /* t0 = M * t3 */
     bnCopy(R->y, curve->t0);
-    bnSubMod_(R->y, curve->t2, curve->p);                 /* Y' = t0 - t2 */
+    bnSubMod_(R->y, curve->t2, curve->p);                        /* Y' = t0 - t2 */
 
     /* Z' = 2*Y*Z */
-    bnMulMod_(curve->t0, ptP->y, mpiTwo, curve->p);       /* t0 = 2 * Y */
-    bnMulMod_(R->z, curve->t0, ptP->z, curve->p);         /* Z' = to * Z */
+    bnMulMod_(curve->t0, ptP->y, mpiTwo, curve->p, curve);       /* t0 = 2 * Y */
+    bnMulMod_(R->z, curve->t0, ptP->z, curve->p, curve);         /* Z' = to * Z */
 
     if (P == R)
         FREE_EC_POINT(&tP);
@@ -371,8 +622,72 @@
     return ret;
 }
 
+static int ecDoublePointEd(const EcCurve *curve, EcPoint *R, const EcPoint *P)
+{
+    EcPoint tP;
+    const EcPoint *ptP = 0;
+
+    /* Check for overlapping arguments, copy if necessary and set pointer */
+    if (P == R) {
+        INIT_EC_POINT(&tP);
+        ptP = &tP;
+        bnCopy(tP.x, P->x);
+        bnCopy(tP.y, P->y);
+        bnCopy(tP.z, P->z);
+    }
+    else 
+        ptP = P;
+
+    /* Compute B, C, D, H, E */
+    bnCopy(curve->t1, ptP->x);
+    bnAddMod_(curve->t1, ptP->y, curve->p);
+    bnSquareMod_(curve->t0, curve->t1, curve->p, curve);       /* t0 -> B */
+
+    bnSquareMod_(R->x, ptP->x, curve->p, curve);               /* Rx -> C */
+
+    bnSquareMod_(R->y, ptP->y, curve->p, curve);               /* Ry -> D */
+
+    bnSquareMod_(R->z, ptP->z, curve->p, curve);               /* Rz -> H */
+    bnAddMod_(R->z, R->z, curve->p);                           /* Rz -> 2H */
+
+    bnCopy(curve->t1, R->x);
+    bnAddMod_(curve->t1, R->y, curve->p);                      /* t1 -> E */
+
+    /* Compute Ry */
+    bnCopy(curve->t2, R->x);
+    bnSubMod_(curve->t2, R->y, curve->p);                      /* C - D */
+    bnMulMod_(R->y, curve->t1, curve->t2, curve->p, curve);    /* E * t3; Ry */
+
+    /* Compute Rx */
+    bnSubMod_(curve->t0, curve->t1, curve->p);                 /* B - E; sub result */
+    bnCopy(curve->t2, curve->t1);
+    bnSubMod_(curve->t2, R->z, curve->p);                      /* t2 -> J; (E - 2H) */
+    bnMulMod_(R->x, curve->t2, curve->t0, curve->p, curve);    /* J * t0 */
+
+    /* Compute Rz */
+    bnMulMod_(R->z, curve->t2, curve->t1, curve->p, curve);    /* J * E */
+
+    if (P == R)
+        FREE_EC_POINT(&tP);
+
+    return 0;
+}
+
+/* 
+ * Curve25519 has no specific Double Point function, all inside curve25519-donna
+ */
+static int ecDoublePoint25519(const EcCurve *curve, EcPoint *R, const EcPoint *P)
+{
+    return -2;
+}
+
 /* Add two elliptic curve points. Any of them may be the same object. */
-int ecAddPoint(const NistECpCurve *curve, EcPoint *R, const EcPoint *P, const EcPoint *Q)
+int ecAddPoint(const EcCurve *curve, EcPoint *R, const EcPoint *P, const EcPoint *Q)
+{
+    return curve->addOp(curve, R, P, Q);
+}
+
+static int ecAddPointNist(const EcCurve *curve, EcPoint *R, const EcPoint *P, const EcPoint *Q)
 {
     int ret = 0;
 
@@ -424,23 +739,23 @@
         ptQ = Q;
 
     /* U1 = X1*Z2^2, where X1: P->x, Z2: Q->z */
-    bnMulMod_(curve->t1, ptQ->z, ptQ->z, curve->p);    /* t1 = Z2^2 */
-    bnMulMod_(curve->U1, ptP->x, curve->t1, curve->p); /* U1 = X1 * z_2 */
+    bnMulMod_(curve->t1, ptQ->z, ptQ->z, curve->p, curve);    /* t1 = Z2^2 */
+    bnMulMod_(curve->U1, ptP->x, curve->t1, curve->p, curve); /* U1 = X1 * z_2 */
 
     /* S1 = Y1*Z2^3, where Y1: P->y */
-    bnMulMod_(curve->t1, curve->t1, ptQ->z, curve->p); /* t1 = Z2^3 */
-    bnMulMod_(curve->S1, ptP->y, curve->t1, curve->p); /* S1 = Y1 * z_2 */
+    bnMulMod_(curve->t1, curve->t1, ptQ->z, curve->p, curve); /* t1 = Z2^3 */
+    bnMulMod_(curve->S1, ptP->y, curve->t1, curve->p, curve); /* S1 = Y1 * z_2 */
 
     /* U2 = X2*Z1^2, where X2: Q->x, Z1: P->z */
-    bnMulMod_(curve->t1, ptP->z, ptP->z, curve->p);    /* t1 = Z1^2 */
-    bnMulMod_(curve->H, ptQ->x, curve->t1, curve->p);  /* H = X2 * t1 (store U2 in H) */
+    bnMulMod_(curve->t1, ptP->z, ptP->z, curve->p, curve);    /* t1 = Z1^2 */
+    bnMulMod_(curve->H, ptQ->x, curve->t1, curve->p, curve);  /* H = X2 * t1 (store U2 in H) */
 
     /* H = U2 - U1 */
     bnSubMod_(curve->H, curve->U1, curve->p);
 
     /* S2 = Y2*Z1^3, where Y2: Q->y */
-    bnMulMod_(curve->t1, curve->t1, ptP->z, curve->p); /* t1 = Z1^3 */
-    bnMulMod_(curve->R, ptQ->y, curve->t1, curve->p);  /* R = Y2 * t1 (store S2 in R) */
+    bnMulMod_(curve->t1, curve->t1, ptP->z, curve->p, curve); /* t1 = Z1^3 */
+    bnMulMod_(curve->R, ptQ->y, curve->t1, curve->p, curve);  /* R = Y2 * t1 (store S2 in R) */
 
     /* R = S2 - S1 */
     bnSubMod_(curve->R, curve->S1, curve->p);
@@ -458,26 +773,26 @@
         return ecDoublePoint(curve, R, P);
     }
     /* X3 = R^2 - H^3 - 2*U1*H^2, where X3: R->x */
-    bnMulMod_(curve->t0, curve->H, curve->H, curve->p);   /* t0 = H^2 */
-    bnMulMod_(curve->t1, curve->U1, curve->t0, curve->p); /* t1 = U1 * t0, (hold t1) */
-    bnMulMod_(curve->t0, curve->t0, curve->H, curve->p);  /* t0 = H^3, (hold t0) */
-    bnMulMod_(curve->t2, curve->R, curve->R, curve->p);   /* t2 = R^2 */
+    bnMulMod_(curve->t0, curve->H, curve->H, curve->p, curve);   /* t0 = H^2 */
+    bnMulMod_(curve->t1, curve->U1, curve->t0, curve->p, curve); /* t1 = U1 * t0, (hold t1) */
+    bnMulMod_(curve->t0, curve->t0, curve->H, curve->p, curve);  /* t0 = H^3, (hold t0) */
+    bnMulMod_(curve->t2, curve->R, curve->R, curve->p, curve);   /* t2 = R^2 */
     bnCopy(curve->t3, curve->t2);
-    bnSubMod_(curve->t3, curve->t0, curve->p);            /* t3 = t2 - t0, (-H^3)*/
-    bnMulMod_(curve->t2, mpiTwo, curve->t1, curve->p);    /* t2 = 2 * t1 */
+    bnSubMod_(curve->t3, curve->t0, curve->p);                   /* t3 = t2 - t0, (-H^3)*/
+    bnMulMod_(curve->t2, mpiTwo, curve->t1, curve->p, curve);    /* t2 = 2 * t1 */
     bnCopy(R->x, curve->t3);
-    bnSubMod_(R->x, curve->t2, curve->p);                 /* X3 = t3 - t2 */
+    bnSubMod_(R->x, curve->t2, curve->p);                        /* X3 = t3 - t2 */
 
     /* Y3 = R*(U1*H^2 - X3) - S1*H^3, where Y3: R->y */
-    bnSubMod_(curve->t1, R->x, curve->p);                 /* t1 = t1 - X3, overwrites t1 now */
-    bnMulMod_(curve->t2, curve->R, curve->t1, curve->p);  /* t2 = R * z_2 */
-    bnMulMod_(curve->S1, curve->S1, curve->t0, curve->p); /* S1 = S1 * t0, (t0 has H^3) */
+    bnSubMod_(curve->t1, R->x, curve->p);                        /* t1 = t1 - X3, overwrites t1 now */
+    bnMulMod_(curve->t2, curve->R, curve->t1, curve->p, curve);  /* t2 = R * z_2 */
+    bnMulMod_(curve->S1, curve->S1, curve->t0, curve->p, curve); /* S1 = S1 * t0, (t0 has H^3) */
     bnCopy(R->y, curve->t2);
-    bnSubMod_(R->y, curve->S1, curve->p);                 /* Y3 = t2 - S1 */
+    bnSubMod_(R->y, curve->S1, curve->p);                        /* Y3 = t2 - S1 */
 
     /* Z3 = H*Z1*Z2, where Z1: P->z, Z2: Q->z, Z3: R->z */
-    bnMulMod_(curve->t2, curve->H, P->z, curve->p);       /* t2 = H * Z1 */
-    bnMulMod_(R->z, curve->t2, Q->z, curve->p);           /* Z3 = t2 * Z2 */
+    bnMulMod_(curve->t2, curve->H, P->z, curve->p, curve);       /* t2 = H * Z1 */
+    bnMulMod_(R->z, curve->t2, Q->z, curve->p, curve);           /* Z3 = t2 * Z2 */
 
     if (P == R)
         FREE_EC_POINT(&tP);
@@ -486,10 +801,117 @@
     return ret;
 }
 
-int ecMulPointScalar(const NistECpCurve *curve, EcPoint *R, const EcPoint *P, const BigNum *scalar)
+/*
+ * Refer to the document: Faster addition and doubling on elliptic curves; Daniel J. Bernstein and Tanja Lange
+ * section 4.
+ *
+ * This function is a variant of the 'addition'. The function returns the result in an own curve point
+ * and does not overwrite its input parameters.
+ */
+static int ecAddPointEd(const EcCurve *curve, EcPoint *R, const EcPoint *P, const EcPoint *Q)
 {
+    EcPoint tP, tQ;
+    const EcPoint *ptP = 0;
+    const EcPoint *ptQ = 0;
 
-    /* MPI_CHK below macro requires a 'ret' variable and a cleanup label */
+    /* if P is (@,@), R = Q */
+    if (!bnCmp(P->z, mpiZero)) {
+        bnCopy(R->x, Q->x);
+        bnCopy(R->y, Q->y);
+        bnCopy(R->z, Q->z);
+        return 0;
+    }
+
+    /* if Q is (@,@), R = P */
+    if (!bnCmp(Q->z, mpiZero)) {
+        bnCopy(R->x, P->x);
+        bnCopy(R->y, P->y);
+        bnCopy(R->z, P->z);
+        return 0;
+    }
+
+    /* Check for overlapping arguments, copy if necessary and set pointers */
+    if (P == R) {
+        INIT_EC_POINT(&tP);
+        ptP = &tP;
+        bnCopy(tP.x, P->x);
+        bnCopy(tP.y, P->y);
+        bnCopy(tP.z, P->z);
+    }
+    else 
+        ptP = P;
+
+    if (Q == R) {
+        INIT_EC_POINT(&tQ);
+        ptQ = &tQ;
+        bnCopy(tQ.x, Q->x);
+        bnCopy(tQ.y, Q->y);
+        bnCopy(tQ.z, Q->z);
+    }
+    else
+        ptQ = Q;
+
+    /* Compute A, C, D first */
+    bnMulMod_(R->z, ptP->z, ptQ->z, curve->p, curve);            /* Rz -> A; (Z1 * z2); Rz becomes R3 */
+    bnMulMod_(R->x, ptP->x, ptQ->x, curve->p, curve);            /* Rx -> C; (X1 * X2); Rx becomes R1 */
+    bnMulMod_(R->y, ptP->y, ptQ->y, curve->p, curve);            /* Ry -> D; (Y1 * Y2); Ry becomes R2 */
+
+    /* Compute large parts of X3 equation, sub result in t0 */
+    bnCopy(curve->t0, ptP->x);
+    bnAddMod_(curve->t0, ptP->y, curve->p);                      /* t0 -> X1 + Y1 */
+    bnCopy(curve->t1, ptQ->x);
+    bnAddMod_(curve->t1, ptQ->y, curve->p);                      /* t1 -> X2 + Y2 */
+    bnMulMod_(curve->t2, curve->t0, curve->t1, curve->p, curve); /* t2 = t0 * t1 */
+    bnSubMod_(curve->t2, R->x, curve->p);                        /* t2 - C */
+    bnSubMod_(curve->t2, R->y, curve->p);                        /* t2 - D */
+    bnMulMod_(curve->t0, curve->t2, R->z, curve->p, curve);      /* t0 -> R7; (t2 * A); sub result */
+
+    /* Compute E */
+    bnMulMod_(curve->t2, R->x, R->y, curve->p, curve);           /* t2 = C * D */
+    bnMulMod_(curve->t1, curve->t2, curve->a, curve->p, curve);  /* t1 -> E; t1 new R8 */
+
+    /* Compute part of Y3 equation, sub result in t2 */
+    bnSubMod_(R->y, R->x, curve->p);                             /* Ry = D - C; sub result */
+    bnMulMod_(curve->t2, R->y, R->z, curve->p, curve);           /* t2 = Ry * A; sub result */
+
+    /* Compute B */
+    bnSquareMod_(R->z, R->z, curve->p, curve);                   /* Rz -> B; (A^2) */
+
+    /* Compute F */
+    bnCopy(curve->t3, R->z);
+    bnSubMod_(curve->t3, curve->t1, curve->p);                   /* t3 -> F; (B - E) */
+
+    /* Compute G */
+    bnAddMod_(R->z, curve->t1, curve->p);                        /* Rz -> G; (B + E) */
+
+    /* Compute, X, Y, Z results */
+    bnMulMod_(R->x, curve->t3, curve->t0, curve->p, curve);      /* Rx = F * t0 */
+    bnMulMod_(R->y, curve->t2, R->z, curve->p, curve);           /* Ry = t2 * G */
+    bnMulMod_(R->z, curve->t3, R->z, curve->p, curve);           /* Rz = F * G */
+
+    if (P == R)
+        FREE_EC_POINT(&tP);
+    if (Q == R)
+        FREE_EC_POINT(&tQ);
+
+    return 0;
+}
+
+/* 
+ * Curve25519 has no specific Add Point function, all inside curve25519-donna
+ */
+static int ecAddPoint25519(const EcCurve *curve, EcPoint *R, const EcPoint *P, const EcPoint *Q)
+{
+    return -2;
+}
+
+int ecMulPointScalar(const EcCurve *curve, EcPoint *R, const EcPoint *P, const BigNum *scalar)
+{
+    return curve->mulScalar(curve, R, P, scalar);
+}
+
+static int ecMulPointScalarNormal(const EcCurve *curve, EcPoint *R, const EcPoint *P, const BigNum *scalar)
+{
     int ret = 0;
     int i;
     int bits = bnBits(scalar);
@@ -515,7 +937,26 @@
     return ret;
 }
 
+/* 
+ * This function uses BigNumber only as containers to transport the 32 byte data.
+ * This makes it compliant to the other functions and thus higher-level API does not change.
+ * 
+ * curve25519_donna function uses data in little endian format.
+ */
+static int ecMulPointScalar25519(const EcCurve *curve, EcPoint *R, const EcPoint *P, const BigNum *scalar)
+{
+    uint8_t basepoint[32], secret[32], result[32];
+
+    bnExtractLittleBytes(P->x, basepoint, 0, 32);  /* 25519 function requires the X coordinate only (compressed) */
+    bnExtractLittleBytes(scalar, secret, 0, 32);
+    curve25519_donna(result, secret, basepoint);
+    bnInsertLittleBytes(R->x, result, 0, 32);
+    return 0;
+}
+
 #ifdef WEAKRANDOM
+#include <fcntl.h>
+
 /*
  * A standard random number generator that uses the portable random() system function.
  *
@@ -523,10 +964,15 @@
  */
 static int _random(unsigned char *output, size_t len)
 {
-    size_t i;
+    size_t num = 0;
 
-    for(i = 0; i < len; ++i )
-        output[i] = random();
+    int rnd = open("/dev/urandom", O_RDONLY);
+    if (rnd >= 0) {
+        num = read(rnd, output, len);
+        close(rnd);
+    }
+    else
+        return num;
 
     return( 0 );
 }
@@ -538,7 +984,12 @@
 }
 #endif
 
-int ecGenerateRandomNumber(const NistECpCurve *curve, BigNum *d)
+int ecGenerateRandomNumber(const EcCurve *curve, BigNum *d)
+{
+    return curve->randomOp(curve, d);
+}
+
+static int ecGenerateRandomNumberNist(const EcCurve *curve, BigNum *d)
 {
     BigNum c, nMinusOne;
 
@@ -568,3 +1019,677 @@
 
     return 0;
 }
+
+static int ecGenerateRandomNumber3617(const EcCurve *curve, BigNum *d)
+{
+    unsigned char random[52];
+    _random(random, 52);
+
+    /* prepare the secret random data: clear bottom 3 bits. Clearing top 2 bits
+     * makes is a 414 bit value
+     */
+    random[51] &= ~0x7;
+    random[0] &= 0x3f;
+    /* convert the random data into big numbers */
+    bnInsertBigBytes(d, random, 0, 52);
+    return 0;
+}
+
+static int ecGenerateRandomNumber25519(const EcCurve *curve, BigNum *d)
+{
+    unsigned char random[32];
+    _random(random, 32);
+
+    /* No specific preparation. The curve25519_donna functions prepares the data.
+     *
+     * convert the random data into big numbers. the bigNumber is a container only.
+     * we don not use the big number for any arithmetic
+     */
+    bnInsertLittleBytes(d, random, 0, 32);
+    return 0;
+
+}
+
+int ecCheckPubKey(const EcCurve *curve, const EcPoint *pub)
+{
+    return curve->checkPubOp(curve, pub);
+}
+
+static int ecCheckPubKeyNist(const NistECpCurve *curve, const EcPoint *pub)
+{
+    /* Represent point at infinity by (0, 0), make sure it's not that */
+    if (bnCmpQ(pub->x, 0) == 0 && bnCmpQ(pub->y, 0) == 0) {
+        return 0;
+    }
+    /* Check that coordinates are within range */
+    if (bnCmpQ(pub->x, 0) < 0 || bnCmp(pub->x, curve->p) >= 0) {
+        return 0;
+    }
+    if (bnCmpQ(pub->y, 0) < 0 || bnCmp(pub->y, curve->p) >= 0) {
+        return 0;
+    }
+    /* Check that point satisfies EC equation y^2 = x^3 - 3x + b, mod P */
+    bnSquareMod_(curve->t1, pub->y, curve->p, curve);
+    bnSquareMod_(curve->t2, pub->x, curve->p, curve);
+    bnSubQMod_(curve->t2, 3, curve->p);
+    bnMulMod_(curve->t2, curve->t2, pub->x, curve->p, curve);
+    bnAddMod_(curve->t2, curve->b, curve->p);
+    if (bnCmp (curve->t1, curve->t2) != 0) {
+        return 0;
+    }
+    return 1;
+
+}
+
+static int ecCheckPubKey3617(const EcCurve *curve, const EcPoint *pub)
+{
+    /* Represent point at infinity by (0, 0), make sure it's not that */
+    if (bnCmpQ(pub->x, 0) == 0 && bnCmpQ(pub->y, 0) == 0) {
+        return 0;
+    }
+    /* Check that coordinates are within range */
+    if (bnCmpQ(pub->x, 0) < 0 || bnCmp(pub->x, curve->p) >= 0) {
+        return 0;
+    }
+    if (bnCmpQ(pub->y, 0) < 0 || bnCmp(pub->y, curve->p) >= 0) {
+        return 0;
+    }
+    /* Check that point satisfies EC equation x^2+y^2 = 1+3617x^2y^2, mod P */
+    bnSquareMod_(curve->t1, pub->y, curve->p, curve);
+    bnSquareMod_(curve->t2, pub->x, curve->p, curve);
+    bnCopy(curve->t3, curve->t1);                                /* Load t3 */
+    bnAddMod_(curve->t3, curve->t2, curve->p);                   /* t3 = t1 + t2, (x^2+y^2)*/
+
+    bnMulMod_(curve->t0, curve->a, curve->t1, curve->p, curve);  /* t0 = a * t1,  (3617 * x^2) */
+    bnMulMod_(curve->t0, curve->t0, curve->t2, curve->p, curve); /* t0 = t0 * t1, (3617 * x^2 * y^2) */
+    bnAddMod_(curve->t0, mpiOne, curve->p);                      /* t0 = t0 + 1,  (3617 * x^2 * y^2 + 1) */
+
+    if (bnCmp (curve->t0, curve->t3) != 0) {
+        return 0;
+    }
+    return 1;
+}
+
+/**
+ * According to http://cr.yp.to/ecdh.html#validate no validation is required if used for Diffie-Hellman
+ * thus always return success.
+ */
+static int ecCheckPubKey25519(const EcCurve *curve, const EcPoint *pub)
+{
+    return 1;
+}
+
+static int mod3617(BigNum *r, const BigNum *a, const BigNum *modulo)
+{
+    unsigned char buffer[52] = {0};
+    int cmp;
+    BigNum tmp;
+
+    bnBegin(&tmp);
+    cmp = bnCmp(modulo, a);
+    if (cmp == 0) {             /* a is equal modulo, set resul to zero */
+        bnSetQ(r, 0);
+        return 0;
+    }
+    if (cmp > 0) {              /* modulo is greater than a - copy a to r and return it */
+        bnCopy(r, a);
+        return 0;
+    }
+    bnExtractLittleBytes(a, buffer, 0, 52);
+    buffer[51] &= 0x3f;
+
+    bnCopy(&tmp, a);
+    bnRShift(&tmp, 414);
+    bnCopy(r, &tmp);
+    bnLShift(r, 4);
+    bnAdd(r, &tmp);
+
+    bnInsertLittleBytes(&tmp, buffer, 0, 52);
+
+    bnAdd(r, &tmp);
+    while (bnCmp(r, modulo) >= 0) {
+        bnSub(r, modulo);
+    }
+    bnEnd(&tmp);
+    return 0;
+}
+
+/* 
+ * Curve25519 has no specific modulo function, all inside curve25519-donna
+ */
+static int mod25519(BigNum *r, const BigNum *a, const BigNum *modulo)
+{
+    return -2;
+}
+
+/*
+ * Beware: Here are the dragons.
+ *
+ * The modulo implementations for the NIST curves. For more detailled information see
+ * FIPS 186-3, chapter D.2 and other papers about Generailzed Mersenne numbers.
+ *
+ * I use byte operations to perfom the additions with carry. On a little endian machine
+ * this saves conversion from/to big endian format if I would use integers for example. Also
+ * using byte addition into a short carry accumulator works on every word size and avoids
+ * complex testing and handling of wordsizes and big/little endian stuff.
+ *
+ */
+
+/* new modulo for 192bit curve */
+static int newMod192(BigNum *r, const BigNum *a, const BigNum *modulo)
+{
+    unsigned char buffer[200] = {0};
+    unsigned char *pt;
+    unsigned char *ps1;
+    unsigned char *ps2;
+    unsigned char *ps3;
+    short ac;
+    int cmp;
+
+    /* Binary big number representation in PolarSSL is always big endian
+     *
+     * the least significant 64bit large word starts at byte offset 40,
+     * the least significant 32bit word starts at byte offset 44
+     * the least significant byte starts at byte offset 47
+     *
+     *           S3    S2   S1          T
+     *                            /-----^------\
+     *           A5    A4   A3    A2    A1    A0
+     * 64bit  0     1     2     3     4     5
+     *        |--+--|--+--|--+--|--+--|--+--|--+--|
+     * 32bit  0  1  2  3  4  5  6  7  8  9 10 11
+     *
+     * perform T + S1 + S2 + S3 mod p
+
+     * where T  = (A2 || A1 || A0)
+     *     + S1 = ( 0 || A3 || A3)
+     *     + S2 = (A4 || A4 ||  0)
+     *     + S3 = (A5 || A5 || A5)
+     *
+     * TODO: error check if input variable is > modulo^2 (do normal mpi_mod_mpi),
+     */
+
+    /* TODO: check if a is > modulo^2 */
+    cmp = bnCmp(modulo, a);
+    if (cmp == 0) {             /* a is equal modulo, set resul to zero */
+        bnSetQ(r, 0);
+        return 0;
+    }
+    if (cmp > 0) {              /* modulo is greater than a - copy a to r and return it */
+        bnCopy(r, a);
+        return 0;
+    }
+    bnExtractBigBytes(a, buffer, 0, bnBytes(modulo)*2);
+
+    /* 6 'A' words, each word is 8 byte. Compute offset to least significant byte of word X */
+#define A(X) buffer + (((6-X)*8)-1)
+
+    ac = 0;
+
+    pt = A(0);      /* pt points to least significant byte of A0  */
+
+    /* Add up first 8 byte word, no need to add ps2 */
+    ps1 = A(3);        /* ps1 points to least significant byte of S1 (A3) */
+    ps3 = A(5);        /* ps3 points to least significant byte of S3 (A5)*/
+
+    /* Each block processes one 32 bit word, big endian, using byte operations */
+    ac += *pt + *ps1--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+
+    ac += *pt + *ps1--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+
+    /* Add up second 8 byte word, all three S words are used here */
+    ps1 = A(3); ps2 = A(4); ps3 = A(5);
+
+    ac += *pt + *ps1--; ac += *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1--; ac += *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1--; ac += *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1--; ac += *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+
+    ac += *pt + *ps1--; ac += *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1--; ac += *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1--; ac += *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1--; ac += *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+
+    /* Add up third 8 byte word, no need to add S1 word */
+    ps2 = A(4); ps3 = A(5);
+
+    ac += *pt + *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+
+    ac += *pt + *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; *pt-- = ac; ac >>= 8;
+
+    /* In this function we cannot have a negative carry and at most a carry of 2
+     * thus just subtract the modulo until we are less than modulo
+     */
+    bnSetQ(r, 0);
+
+    *(A(3)) = ac;      /* Store the carry */
+    bnInsertBigBytes(r, A(3), 0, 25);  /* 25: 3 * 8 byte words + 1 carry byte */
+    while (bnCmp(r, modulo) >= 0) {
+        bnSub(r, modulo);
+    }
+    return 0;
+}
+#undef A
+
+/* new modulo for 256bit curve */
+static int newMod256(BigNum *r, const BigNum *a, const BigNum *modulo)
+{
+    unsigned char buffer[200] = {0};
+    unsigned char *pt;
+    unsigned char *ps1;
+    unsigned char *ps2;
+    unsigned char *ps3;
+    unsigned char *ps4;
+
+    unsigned char *pd1;
+    unsigned char *pd2;
+    unsigned char *pd3;
+    unsigned char *pd4;
+    short ac;
+    int cmp;
+
+    /* Binary big number representation in PolarSSL is always big endian
+     *
+     * the least significant byte starts at byte offset 63
+     *
+     *                                                                    T
+     *                                                  /-----------------^------------------\
+     *          A15  A14  A13  A12  A11  A10  A9   A8   A7   A6   A5   A4   A3   A2   A1   A0
+     *        |----+----|----+----|----+----|----+----|----+----|----+----|----+----|----+----|
+     * offset 0    4    8   12   16   20   24   28   32   36   40   44   48   52   56   60    64
+     *
+     * T  = (  A7 ||  A6 ||  A5 ||  A4 ||  A3 ||  A2 ||  A1 ||  A0 )
+     *
+     * S1 = ( A15 || A14 || A13 || A12 || A11 ||  00 ||  00 ||  00 )
+     * S2 = (  00 || A15 || A14 || A13 || A12 ||  00 ||  00 ||  00 )
+     * S3 = ( A15 || A14 ||  00 ||  00 ||  00 || A10 ||  A9 ||  A8 )
+     * S4 = (  A8 || A13 || A15 || A14 || A13 || A11 || A10 ||  A9 )
+     * D1 = ( A10 ||  A8 ||  00 ||  00 ||  00 || A13 || A12 || A11 )
+     * D2 = ( A11 ||  A9 ||  00 ||  00 || A15 || A14 || A13 || A12 )
+     * D3 = ( A12 ||  00 || A10 ||  A9 ||  A8 || A15 || A14 || A13 )
+     * D4 = ( A13 ||  00 || A11 || A10 ||  A9 ||  00 || A15 || A14 )
+     *
+     * perform B = T + 2*S1 + 2*S2 + S3 + S4 - D1 - D2 - D3 - D4 mod p
+     *
+     * TODO: error check if input variable is > modulo^2 (do normal mpi_mod_mpi),
+     */
+
+    cmp = bnCmp(modulo, a);
+    if (cmp == 0) {             /* a is equal modulo, set resul to zero */
+        bnSetQ(r, 0);
+        return 0;
+    }
+    if (cmp > 0) {              /* modulo is greater than a - copya to r and return it */
+        bnCopy(r, a);
+        return 0;
+    }
+    bnExtractBigBytes(a, buffer, 0, bnBytes(modulo)*2);
+
+    /* 16 'A' words, each word is 4 byte. Compute offset to least significant byte of word X */
+#define A(X) buffer + (((16-X)*4)-1)
+
+    ac = 0;
+
+    pt = A(0);          /* pt points to least significant byte of A0  */
+
+    /* Set up to add up data that goes into A0 (right-most column abover); S1, S2 not used */
+    ps3 = A(8);         /* ps3 points to least significant byte of S3 */
+    ps4 = A(9);         /* ps4 points to least significant byte of S4 */
+    pd1 = A(11);        /* pd1 points to least significant byte of D1 */
+    pd2 = A(12);        /* pd2 points to least significant byte of D2 */
+    pd3 = A(13);        /* pd3 points to least significant byte of D3 */
+    pd4 = A(14);        /* pd4 points to least significant byte of D4 */
+
+    /* Each block processes one 32 bit word, big endian, using byte operations */
+    ac += *pt + *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add up data that goes into A1; S1, S2 not used */
+    ps3 = A(9);  ps4 = A(10); pd1 = A(12); pd2 = A(13); pd3 = A(14); pd4 = A(15);
+    ac += *pt + *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add up data that goes into A2; S1, S2, D4 not used */
+    ps3 = A(10); ps4 = A(11); pd1 = A(13); pd2 = A(14); pd3 = A(15);
+    ac += *pt + *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add up data that goes into A3; S3, D1 not used */
+    ps1 = A(11); ps2 = A(12); ps4 = A(13); pd2 = A(15); pd3 = A(8); pd4 = A(9);
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps4--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps4--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps4--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps4--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add up data that goes into A4; S3, D1, D2 not used */
+    ps1 = A(12); ps2 = A(13); ps4 = A(14); pd3 = A(9); pd4 = A(10);
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps4--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps4--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps4--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps4--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add up data that goes into A5; S3, D1, D2 not used */
+    ps1 = A(13); ps2 = A(14); ps4 = A(15); pd3 = A(10); pd4 = A(11);
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps4--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps4--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps4--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps4--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add up data that goes into A6; D3, D4 not used */
+    ps1 = A(14); ps2 = A(15); ps3 = A(14); ps4 = A(13); pd1 = A(8); pd2 = A(9);
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps2;ac += *ps2--;  ac += *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add up data that goes into A7; S2 not used */
+    ps1 = A(15); ps3 = A(15); ps4 = A(8); pd1 = A(10); pd2 = A(11); pd3 = A(12); pd4 = A(13);
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--;  ac += *ps3--; ac += *ps4--; ac -= *pd1--; ac -= *pd2--; ac -= *pd3--; ac -= *pd4--; *pt-- = ac; ac >>= 8;
+
+    bnSetQ(r, 0);
+    if (ac > 0) {
+        *(A(8)) = ac;      /* Store the carry */
+        bnInsertBigBytes(r, A(8), 0, 33);  /* 33: 8 * 4 byte words + 1 carry byte */
+    }
+    /* Negative carry requires that we add the modulo (carry * -1) times to make
+     * the result positive. Then get the result mod(256).
+     */
+    else if (ac < 0) {
+        int msb, maxMsb;
+
+        *(A(8)) = 0;
+        bnInsertBigBytes(r, A(8), 0, 33);  /* 33: 8 * 4 byte words + 1 carry byte */
+        ac *= -1;
+        while (ac--) {
+            bnAdd(r, modulo);
+        }
+        maxMsb =  bnBits(modulo);
+        msb = bnBits(r) - maxMsb;
+        /* clear all bits above bit length of modulo. This length is 256 here, thus
+         * we effectiviely doing a mod(256)
+         */
+        if (msb > 0) {
+            BigNum tmp;
+            bnBegin(&tmp);
+            bnSetQ (&tmp, 1);
+            bnLShift (&tmp, maxMsb);
+            bnMod(r, r, &tmp);
+            bnEnd(&tmp);
+        }
+    }
+    else {
+        *(A(8)) = 0;
+        bnInsertBigBytes(r, A(8), 0, 33);  /* 33: 8 * 4 byte words + 1 carry byte */
+    }
+    while (bnCmp(r, modulo) >= 0) {
+        bnSub(r, modulo);
+    }
+    return 0;
+}
+#undef A
+
+
+/* new modulo for 384bit curve */
+static int newMod384(BigNum *r, const BigNum *a, const BigNum *modulo)
+{
+    unsigned char buffer[200] = {0};
+    unsigned char *pt;
+    unsigned char *ps1;
+    unsigned char *ps2;
+    unsigned char *ps3;
+    unsigned char *ps4;
+    unsigned char *ps5;
+    unsigned char *ps6;
+
+    unsigned char *pd1;
+    unsigned char *pd2;
+    unsigned char *pd3;
+    short ac;
+    int cmp;
+
+    /*
+     *
+     * the least significant byte starts at byte offset 97
+     *
+     *                                                                    T
+     *                                        /---------------------------^----------------------------\
+     *      A23 ......... A15  A14  A13  A12  A11  A10  A9   A8   A7   A6   A5   A4   A3   A2   A1   A0
+     *    |----+ ...... |----+----|----+----|----+----|----+----|----+----|----+----|----+----|----+----|
+     *
+     * T  = (A11 || A10 ||  A9 ||  A8 ||  A7 ||  A6 ||  A5 ||  A4 ||  A3 ||  A2 ||  A1 ||  A0)
+
+     * S1 = ( 00 ||  00 ||  00 ||  00 ||  00 || A23 || A22 || A21 ||  00 ||  00 ||  00 ||  00)
+     * S2 = (A23 || A22 || A21 || A20 || A19 || A18 || A17 || A16 || A15 || A14 || A13 || A12)
+     * S3 = (A20 || A19 || A18 || A17 || A16 || A15 || A14 || A13 || A12 || A23 || A22 || A21)
+     * S4 = (A19 || A18 || A17 || A16 || A15 || A14 || A13 || A12 || A20 ||  00 || A23 ||  00)
+     * S5 = ( 00 ||  00 ||  00 ||  00 || A23 || A22 || A21 || A20 ||  00 ||  00 ||  00 ||  00)
+     * S6 = ( 00 ||  00 ||  00 ||  00 ||  00 ||  00 || A23 || A22 || A21 ||  00 ||  00 || A20)
+     * D1 = (A22 || A21 || A20 || A19 || A18 || A17 || A16 || A15 || A14 || A13 || A12 || A23)
+     * D2 = ( 00 ||  00 ||  00 ||  00 ||  00 ||  00 ||  00 || A23 || A22 || A21 || A20 ||  00)
+     * D3 = ( 00 ||  00 ||  00 ||  00 ||  00 ||  00 ||  00 || A23 || A23 ||  00 ||  00 ||  00)
+     *
+     * perform B =  T + 2S1 + S2 + S3 + S4 + S5 + S6 – D1 – D2 – D3 mod p
+     *
+     * TODO: error check if input variable is > modulo^2 (do normal mpi_mod_mpi),
+     *       optimize if input is already < modulo  (just copy over in this case).
+     */
+
+    cmp = bnCmp(modulo, a);
+    if (cmp == 0) {             /* a is equal modulo, set resul to zero */
+        bnSetQ(r, 0);
+        return 0;
+    }
+    if (cmp > 0) {              /* modulo is greater than a - copy a to r and return it */
+        bnCopy(r, a);
+        return 0;
+    }
+
+    bnExtractBigBytes(a, buffer, 0, bnBytes(modulo)*2);
+
+    /* 24 'A' words, each word is 4 byte. Compute offset to least significant byte of word X */
+#define A(X) buffer + (((24-X)*4)-1)
+
+    ac = 0;
+
+    pt = A(0);      /* pt points to least significant byte of A0  */
+
+    /* Set up to add data that goes into A0; S1, S4, S5, D2, D3 not used */
+    ps2 = A(12); ps3 = A(21); ps6 = A(20); pd1 = A(23);
+
+    /* Each block processes one 32 bit word, big endian, using byte operations */
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps6--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps6--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps6--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps6--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add data that goes into A1; S1, S5, S6, D3 not used */
+    ps2 = A(13); ps3 = A(22); ps4 = A(23); pd1= A(12); pd2 = A(20);
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--;  ac -= *pd2--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--;  ac -= *pd2--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--;  ac -= *pd2--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--;  ac -= *pd2--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add data that goes into A2; S1, S4, S5, S6, D3 not used */
+    ps2 = A(14); ps3 = A(23); pd1 = A(13); pd2 = A(21);
+    ac += *pt + *ps2--; ac += *ps3--; ac -= *pd1--;  ac -= *pd2--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac -= *pd1--;  ac -= *pd2--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac -= *pd1--;  ac -= *pd2--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac -= *pd1--;  ac -= *pd2--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add data that goes into A3; S1, S5, S6 not used */
+    ps2 = A(15); ps3 = A(12); ps4 = A(20); ps6 = A(21); pd1 = A(14); pd2 = A(22); pd3 = A(23);
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps6--; ac -= *pd1--;  ac -= *pd2--; ac -= *pd3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps6--; ac -= *pd1--;  ac -= *pd2--; ac -= *pd3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps6--; ac -= *pd1--;  ac -= *pd2--; ac -= *pd3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps6--; ac -= *pd1--;  ac -= *pd2--; ac -= *pd3--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add data that goes into A4 */
+    ps1 = A(21); ps2 = A(16); ps3 = A(13); ps4 = A(12); ps5 = A(20); ps6 = A(22); pd1 = A(15); pd2 = A(23), pd3 = A(23);
+    ac += *pt + *ps1;ac += *ps1--; ac += *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac += *ps6--; ac -= *pd1--;  ac -= *pd2--; ac -= *pd3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--; ac += *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac += *ps6--; ac -= *pd1--;  ac -= *pd2--; ac -= *pd3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--; ac += *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac += *ps6--; ac -= *pd1--;  ac -= *pd2--; ac -= *pd3--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--; ac += *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac += *ps6--; ac -= *pd1--;  ac -= *pd2--; ac -= *pd3--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add data that goes into A5; D2, D3 not used */
+    ps1 = A(22); ps2 = A(17); ps3 = A(14); ps4 = A(13); ps5 = A(21); ps6 = A(23); pd1 = A(16);
+    ac += *pt + *ps1;ac += *ps1--; ac += *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac += *ps6--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--; ac += *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac += *ps6--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--; ac += *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac += *ps6--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--; ac += *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac += *ps6--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add data that goes into A6; S6, D2, D3 not used */
+    ps1 = A(23); ps2 = A(18); ps3 = A(15); ps4 = A(14); ps5 = A(22); pd1 = A(17);
+    ac += *pt + *ps1;ac += *ps1--; ac += *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--; ac += *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--; ac += *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps1;ac += *ps1--; ac += *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add data that goes into A7; S1, S6, D2, D3 not used */
+    ps2 = A(19); ps3 = A(16); ps4 = A(15); ps5 = A(23); pd1 = A(18);
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac += *ps5--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add data that goes into A8; S1, S5, S6, D2, D3 not used */
+    ps2 = A(20); ps3 = A(17); ps4 = A(16); pd1 = A(19);
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add data that goes into A9; S1, S5, S6, D2, D3 not used */
+    ps2 = A(21); ps3 = A(18); ps4 = A(17); pd1 = A(20);
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add data that goes into A10; S1, S5, S6, D2, D3 not used */
+    ps2 = A(22); ps3 = A(19); ps4 = A(18); pd1 = A(21);
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+
+    /* Set up to add data that goes into A10; S1, S5, S6, D2, D3 not used */
+    ps2 = A(23); ps3 = A(20); ps4 = A(19); pd1 = A(22);
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+    ac += *pt + *ps2--; ac += *ps3--; ac += *ps4--; ac -= *pd1--; *pt-- = ac; ac >>= 8;
+
+    bnSetQ(r, 0);
+    if (ac > 0) {
+        *(A(12)) = ac;      /* Store the carry */
+        bnInsertBigBytes(r, A(12), 0, 49);  /* 49: 12 * 4 byte words + 1 carry byte */
+    }
+    /* Negative carry requires that we add the modulo (carry * -1) times to make
+     * the result positive. Then get the result mod(256).
+     */
+    else if (ac < 0) {
+        int msb, maxMsb;
+
+        *(A(12)) = 0;
+        bnInsertBigBytes(r, A(12), 0, 49);  /* 49: 12 * 4 byte words + 1 carry byte */
+        ac *= -1;
+        while (ac--) {
+            bnAdd(r, modulo);
+        }
+        maxMsb =  bnBits(modulo);
+        msb = bnBits(r) - maxMsb;
+        /* clear all bits above bit length of modulo. This length is 384 here, thus
+         * we effectiviely doing a mod(384)
+         */
+        if (msb > 0) {
+            BigNum tmp;
+            bnBegin(&tmp);
+            bnSetQ (&tmp, 1);
+            bnLShift (&tmp, maxMsb);
+            bnMod(r, r, &tmp);
+            bnEnd(&tmp);
+        }
+    }
+    else {
+        *(A(12)) = 0;
+        bnInsertBigBytes(r, A(12), 0, 49);  /* 49: 12 * 4 byte words + 1 carry byte */
+    }
+    while (bnCmp(r, modulo) >= 0) {
+        bnSub(r, modulo);
+    }
+    return 0;
+}
+#undef A
+
+
+/* new modulo for 521bit curve, much easier because the prime for 521 is a real Mersenne prime */
+static int newMod521(BigNum *r, const BigNum *a, const BigNum *modulo)
+{
+    unsigned char buf1[200] = {0};
+    unsigned char buf2[200] = {0};
+    unsigned char *p1;
+    unsigned char *p2;
+    size_t modSize;
+    short ac = 0;
+    unsigned int i;
+    int cmp;
+
+    /* TODO: check if a is > modulo^2 */
+#if 0
+    if (a->s < 0)               /* is it a negative value? */
+        return bnMod(r, a, modulo);
+#endif
+    cmp = bnCmp(modulo, a);
+    if (cmp == 0) {             /* a is equal modulo, set resul to zero */
+        bnSetQ(r, 0);
+        return 0;
+    }
+    bnCopy(r, a);
+    if (cmp > 0) {              /* modulo is greater than a - return the prepared r */
+        return 0;
+    }
+    modSize = bnBytes(modulo);
+
+    bnExtractBigBytes(a, buf1, 0, modSize*2); /* a must be less modulo^2 */
+    buf1[modSize] &= 1;                   /* clear all bits except least significat */
+
+    bnRShift(r, 521);
+    bnExtractBigBytes(r, buf2, 0, modSize*2);
+    buf2[modSize] &= 1;
+
+    p1 = &buf2[131];            /* p1 is pointer to A0 */
+    p2 = &buf1[131];            /* p2 is pointer to A1 */
+
+    for (i = 0; i < modSize; i++) {
+        ac += *p1 + *p2--; *p1-- = ac; ac >>= 8;
+    }
+    bnSetQ(r, 0);
+    bnInsertBigBytes(r, p1+1, 0, modSize);
+
+    while (bnCmp(r, modulo) >= 0) {
+        bnSub(r, modulo);
+    }
+    return 0;
+}
+
diff --git a/jni/libzrtp/sources/bnlib/ec/ec.h b/jni/libzrtp/sources/bnlib/ec/ec.h
index b1bb47c..172ffd8 100644
--- a/jni/libzrtp/sources/bnlib/ec/ec.h
+++ b/jni/libzrtp/sources/bnlib/ec/ec.h
@@ -29,18 +29,34 @@
     NIST224P = 2,
     NIST256P = 3,
     NIST384P = 4,
-    NIST521P = 5
-} NistCurves;
+    NIST521P = 5,
+    Curve25519 = 10,
+    Curve3617  = 11
+} Curves;
 
 /**
- * @brief This structure contains the value of NIST EC curves over Prime Fields.
- *
- * The <b>a</b> curve parameter is the constant -3 and is computed during initialization
- * of the curve structure.
- *
- * The field names correspond to the variable names defined in NIST FIPS 186-3, E.1.2
+ * \brief This structure contains the x, y affine coordinates and the z value if we
+ *        use projective coordinates during EC point arithmetic.
  */
-typedef struct {
+typedef struct _EcPoint {
+    BigNum *x, *y, *z;
+    BigNum tx, ty, tz;
+} EcPoint;
+
+/**
+ * @brief This structure contains the value of EC curves over Prime Fields.
+ *
+ * The for NIST curves the field names correspond to the variable names defined in 
+ * NIST FIPS 186-3, E.1.2. The <b>a</b> curve parameter is the constant -3 and is 
+ * computed during initialization of the curve structure.
+ *
+ * For other curves, for example curve3917 we have less parameters to fill in, mostly
+ * the prime number, the base point, etc. Refer to the curve's initialization function
+ * about the use of the fileds.
+ */
+struct EcCurve;
+struct EcCurve {
+    Curves id;
     BigNum _p;
     BigNum _n;
     BigNum _SEED;
@@ -62,16 +78,18 @@
        avoid to much memory allocation/deallocatio0n overhead */
   BigNum _S1, _U1, _H, _R, _t0, _t1, _t2, _t3;
   BigNum *S1, *U1, *H, *R, *t0, *t1, *t2, *t3;
-} NistECpCurve;
+  int (*affineOp)(const struct EcCurve *curve, EcPoint *R, const EcPoint *P);
+  int (*doubleOp)(const struct EcCurve *curve, EcPoint *R, const EcPoint *P);
+  int (*addOp)(const struct EcCurve *curve, EcPoint *R, const EcPoint *P, const EcPoint *Q);
+  int (*modOp)(BigNum *, const BigNum *, const BigNum *);
+  int (*checkPubOp)(const struct EcCurve *curve, const EcPoint *pub);
+  int (*randomOp)(const struct EcCurve *curve, BigNum *d);
+  int (*mulScalar)(const struct EcCurve *curve, EcPoint *R, const EcPoint *P, const BigNum *scalar);
 
-/**
- * \brief This structure contains the x, y affine coordinates and the z value if we
- *        use projective coordinates during EC point arithmetic.
- */
-typedef struct _EcPoint {
-    BigNum *x, *y, *z;
-    BigNum tx, ty, tz;
-} EcPoint;
+};
+
+typedef struct EcCurve EcCurve;
+typedef EcCurve NistECpCurve;
 
 /**
  * \brief          Marco to initialize a EC point structure.
@@ -94,7 +112,16 @@
  *
  * \param P        Address of the EC point structure.
  */
-#define SET_EC_BASE_POINT(C, P) {EcPoint *e = P;  const NistECpCurve *c = C; bnCopy(e->x, c->Gx); bnCopy(e->y, c->Gy); bnSetQ(e->z, 1);}
+#define SET_EC_BASE_POINT(C, P) {EcPoint *e = P;  const EcCurve *c = C; bnCopy(e->x, c->Gx); bnCopy(e->y, c->Gy); bnSetQ(e->z, 1);}
+
+/*
+ * EC point helper functions
+ */
+extern void ecInitPoint(EcPoint *P);
+
+extern void ecFreePoint(EcPoint *P);
+
+extern void ecSetBasePoint(EcCurve *C, EcPoint *P);
 
 /**
  * \brief          Get NIST EC curve parameters.
@@ -104,23 +131,23 @@
  *
  * \param curveId  Which curve to initialize
  *
- * \param curve    Pointer to a NistECpCurve structure
+ * \param curve    Pointer to a EcCurve structure
  *
- * \return         0 if successful, or a POLARSSL_ERR_EC_XXX/ POLARSSL_ERR_MPI_XXX error code.
+ * \return         0 if successful
  *
  * \note           Call ecFreeCurveNistECp to return allocated memory.
  */
-int ecGetCurveNistECp(NistCurves curveId, NistECpCurve *curve);
+int ecGetCurveNistECp(Curves curveId, NistECpCurve *curve);
 
 
 /**
- * \brief          Free NIST EC curve parameters.
+ * \brief          Free EC curve parameters.
  *
- * \param curve    Pointer to a NistECpCurve structure
+ * \param curve    Pointer to a EcCurve structure
  *
  * \note           Curve parameters must be initialized calling ecGetCurveNistECp.
  */
-void ecFreeCurveNistECp(NistECpCurve *curve);
+void ecFreeCurveNistECp(EcCurve *curve);
 
 /**
  * \brief          Double an EC point.
@@ -129,13 +156,13 @@
  *                 further reference see RFC 6090 or the standard work <i>Guide to Elliptic
  *                 Curve Cryptography</i>.
  *
- * \param          curve  Address of Nist EC curve structure
+ * \param          curve  Address of EC curve structure
  * \param          R      Address of resulting EC point structure
  * \param          P      Address of the EC point structure
  *
- * \return         0 if successful, or a POLARSSL_ERR_EC_XXX / POLARSSL_ERR_MPI_XXX error code.
+ * \return         0 if successful
  */
-int ecDoublePoint(const NistECpCurve *curve, EcPoint *R, const EcPoint *P);
+int ecDoublePoint(const EcCurve *curve, EcPoint *R, const EcPoint *P);
 
 /**
  * \brief          Add two EC points.
@@ -144,43 +171,43 @@
  *                 further reference see RFC 6090 or the standard work <i>Guide to Elliptic
  *                 Curve Cryptography</i>.
  *
- * \param          curve  Address of Nist EC curve structure
+ * \param          curve  Address of EC curve structure
  * \param          R      Address of resulting EC point structure
  * \param          P      Address of the first EC point structure
  * \param          Q      Address of the second EC point structure
  *
- * \return         0 if successful, or a POLARSSL_ERR_EC_XXX / POLARSSL_ERR_MPI_XXX error code.
+ * \return         0 if successful
  */
-int ecAddPoint(const NistECpCurve *curve, EcPoint *R, const EcPoint *P, const EcPoint *Q);
+int ecAddPoint(const EcCurve *curve, EcPoint *R, const EcPoint *P, const EcPoint *Q);
 
 /**
  * \brief          Mulitply an EC point with a scalar value.
  *
- * \param          curve  Address of Nist EC curve structure
+ * \param          curve  Address of EC curve structure
  * \param          R      Address of resulting EC point structure
  * \param          P      Address of the EC point structure
  * \param          scalar Address of the scalar multi-precision integer value
  *
- * \return         0 if successful, or a POLARSSL_ERR_EC_XXX / POLARSSL_ERR_MPI_XXX error code.
+ * \return         0 if successful
  */
-int ecMulPointScalar(const NistECpCurve *curve, EcPoint *R, const EcPoint *P, const BigNum *scalar);
+int ecMulPointScalar(const EcCurve *curve, EcPoint *R, const EcPoint *P, const BigNum *scalar);
 
 /**
  * \brief          Convert an EC point from Jacobian projective coordinates to normal affine x/y coordinates.
  *
- * \param          curve  Address of Nist EC curve structure
+ * \param          curve  Address of EC curve structure
  * \param          R      Address of EC point structure that receives the x/y coordinates
  * \param          P      Address of the EC point structure that contains the jacobian x/y/z coordinates.
  *
- * \return         0 if successful, or a POLARSSL_ERR_EC_XXX / POLARSSL_ERR_MPI_XXX error code.
+ * \return         0 if successful
  */
-int ecGetAffine(const NistECpCurve *curve, EcPoint *R, const EcPoint *P);
+int ecGetAffine(const EcCurve *curve, EcPoint *R, const EcPoint *P);
 
 /**
  * @brief Generate a random number.
  *
  * The method generates a random number and checks if it matches the curve restricitions.
- * Use this number to generate a ECDH public key.
+ * Use this number as ECDH private key.
  *
  * @param curve the NIST curve to use.
  *
@@ -188,6 +215,35 @@
  */
 int ecGenerateRandomNumber(const NistECpCurve *curve, BigNum *d);
 
+/**
+ * @brief Check a public key.
+ *
+ * The method checks if a public key is valid. For NIST curves it uses the
+ * ECC Partial Validation, NIST SP800-56A section 5.6.2.6
+ * 
+ * For other curves it computes the equation and compares the left hand and 
+ * the right handresults. If they are equal the point is on the curve.
+ *
+ * @param curve the curve to use.
+ *
+ * @param pub the public key to check.
+ *
+ * @returns true (!0) if the check was ok, false (0) otherwise.
+ *
+ * @note The function uses some scratch pad variable of the NistECpCurve structure.
+ */
+int ecCheckPubKey(const EcCurve *curve, const EcPoint *pub);
+
+int ecGetCurvesCurve(Curves curveId, EcCurve *curve);
+
+void ecFreeCurvesCurve(EcCurve *curve);
+
+/**
+ * This is a special function for DJB's curve 25519. Actually it's the scalar multiplication
+ * mypublic = basepoint * secret
+ */
+int curve25519_donna(unsigned char *mypublic, const unsigned char *secret, const unsigned char *basepoint);
+
 /*
  * Some additional functions that are not available in bnlib
  */
@@ -199,11 +255,11 @@
 
 int bnSubQMod_ (struct BigNum *rslt, unsigned n1, struct BigNum *mod);
 
-int bnMulMod_ (struct BigNum *rslt, struct BigNum *n1, struct BigNum *n2, struct BigNum *mod);
+int bnMulMod_ (struct BigNum *rslt, struct BigNum *n1, struct BigNum *n2, struct BigNum *mod, const EcCurve *curve);
 
-int bnMulQMod_ (struct BigNum *rslt, struct BigNum *n1, unsigned n2, struct BigNum *mod);
+int bnMulQMod_ (struct BigNum *rslt, struct BigNum *n1, unsigned n2, struct BigNum *mod, const EcCurve *curve);
 
-int bnSquareMod_ (struct BigNum *rslt, struct BigNum *n1, struct BigNum *mod);
+int bnSquareMod_ (struct BigNum *rslt, struct BigNum *n1, struct BigNum *mod, const EcCurve *curve);
 
 #ifdef __cplusplus
 }
diff --git a/jni/libzrtp/sources/bnlib/ec/ecdh.c b/jni/libzrtp/sources/bnlib/ec/ecdh.c
index 25ea8cd..8d1bc23 100644
--- a/jni/libzrtp/sources/bnlib/ec/ecdh.c
+++ b/jni/libzrtp/sources/bnlib/ec/ecdh.c
@@ -9,7 +9,7 @@
 #include <ec/ec.h>
 #include <ec/ecdh.h>
 
-int ecdhGeneratePublic(const NistECpCurve *curve, EcPoint *Q, const BigNum *d)
+int ecdhGeneratePublic(const EcCurve *curve, EcPoint *Q, const BigNum *d)
 {
     EcPoint G;
 
@@ -21,10 +21,10 @@
 
     FREE_EC_POINT(&G);
 
-    return 0;
+    return ecCheckPubKey(curve, Q);
 }
 
-int ecdhComputeAgreement(const NistECpCurve *curve, BigNum *agreement, const EcPoint *Q, const BigNum *d)
+int ecdhComputeAgreement(const EcCurve *curve, BigNum *agreement, const EcPoint *Q, const BigNum *d)
 {
     EcPoint t0;
 
diff --git a/jni/libzrtp/sources/bnlib/ec/ecdh.h b/jni/libzrtp/sources/bnlib/ec/ecdh.h
index 0cf083a..7ec32ad 100644
--- a/jni/libzrtp/sources/bnlib/ec/ecdh.h
+++ b/jni/libzrtp/sources/bnlib/ec/ecdh.h
@@ -22,15 +22,17 @@
 /**
  * @brief Takes a secret large random number and computes the public EC point.
  *
- * @param curve is the NIST curve to use.
+ * @param curve is the curve to use.
  *
  * @param Q the functions writes the computed public point in this parameter.
  *
  * @param d is the secret random number.
  *
+ * @return @c true (!0) if public key was computed, @c false otherwise.
+ *
  * @sa ecGenerateRandomNumber
  */
-int ecdhGeneratePublic(const NistECpCurve *curve, EcPoint *Q, const BigNum *d);
+int ecdhGeneratePublic(const EcCurve *curve, EcPoint *Q, const BigNum *d);
 
 /**
  * @brief Computes the key agreement value.
@@ -38,7 +40,7 @@
  * Takes the public EC point of the other party and applies the EC DH algorithm
  * to compute the agreed value.
  *
- * @param curve is the NIST curve to use, must be the same curve as used in
+ * @param curve is the curve to use, must be the same curve as used in
  *              @c ecdhGeneratePublic.
  *
  * @param agreemtn the functions writes the computed agreed value in this parameter.
@@ -47,7 +49,7 @@
  *
  * @param d is the secret random number.
  */
-int ecdhComputeAgreement(const NistECpCurve *curve, BigNum *agreement, const EcPoint *Q, const BigNum *d);
+int ecdhComputeAgreement(const EcCurve *curve, BigNum *agreement, const EcPoint *Q, const BigNum *d);
 
 #ifdef __cplusplus
 }