Move all BigInt stuff into bigint/. Currently all asm modules are disabled;

configure.pl doesn't understand how to handle this yet (replace logic only
understands stuff in src, not how one module can replace another modules
src, or anything about prioritizing).

Move some hex and base64 stuff out of charset.cpp and into their
codec directories.

42 files changed, 5212 insertions, 0 deletions

diff --git a/src/bigint/asm_amd64/modinfo.txt b/src/bigint/asm_amd64/xxxinfo.txt
index 2a8f9fe5b..2a8f9fe5b 100644
--- a/src/bigint/asm_amd64/modinfo.txt
+++ b/src/bigint/asm_amd64/xxxinfo.txt
diff --git a/src/bigint/asm_ia32/modinfo.txt b/src/bigint/asm_ia32/xxxinfo.txt
index 12c8cd96d..12c8cd96d 100644
--- a/src/bigint/asm_ia32/modinfo.txt
+++ b/src/bigint/asm_ia32/xxxinfo.txt
diff --git a/src/bigint/big_code.cpp b/src/bigint/big_code.cpp
new file mode 100644
index 000000000..824cbb63e
--- /dev/null
+++ b/src/bigint/big_code.cpp
@@ -0,0 +1,152 @@
+/*************************************************
+* BigInt Encoding/Decoding Source File           *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/bigint.h>
+#include <botan/numthry.h>
+#include <botan/charset.h>
+#include <botan/hex.h>
+
+namespace Botan {
+
+/*************************************************
+* Encode a BigInt                                *
+*************************************************/
+void BigInt::encode(byte output[], const BigInt& n, Base base)
+   {
+   if(base == Binary)
+      n.binary_encode(output);
+   else if(base == Hexadecimal)
+      {
+      SecureVector<byte> binary(n.encoded_size(Binary));
+      n.binary_encode(binary);
+      for(u32bit j = 0; j != binary.size(); ++j)
+         Hex_Encoder::encode(binary[j], output + 2*j);
+      }
+   else if(base == Octal)
+      {
+      BigInt copy = n;
+      const u32bit output_size = n.encoded_size(Octal);
+      for(u32bit j = 0; j != output_size; ++j)
+         {
+         output[output_size - 1 - j] = Charset::digit2char(copy % 8);
+         copy /= 8;
+         }
+      }
+   else if(base == Decimal)
+      {
+      BigInt copy = n;
+      BigInt remainder;
+      copy.set_sign(Positive);
+      const u32bit output_size = n.encoded_size(Decimal);
+      for(u32bit j = 0; j != output_size; ++j)
+         {
+         divide(copy, 10, copy, remainder);
+         output[output_size - 1 - j] =
+            Charset::digit2char(remainder.word_at(0));
+         if(copy.is_zero())
+            break;
+         }
+      }
+   else
+      throw Invalid_Argument("Unknown BigInt encoding method");
+   }
+
+/*************************************************
+* Encode a BigInt                                *
+*************************************************/
+SecureVector<byte> BigInt::encode(const BigInt& n, Base base)
+   {
+   SecureVector<byte> output(n.encoded_size(base));
+   encode(output, n, base);
+   if(base != Binary)
+      for(u32bit j = 0; j != output.size(); ++j)
+         if(output[j] == 0)
+            output[j] = '0';
+   return output;
+   }
+
+/*************************************************
+* Encode a BigInt, with leading 0s if needed     *
+*************************************************/
+SecureVector<byte> BigInt::encode_1363(const BigInt& n, u32bit bytes)
+   {
+   const u32bit n_bytes = n.bytes();
+   if(n_bytes > bytes)
+      throw Encoding_Error("encode_1363: n is too large to encode properly");
+
+   const u32bit leading_0s = bytes - n_bytes;
+
+   SecureVector<byte> output(bytes);
+   encode(output + leading_0s, n, Binary);
+   return output;
+   }
+
+/*************************************************
+* Decode a BigInt                                *
+*************************************************/
+BigInt BigInt::decode(const MemoryRegion<byte>& buf, Base base)
+   {
+   return BigInt::decode(buf, buf.size(), base);
+   }
+
+/*************************************************
+* Decode a BigInt                                *
+*************************************************/
+BigInt BigInt::decode(const byte buf[], u32bit length, Base base)
+   {
+   BigInt r;
+   if(base == Binary)
+      r.binary_decode(buf, length);
+   else if(base == Hexadecimal)
+      {
+      SecureVector<byte> hex;
+      for(u32bit j = 0; j != length; ++j)
+         if(Hex_Decoder::is_valid(buf[j]))
+            hex.append(buf[j]);
+
+      u32bit offset = (hex.size() % 2);
+      SecureVector<byte> binary(hex.size() / 2 + offset);
+
+      if(offset)
+         {
+         byte temp[2] = { '0', hex[0] };
+         binary[0] = Hex_Decoder::decode(temp);
+         }
+
+      for(u32bit j = offset; j != binary.size(); ++j)
+         binary[j] = Hex_Decoder::decode(hex+2*j-offset);
+      r.binary_decode(binary, binary.size());
+      }
+   else if(base == Decimal || base == Octal)
+      {
+      const u32bit RADIX = ((base == Decimal) ? 10 : 8);
+      for(u32bit j = 0; j != length; ++j)
+         {
+         if(Charset::is_space(buf[j]))
+            continue;
+
+         if(!Charset::is_digit(buf[j]))
+            throw Invalid_Argument("BigInt::decode: "
+                                   "Invalid character in decimal input");
+
+         byte x = Charset::char2digit(buf[j]);
+         if(x >= RADIX)
+            {
+            if(RADIX == 10)
+               throw Invalid_Argument("BigInt: Invalid decimal string");
+            else
+               throw Invalid_Argument("BigInt: Invalid octal string");
+            }
+
+         r *= RADIX;
+         r += x;
+         }
+      }
+   else
+      throw Invalid_Argument("Unknown BigInt decoding method");
+   return r;
+   }
+
+}
diff --git a/src/bigint/big_io.cpp b/src/bigint/big_io.cpp
new file mode 100644
index 000000000..3c201e8b2
--- /dev/null
+++ b/src/bigint/big_io.cpp
@@ -0,0 +1,53 @@
+/*************************************************
+* BigInt Input/Output Source File                *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/bigint.h>
+#include <iostream>
+
+namespace Botan {
+
+/*************************************************
+* Write the BigInt into a stream                 *
+*************************************************/
+std::ostream& operator<<(std::ostream& stream, const BigInt& n)
+   {
+   BigInt::Base base = BigInt::Decimal;
+   if(stream.flags() & std::ios::hex)
+      base = BigInt::Hexadecimal;
+   else if(stream.flags() & std::ios::oct)
+      base = BigInt::Octal;
+
+   if(n == 0)
+      stream.write("0", 1);
+   else
+      {
+      if(n < 0)
+         stream.write("-", 1);
+      SecureVector<byte> buffer = BigInt::encode(n, base);
+      u32bit skip = 0;
+      while(buffer[skip] == '0' && skip < buffer.size())
+         ++skip;
+      stream.write(reinterpret_cast<const char*>(buffer.begin()) + skip,
+                   buffer.size() - skip);
+      }
+   if(!stream.good())
+      throw Stream_IO_Error("BigInt output operator has failed");
+   return stream;
+   }
+
+/*************************************************
+* Read the BigInt from a stream                  *
+*************************************************/
+std::istream& operator>>(std::istream& stream, BigInt& n)
+   {
+   std::string str;
+   std::getline(stream, str);
+   if(stream.bad() || (stream.fail() && !stream.eof()))
+      throw Stream_IO_Error("BigInt input operator has failed");
+   n = BigInt(str);
+   return stream;
+   }
+
+}
diff --git a/src/bigint/big_ops2.cpp b/src/bigint/big_ops2.cpp
new file mode 100644
index 000000000..ef083f394
--- /dev/null
+++ b/src/bigint/big_ops2.cpp
@@ -0,0 +1,222 @@
+/*************************************************
+* BigInt Assignment Operators Source File        *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/bigint.h>
+#include <botan/mp_core.h>
+#include <botan/bit_ops.h>
+#include <algorithm>
+
+namespace Botan {
+
+/*************************************************
+* Addition Operator                              *
+*************************************************/
+BigInt& BigInt::operator+=(const BigInt& y)
+   {
+   const u32bit x_sw = sig_words(), y_sw = y.sig_words();
+
+   const u32bit reg_size = std::max(x_sw, y_sw) + 1;
+   grow_to(reg_size);
+
+   if(sign() == y.sign())
+      bigint_add2(get_reg(), reg_size - 1, y.data(), y_sw);
+   else
+      {
+      s32bit relative_size = bigint_cmp(data(), x_sw, y.data(), y_sw);
+
+      if(relative_size < 0)
+         {
+         SecureVector<word> z(reg_size - 1);
+         bigint_sub3(z, y.data(), reg_size - 1, data(), x_sw);
+         copy_mem(get_reg().begin(), z.begin(), z.size());
+         set_sign(y.sign());
+         }
+      else if(relative_size == 0)
+         {
+         get_reg().clear();
+         set_sign(Positive);
+         }
+      else if(relative_size > 0)
+         bigint_sub2(get_reg(), x_sw, y.data(), y_sw);
+      }
+
+   return (*this);
+   }
+
+/*************************************************
+* Subtraction Operator                           *
+*************************************************/
+BigInt& BigInt::operator-=(const BigInt& y)
+   {
+   const u32bit x_sw = sig_words(), y_sw = y.sig_words();
+
+   s32bit relative_size = bigint_cmp(data(), x_sw, y.data(), y_sw);
+
+   const u32bit reg_size = std::max(x_sw, y_sw) + 1;
+   grow_to(reg_size);
+
+   if(relative_size < 0)
+      {
+      if(sign() == y.sign())
+         {
+         SecureVector<word> z(reg_size - 1);
+         bigint_sub3(z, y.data(), reg_size - 1, data(), x_sw);
+         copy_mem(get_reg().begin(), z.begin(), z.size());
+         }
+      else
+         bigint_add2(get_reg(), reg_size - 1, y.data(), y_sw);
+
+      set_sign(y.reverse_sign());
+      }
+   else if(relative_size == 0)
+      {
+      if(sign() == y.sign())
+         {
+         get_reg().clear();
+         set_sign(Positive);
+         }
+      else
+         bigint_shl1(get_reg(), x_sw, 0, 1);
+      }
+   else if(relative_size > 0)
+      {
+      if(sign() == y.sign())
+         bigint_sub2(get_reg(), x_sw, y.data(), y_sw);
+      else
+         bigint_add2(get_reg(), reg_size - 1, y.data(), y_sw);
+      }
+
+   return (*this);
+   }
+
+/*************************************************
+* Multiplication Operator                        *
+*************************************************/
+BigInt& BigInt::operator*=(const BigInt& y)
+   {
+   const u32bit x_sw = sig_words(), y_sw = y.sig_words();
+   set_sign((sign() == y.sign()) ? Positive : Negative);
+
+   if(x_sw == 0 || y_sw == 0)
+      {
+      get_reg().clear();
+      set_sign(Positive);
+      }
+   else if(x_sw == 1 && y_sw)
+      {
+      grow_to(y_sw + 2);
+      bigint_linmul3(get_reg(), y.data(), y_sw, word_at(0));
+      }
+   else if(y_sw == 1 && x_sw)
+      {
+      grow_to(x_sw + 2);
+      bigint_linmul2(get_reg(), x_sw, y.word_at(0));
+      }
+   else
+      {
+      grow_to(size() + y.size());
+
+      SecureVector<word> z(data(), x_sw);
+      SecureVector<word> workspace(size());
+
+      bigint_mul(get_reg(), size(), workspace,
+                 z, z.size(), x_sw,
+                 y.data(), y.size(), y_sw);
+      }
+
+   return (*this);
+   }
+
+/*************************************************
+* Division Operator                              *
+*************************************************/
+BigInt& BigInt::operator/=(const BigInt& y)
+   {
+   if(y.sig_words() == 1 && power_of_2(y.word_at(0)))
+      (*this) >>= (y.bits() - 1);
+   else
+      (*this) = (*this) / y;
+   return (*this);
+   }
+
+/*************************************************
+* Modulo Operator                                *
+*************************************************/
+BigInt& BigInt::operator%=(const BigInt& mod)
+   {
+   return (*this = (*this) % mod);
+   }
+
+/*************************************************
+* Modulo Operator                                *
+*************************************************/
+word BigInt::operator%=(word mod)
+   {
+   if(mod == 0)
+      throw BigInt::DivideByZero();
+   if(power_of_2(mod))
+       {
+       word result = (word_at(0) & (mod - 1));
+       clear();
+       grow_to(2);
+       get_reg()[0] = result;
+       return result;
+       }
+
+   word remainder = 0;
+
+   for(u32bit j = sig_words(); j > 0; --j)
+      remainder = bigint_modop(remainder, word_at(j-1), mod);
+   clear();
+   grow_to(2);
+
+   if(remainder && sign() == BigInt::Negative)
+      get_reg()[0] = mod - remainder;
+   else
+      get_reg()[0] = remainder;
+
+   set_sign(BigInt::Positive);
+
+   return word_at(0);
+   }
+
+/*************************************************
+* Left Shift Operator                            *
+*************************************************/
+BigInt& BigInt::operator<<=(u32bit shift)
+   {
+   if(shift)
+      {
+      const u32bit shift_words = shift / MP_WORD_BITS,
+                   shift_bits  = shift % MP_WORD_BITS,
+                   words = sig_words();
+
+      grow_to(words + shift_words + (shift_bits ? 1 : 0));
+      bigint_shl1(get_reg(), words, shift_words, shift_bits);
+      }
+
+   return (*this);
+   }
+
+/*************************************************
+* Right Shift Operator                           *
+*************************************************/
+BigInt& BigInt::operator>>=(u32bit shift)
+   {
+   if(shift)
+      {
+      const u32bit shift_words = shift / MP_WORD_BITS,
+                   shift_bits  = shift % MP_WORD_BITS;
+
+      bigint_shr1(get_reg(), sig_words(), shift_words, shift_bits);
+
+      if(is_zero())
+         set_sign(Positive);
+      }
+
+   return (*this);
+   }
+
+}
diff --git a/src/bigint/big_ops3.cpp b/src/bigint/big_ops3.cpp
new file mode 100644
index 000000000..7f412f6db
--- /dev/null
+++ b/src/bigint/big_ops3.cpp
@@ -0,0 +1,188 @@
+/*************************************************
+* BigInt Binary Operators Source File            *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/bigint.h>
+#include <botan/numthry.h>
+#include <botan/mp_core.h>
+#include <botan/bit_ops.h>
+#include <algorithm>
+
+namespace Botan {
+
+/*************************************************
+* Addition Operator                              *
+*************************************************/
+BigInt operator+(const BigInt& x, const BigInt& y)
+   {
+   const u32bit x_sw = x.sig_words(), y_sw = y.sig_words();
+
+   BigInt z(x.sign(), std::max(x_sw, y_sw) + 1);
+
+   if((x.sign() == y.sign()))
+      bigint_add3(z.get_reg(), x.data(), x_sw, y.data(), y_sw);
+   else
+      {
+      s32bit relative_size = bigint_cmp(x.data(), x_sw, y.data(), y_sw);
+
+      if(relative_size < 0)
+         {
+         bigint_sub3(z.get_reg(), y.data(), y_sw, x.data(), x_sw);
+         z.set_sign(y.sign());
+         }
+      else if(relative_size == 0)
+         z.set_sign(BigInt::Positive);
+      else if(relative_size > 0)
+         bigint_sub3(z.get_reg(), x.data(), x_sw, y.data(), y_sw);
+      }
+
+   return z;
+   }
+
+/*************************************************
+* Subtraction Operator                           *
+*************************************************/
+BigInt operator-(const BigInt& x, const BigInt& y)
+   {
+   const u32bit x_sw = x.sig_words(), y_sw = y.sig_words();
+
+   s32bit relative_size = bigint_cmp(x.data(), x_sw, y.data(), y_sw);
+
+   BigInt z(BigInt::Positive, std::max(x_sw, y_sw) + 1);
+
+   if(relative_size < 0)
+      {
+      if(x.sign() == y.sign())
+         bigint_sub3(z.get_reg(), y.data(), y_sw, x.data(), x_sw);
+      else
+         bigint_add3(z.get_reg(), x.data(), x_sw, y.data(), y_sw);
+      z.set_sign(y.reverse_sign());
+      }
+   else if(relative_size == 0)
+      {
+      if(x.sign() != y.sign())
+         bigint_shl2(z.get_reg(), x.data(), x_sw, 0, 1);
+      }
+   else if(relative_size > 0)
+      {
+      if(x.sign() == y.sign())
+         bigint_sub3(z.get_reg(), x.data(), x_sw, y.data(), y_sw);
+      else
+         bigint_add3(z.get_reg(), x.data(), x_sw, y.data(), y_sw);
+      z.set_sign(x.sign());
+      }
+   return z;
+   }
+
+/*************************************************
+* Multiplication Operator                        *
+*************************************************/
+BigInt operator*(const BigInt& x, const BigInt& y)
+   {
+   const u32bit x_sw = x.sig_words(), y_sw = y.sig_words();
+
+   BigInt z(BigInt::Positive, x.size() + y.size());
+
+   if(x_sw == 1 && y_sw)
+      bigint_linmul3(z.get_reg(), y.data(), y_sw, x.word_at(0));
+   else if(y_sw == 1 && x_sw)
+      bigint_linmul3(z.get_reg(), x.data(), x_sw, y.word_at(0));
+   else if(x_sw && y_sw)
+      {
+      SecureVector<word> workspace(z.size());
+      bigint_mul(z.get_reg(), z.size(), workspace,
+                 x.data(), x.size(), x_sw,
+                 y.data(), y.size(), y_sw);
+      }
+
+   if(x_sw && y_sw && x.sign() != y.sign())
+      z.flip_sign();
+   return z;
+   }
+
+/*************************************************
+* Division Operator                              *
+*************************************************/
+BigInt operator/(const BigInt& x, const BigInt& y)
+   {
+   BigInt q, r;
+   divide(x, y, q, r);
+   return q;
+   }
+
+/*************************************************
+* Modulo Operator                                *
+*************************************************/
+BigInt operator%(const BigInt& n, const BigInt& mod)
+   {
+   if(mod.is_zero())
+      throw BigInt::DivideByZero();
+   if(mod.is_negative())
+      throw Invalid_Argument("BigInt::operator%: modulus must be > 0");
+   if(n.is_positive() && mod.is_positive() && n < mod)
+      return n;
+
+   BigInt q, r;
+   divide(n, mod, q, r);
+   return r;
+   }
+
+/*************************************************
+* Modulo Operator                                *
+*************************************************/
+word operator%(const BigInt& n, word mod)
+   {
+   if(mod == 0)
+      throw BigInt::DivideByZero();
+   if(power_of_2(mod))
+      return (n.word_at(0) & (mod - 1));
+
+   word remainder = 0;
+
+   for(u32bit j = n.sig_words(); j > 0; --j)
+      remainder = bigint_modop(remainder, n.word_at(j-1), mod);
+
+   if(remainder && n.sign() == BigInt::Negative)
+      return mod - remainder;
+   return remainder;
+   }
+
+/*************************************************
+* Left Shift Operator                            *
+*************************************************/
+BigInt operator<<(const BigInt& x, u32bit shift)
+   {
+   if(shift == 0)
+      return x;
+
+   const u32bit shift_words = shift / MP_WORD_BITS,
+                shift_bits  = shift % MP_WORD_BITS;
+
+   const u32bit x_sw = x.sig_words();
+
+   BigInt y(x.sign(), x_sw + shift_words + (shift_bits ? 1 : 0));
+   bigint_shl2(y.get_reg(), x.data(), x_sw, shift_words, shift_bits);
+   return y;
+   }
+
+/*************************************************
+* Right Shift Operator                           *
+*************************************************/
+BigInt operator>>(const BigInt& x, u32bit shift)
+   {
+   if(shift == 0)
+      return x;
+   if(x.bits() <= shift)
+      return 0;
+
+   const u32bit shift_words = shift / MP_WORD_BITS,
+                shift_bits  = shift % MP_WORD_BITS,
+                x_sw = x.sig_words();
+
+   BigInt y(x.sign(), x_sw - shift_words);
+   bigint_shr2(y.get_reg(), x.data(), x_sw, shift_words, shift_bits);
+   return y;
+   }
+
+}
diff --git a/src/bigint/big_rand.cpp b/src/bigint/big_rand.cpp
new file mode 100644
index 000000000..b8cad3a4c
--- /dev/null
+++ b/src/bigint/big_rand.cpp
@@ -0,0 +1,76 @@
+/*************************************************
+* BigInt Random Generation Source File           *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/bigint.h>
+#include <botan/parsing.h>
+#include <botan/numthry.h>
+
+namespace Botan {
+
+/*************************************************
+* Construct a BigInt of a specific form          *
+*************************************************/
+BigInt::BigInt(NumberType type, u32bit bits)
+   {
+   set_sign(Positive);
+
+   if(type == Power2)
+      set_bit(bits);
+   else
+      throw Invalid_Argument("BigInt(NumberType): Unknown type");
+   }
+
+/*************************************************
+* Randomize this number                          *
+*************************************************/
+void BigInt::randomize(RandomNumberGenerator& rng,
+                       u32bit bitsize)
+   {
+   set_sign(Positive);
+
+   if(bitsize == 0)
+      clear();
+   else
+      {
+      SecureVector<byte> array((bitsize + 7) / 8);
+      rng.randomize(array, array.size());
+      if(bitsize % 8)
+         array[0] &= 0xFF >> (8 - (bitsize % 8));
+      array[0] |= 0x80 >> ((bitsize % 8) ? (8 - bitsize % 8) : 0);
+      binary_decode(array, array.size());
+      }
+   }
+
+/*************************************************
+* Generate a random integer within given range   *
+*************************************************/
+BigInt random_integer(RandomNumberGenerator& rng,
+                      const BigInt& min, const BigInt& max)
+   {
+   BigInt range = max - min;
+
+   if(range <= 0)
+      throw Invalid_Argument("random_integer: invalid min/max values");
+
+   return (min + (BigInt(rng, range.bits() + 2) % range));
+   }
+
+/*************************************************
+* Generate a random safe prime                   *
+*************************************************/
+BigInt random_safe_prime(RandomNumberGenerator& rng, u32bit bits)
+   {
+   if(bits <= 64)
+      throw Invalid_Argument("random_safe_prime: Can't make a prime of " +
+                             to_string(bits) + " bits");
+
+   BigInt p;
+   do
+      p = (random_prime(rng, bits - 1) << 1) + 1;
+   while(!is_prime(p, rng));
+   return p;
+   }
+
+}
diff --git a/src/bigint/bigint.cpp b/src/bigint/bigint.cpp
new file mode 100644
index 000000000..e3c7931e6
--- /dev/null
+++ b/src/bigint/bigint.cpp
@@ -0,0 +1,367 @@
+/*************************************************
+* BigInt Base Source File                        *
+* (C) 1999-2008 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/bigint.h>
+#include <botan/mp_core.h>
+#include <botan/loadstor.h>
+#include <botan/parsing.h>
+#include <botan/util.h>
+
+namespace Botan {
+
+/*************************************************
+* Construct a BigInt from a regular number       *
+*************************************************/
+BigInt::BigInt(u64bit n)
+   {
+   set_sign(Positive);
+
+   if(n == 0)
+      return;
+
+   const u32bit limbs_needed = sizeof(u64bit) / sizeof(word);
+
+   reg.create(4*limbs_needed);
+   for(u32bit j = 0; j != limbs_needed; ++j)
+      reg[j] = ((n >> (j*MP_WORD_BITS)) & MP_WORD_MASK);
+   }
+
+/*************************************************
+* Construct a BigInt of the specified size       *
+*************************************************/
+BigInt::BigInt(Sign s, u32bit size)
+   {
+   reg.create(round_up(size, 8));
+   signedness = s;
+   }
+
+/*************************************************
+* Construct a BigInt from a "raw" BigInt         *
+*************************************************/
+BigInt::BigInt(const BigInt& b)
+   {
+   const u32bit b_words = b.sig_words();
+
+   if(b_words)
+      {
+      reg.create(round_up(b_words, 8));
+      reg.copy(b.data(), b_words);
+      set_sign(b.sign());
+      }
+   else
+      {
+      reg.create(2);
+      set_sign(Positive);
+      }
+   }
+
+/*************************************************
+* Construct a BigInt from a string               *
+*************************************************/
+BigInt::BigInt(const std::string& str)
+   {
+   Base base = Decimal;
+   u32bit markers = 0;
+   bool negative = false;
+   if(str.length() > 0 && str[0] == '-') { markers += 1; negative = true; }
+
+   if(str.length() > markers + 2 && str[markers    ] == '0' &&
+                                    str[markers + 1] == 'x')
+      { markers += 2; base = Hexadecimal; }
+   else if(str.length() > markers + 1 && str[markers] == '0')
+      { markers += 1; base = Octal; }
+
+   *this = decode(reinterpret_cast<const byte*>(str.data()) + markers,
+                  str.length() - markers, base);
+
+   if(negative) set_sign(Negative);
+   else         set_sign(Positive);
+   }
+
+/*************************************************
+* Construct a BigInt from an encoded BigInt      *
+*************************************************/
+BigInt::BigInt(const byte input[], u32bit length, Base base)
+   {
+   set_sign(Positive);
+   *this = decode(input, length, base);
+   }
+
+/*************************************************
+* Construct a BigInt from an encoded BigInt      *
+*************************************************/
+BigInt::BigInt(RandomNumberGenerator& rng, u32bit bits)
+   {
+   set_sign(Positive);
+   randomize(rng, bits);
+   }
+
+/*************************************************
+* Swap this BigInt with another                  *
+*************************************************/
+void BigInt::swap(BigInt& other)
+   {
+   reg.swap(other.reg);
+   std::swap(signedness, other.signedness);
+   }
+
+/*************************************************
+* Grow the internal storage                      *
+*************************************************/
+void BigInt::grow_reg(u32bit n)
+   {
+   reg.grow_to(round_up(size() + n, 8));
+   }
+
+/*************************************************
+* Grow the internal storage                      *
+*************************************************/
+void BigInt::grow_to(u32bit n)
+   {
+   if(n > size())
+      reg.grow_to(round_up(n, 8));
+   }
+
+/*************************************************
+* Comparison Function                            *
+*************************************************/
+s32bit BigInt::cmp(const BigInt& n, bool check_signs) const
+   {
+   if(check_signs)
+      {
+      if(n.is_positive() && this->is_negative()) return -1;
+      if(n.is_negative() && this->is_positive()) return 1;
+      if(n.is_negative() && this->is_negative())
+         return (-bigint_cmp(data(), sig_words(), n.data(), n.sig_words()));
+      }
+   return bigint_cmp(data(), sig_words(), n.data(), n.sig_words());
+   }
+
+/*************************************************
+* Convert this number to a u32bit, if possible   *
+*************************************************/
+u32bit BigInt::to_u32bit() const
+   {
+   if(is_negative())
+      throw Encoding_Error("BigInt::to_u32bit: Number is negative");
+   if(bits() >= 32)
+      throw Encoding_Error("BigInt::to_u32bit: Number is too big to convert");
+
+   u32bit out = 0;
+   for(u32bit j = 0; j != 4; ++j)
+      out = (out << 8) | byte_at(3-j);
+   return out;
+   }
+
+/*************************************************
+* Return byte n of this number                   *
+*************************************************/
+byte BigInt::byte_at(u32bit n) const
+   {
+   const u32bit WORD_BYTES = sizeof(word);
+   u32bit word_num = n / WORD_BYTES, byte_num = n % WORD_BYTES;
+   if(word_num >= size())
+      return 0;
+   else
+      return get_byte(WORD_BYTES - byte_num - 1, reg[word_num]);
+   }
+
+/*************************************************
+* Return bit n of this number                    *
+*************************************************/
+bool BigInt::get_bit(u32bit n) const
+   {
+   return ((word_at(n / MP_WORD_BITS) >> (n % MP_WORD_BITS)) & 1);
+   }
+
+/*************************************************
+* Return bits {offset...offset+length}           *
+*************************************************/
+u32bit BigInt::get_substring(u32bit offset, u32bit length) const
+   {
+   if(length > 32)
+      throw Invalid_Argument("BigInt::get_substring: Substring size too big");
+
+   u64bit piece = 0;
+   for(u32bit j = 0; j != 8; ++j)
+      piece = (piece << 8) | byte_at((offset / 8) + (7-j));
+
+   u64bit mask = (1 << length) - 1;
+   u32bit shift = (offset % 8);
+
+   return static_cast<u32bit>((piece >> shift) & mask);
+   }
+
+/*************************************************
+* Set bit number n                               *
+*************************************************/
+void BigInt::set_bit(u32bit n)
+   {
+   const u32bit which = n / MP_WORD_BITS;
+   const word mask = static_cast<word>(1) << (n % MP_WORD_BITS);
+   if(which >= size()) grow_to(which + 1);
+   reg[which] |= mask;
+   }
+
+/*************************************************
+* Clear bit number n                             *
+*************************************************/
+void BigInt::clear_bit(u32bit n)
+   {
+   const u32bit which = n / MP_WORD_BITS;
+   const word mask = static_cast<word>(1) << (n % MP_WORD_BITS);
+   if(which < size())
+      reg[which] &= ~mask;
+   }
+
+/*************************************************
+* Clear all but the lowest n bits                *
+*************************************************/
+void BigInt::mask_bits(u32bit n)
+   {
+   if(n == 0) { clear(); return; }
+   if(n >= bits()) return;
+
+   const u32bit top_word = n / MP_WORD_BITS;
+   const word mask = (static_cast<word>(1) << (n % MP_WORD_BITS)) - 1;
+
+   if(top_word < size())
+      for(u32bit j = top_word + 1; j != size(); ++j)
+         reg[j] = 0;
+
+   reg[top_word] &= mask;
+   }
+
+/*************************************************
+* Count how many bytes are being used            *
+*************************************************/
+u32bit BigInt::bytes() const
+   {
+   return (bits() + 7) / 8;
+   }
+
+/*************************************************
+* Count how many bits are being used             *
+*************************************************/
+u32bit BigInt::bits() const
+   {
+   if(sig_words() == 0)
+      return 0;
+
+   u32bit full_words = sig_words() - 1, top_bits = MP_WORD_BITS;
+   word top_word = word_at(full_words), mask = MP_WORD_TOP_BIT;
+
+   while(top_bits && ((top_word & mask) == 0))
+      { mask >>= 1; top_bits--; }
+
+   return (full_words * MP_WORD_BITS + top_bits);
+   }
+
+/*************************************************
+* Calcluate the size in a certain base           *
+*************************************************/
+u32bit BigInt::encoded_size(Base base) const
+   {
+   static const double LOG_2_BASE_10 = 0.30102999566;
+
+   if(base == Binary)
+      return bytes();
+   else if(base == Hexadecimal)
+      return 2*bytes();
+   else if(base == Octal)
+      return ((bits() + 2) / 3);
+   else if(base == Decimal)
+      return static_cast<u32bit>((bits() * LOG_2_BASE_10) + 1);
+   else
+      throw Invalid_Argument("Unknown base for BigInt encoding");
+   }
+
+/*************************************************
+* Set the sign                                   *
+*************************************************/
+void BigInt::set_sign(Sign s)
+   {
+   if(is_zero())
+      signedness = Positive;
+   else
+      signedness = s;
+   }
+
+/*************************************************
+* Reverse the value of the sign flag             *
+*************************************************/
+void BigInt::flip_sign()
+   {
+   set_sign(reverse_sign());
+   }
+
+/*************************************************
+* Return the opposite value of the current sign  *
+*************************************************/
+BigInt::Sign BigInt::reverse_sign() const
+   {
+   if(sign() == Positive)
+      return Negative;
+   return Positive;
+   }
+
+/*************************************************
+* Return the negation of this number             *
+*************************************************/
+BigInt BigInt::operator-() const
+   {
+   BigInt x = (*this);
+   x.flip_sign();
+   return x;
+   }
+
+/*************************************************
+* Return the absolute value of this number       *
+*************************************************/
+BigInt BigInt::abs() const
+   {
+   BigInt x = (*this);
+   x.set_sign(Positive);
+   return x;
+   }
+
+/*************************************************
+* Encode this number into bytes                  *
+*************************************************/
+void BigInt::binary_encode(byte output[]) const
+   {
+   const u32bit sig_bytes = bytes();
+   for(u32bit j = 0; j != sig_bytes; ++j)
+      output[sig_bytes-j-1] = byte_at(j);
+   }
+
+/*************************************************
+* Set this number to the value in buf            *
+*************************************************/
+void BigInt::binary_decode(const byte buf[], u32bit length)
+   {
+   const u32bit WORD_BYTES = sizeof(word);
+
+   reg.create(round_up((length / WORD_BYTES) + 1, 8));
+
+   for(u32bit j = 0; j != length / WORD_BYTES; ++j)
+      {
+      u32bit top = length - WORD_BYTES*j;
+      for(u32bit k = WORD_BYTES; k > 0; --k)
+         reg[j] = (reg[j] << 8) | buf[top - k];
+      }
+   for(u32bit j = 0; j != length % WORD_BYTES; ++j)
+      reg[length / WORD_BYTES] = (reg[length / WORD_BYTES] << 8) | buf[j];
+   }
+
+/*************************************************
+* Set this number to the value in buf            *
+*************************************************/
+void BigInt::binary_decode(const MemoryRegion<byte>& buf)
+   {
+   binary_decode(buf, buf.size());
+   }
+
+}
diff --git a/src/bigint/bigint.h b/src/bigint/bigint.h
new file mode 100644
index 000000000..b1286551e
--- /dev/null
+++ b/src/bigint/bigint.h
@@ -0,0 +1,181 @@
+/*************************************************
+* BigInt Header File                             *
+* (C) 1999-2008 Jack Lloyd                       *
+*************************************************/
+
+#ifndef BOTAN_BIGINT_H__
+#define BOTAN_BIGINT_H__
+
+#include <botan/rng.h>
+#include <botan/secmem.h>
+#include <botan/mp_types.h>
+#include <iosfwd>
+
+namespace Botan {
+
+/*************************************************
+* BigInt                                         *
+*************************************************/
+class BOTAN_DLL BigInt
+   {
+   public:
+      enum Base { Octal = 8, Decimal = 10, Hexadecimal = 16, Binary = 256 };
+      enum Sign { Negative = 0, Positive = 1 };
+      enum NumberType { Power2 };
+
+      struct DivideByZero : public Exception
+         { DivideByZero() : Exception("BigInt divide by zero") {} };
+
+      BigInt& operator+=(const BigInt&);
+      BigInt& operator-=(const BigInt&);
+
+      BigInt& operator*=(const BigInt&);
+      BigInt& operator/=(const BigInt&);
+      BigInt& operator%=(const BigInt&);
+      word    operator%=(word);
+      BigInt& operator<<=(u32bit);
+      BigInt& operator>>=(u32bit);
+
+      BigInt& operator++() { return (*this += 1); }
+      BigInt& operator--() { return (*this -= 1); }
+      BigInt  operator++(int) { BigInt x = (*this); ++(*this); return x; }
+      BigInt  operator--(int) { BigInt x = (*this); --(*this); return x; }
+
+      BigInt operator-() const;
+      bool operator !() const { return (!is_nonzero()); }
+
+      s32bit cmp(const BigInt&, bool = true) const;
+      bool is_even() const { return (get_bit(0) == 0); }
+      bool is_odd()  const { return (get_bit(0) == 1); }
+
+      bool is_zero() const
+         {
+         const u32bit sw = sig_words();
+
+         for(u32bit i = 0; i != sw; ++i)
+            if(reg[i])
+               return false;
+         return true;
+         }
+
+      bool is_nonzero() const { return (!is_zero()); }
+
+      void set_bit(u32bit);
+      void clear_bit(u32bit);
+      void mask_bits(u32bit);
+
+      bool get_bit(u32bit) const;
+      u32bit get_substring(u32bit, u32bit) const;
+      byte byte_at(u32bit) const;
+
+      // same as operator[], remove this
+      word word_at(u32bit n) const
+         { return ((n < size()) ? get_reg()[n] : 0); }
+
+      u32bit to_u32bit() const;
+
+      bool is_negative() const { return (sign() == Negative); }
+      bool is_positive() const { return (sign() == Positive); }
+      Sign sign() const { return (signedness); }
+      Sign reverse_sign() const;
+      void flip_sign();
+      void set_sign(Sign);
+      BigInt abs() const;
+
+      u32bit size() const { return get_reg().size(); }
+
+      u32bit sig_words() const
+         {
+         const word* x = reg.begin();
+         u32bit sig = reg.size();
+
+         while(sig && (x[sig-1] == 0))
+            sig--;
+         return sig;
+         }
+
+      u32bit bytes() const;
+      u32bit bits() const;
+
+      const word* data() const { return reg.begin(); }
+      SecureVector<word>& get_reg() { return reg; }
+      const SecureVector<word>& get_reg() const { return reg; }
+
+      void grow_reg(u32bit);
+      void grow_to(u32bit);
+
+      word& operator[](u32bit i) { return reg[i]; }
+      word operator[](u32bit i) const { return reg[i]; }
+      void clear() { get_reg().clear(); }
+
+      void randomize(RandomNumberGenerator& rng, u32bit n);
+
+      void binary_encode(byte[]) const;
+      void binary_decode(const byte[], u32bit);
+      void binary_decode(const MemoryRegion<byte>&);
+      u32bit encoded_size(Base = Binary) const;
+
+      static SecureVector<byte> encode(const BigInt&, Base = Binary);
+      static void encode(byte[], const BigInt&, Base = Binary);
+      static BigInt decode(const byte[], u32bit, Base = Binary);
+      static BigInt decode(const MemoryRegion<byte>&, Base = Binary);
+      static SecureVector<byte> encode_1363(const BigInt&, u32bit);
+
+      void swap(BigInt&);
+
+      BigInt() { signedness = Positive; }
+      BigInt(u64bit);
+      BigInt(const BigInt&);
+      BigInt(const std::string&);
+      BigInt(const byte[], u32bit, Base = Binary);
+      BigInt(RandomNumberGenerator& rng, u32bit bits);
+      BigInt(Sign, u32bit);
+      BigInt(NumberType, u32bit);
+   private:
+      SecureVector<word> reg;
+      Sign signedness;
+   };
+
+/*************************************************
+* Arithmetic Operators                           *
+*************************************************/
+BigInt BOTAN_DLL operator+(const BigInt&, const BigInt&);
+BigInt BOTAN_DLL operator-(const BigInt&, const BigInt&);
+BigInt BOTAN_DLL operator*(const BigInt&, const BigInt&);
+BigInt BOTAN_DLL operator/(const BigInt&, const BigInt&);
+BigInt BOTAN_DLL operator%(const BigInt&, const BigInt&);
+word   BOTAN_DLL operator%(const BigInt&, word);
+BigInt BOTAN_DLL operator<<(const BigInt&, u32bit);
+BigInt BOTAN_DLL operator>>(const BigInt&, u32bit);
+
+/*************************************************
+* Comparison Operators                           *
+*************************************************/
+inline bool operator==(const BigInt& a, const BigInt& b)
+   { return (a.cmp(b) == 0); }
+inline bool operator!=(const BigInt& a, const BigInt& b)
+   { return (a.cmp(b) != 0); }
+inline bool operator<=(const BigInt& a, const BigInt& b)
+   { return (a.cmp(b) <= 0); }
+inline bool operator>=(const BigInt& a, const BigInt& b)
+   { return (a.cmp(b) >= 0); }
+inline bool operator<(const BigInt& a, const BigInt& b)
+   { return (a.cmp(b) < 0); }
+inline bool operator>(const BigInt& a, const BigInt& b)
+   { return (a.cmp(b) > 0); }
+
+/*************************************************
+* I/O Operators                                  *
+*************************************************/
+BOTAN_DLL std::ostream& operator<<(std::ostream&, const BigInt&);
+BOTAN_DLL std::istream& operator>>(std::istream&, BigInt&);
+
+}
+
+namespace std {
+
+inline void swap(Botan::BigInt& a, Botan::BigInt& b) { a.swap(b); }
+
+}
+
+#endif
diff --git a/src/bigint/blinding.cpp b/src/bigint/blinding.cpp
new file mode 100644
index 000000000..740904d10
--- /dev/null
+++ b/src/bigint/blinding.cpp
@@ -0,0 +1,47 @@
+/*************************************************
+* Blinder Source File                            *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/blinding.h>
+#include <botan/numthry.h>
+
+namespace Botan {
+
+/*************************************************
+* Blinder Constructor                            *
+*************************************************/
+Blinder::Blinder(const BigInt& e, const BigInt& d, const BigInt& n)
+   {
+   if(e < 1 || d < 1 || n < 1)
+      throw Invalid_Argument("Blinder: Arguments too small");
+
+   reducer = Modular_Reducer(n);
+   this->e = e;
+   this->d = d;
+   }
+
+/*************************************************
+* Blind a number                                 *
+*************************************************/
+BigInt Blinder::blind(const BigInt& i) const
+   {
+   if(!reducer.initialized())
+      return i;
+
+   e = reducer.square(e);
+   d = reducer.square(d);
+   return reducer.multiply(i, e);
+   }
+
+/*************************************************
+* Unblind a number                               *
+*************************************************/
+BigInt Blinder::unblind(const BigInt& i) const
+   {
+   if(!reducer.initialized())
+      return i;
+   return reducer.multiply(i, d);
+   }
+
+}
diff --git a/src/bigint/blinding.h b/src/bigint/blinding.h
new file mode 100644
index 000000000..958686fb1
--- /dev/null
+++ b/src/bigint/blinding.h
@@ -0,0 +1,32 @@
+/*************************************************
+* Blinder Header File                            *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#ifndef BOTAN_BLINDER_H__
+#define BOTAN_BLINDER_H__
+
+#include <botan/bigint.h>
+#include <botan/reducer.h>
+
+namespace Botan {
+
+/*************************************************
+* Blinding Function Object                       *
+*************************************************/
+class BOTAN_DLL Blinder
+   {
+   public:
+      BigInt blind(const BigInt&) const;
+      BigInt unblind(const BigInt&) const;
+
+      Blinder() {}
+      Blinder(const BigInt&, const BigInt&, const BigInt&);
+   private:
+      Modular_Reducer reducer;
+      mutable BigInt e, d;
+   };
+
+}
+
+#endif
diff --git a/src/bigint/def_powm.h b/src/bigint/def_powm.h
new file mode 100644
index 000000000..c91ff002c
--- /dev/null
+++ b/src/bigint/def_powm.h
@@ -0,0 +1,62 @@
+/*************************************************
+* Modular Exponentiation Header File             *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#ifndef BOTAN_DEFAULT_MODEXP_H__
+#define BOTAN_DEFAULT_MODEXP_H__
+
+#include <botan/pow_mod.h>
+#include <botan/reducer.h>
+#include <vector>
+
+namespace Botan {
+
+/*************************************************
+* Fixed Window Exponentiator                     *
+*************************************************/
+class BOTAN_DLL Fixed_Window_Exponentiator : public Modular_Exponentiator
+   {
+   public:
+      void set_exponent(const BigInt&);
+      void set_base(const BigInt&);
+      BigInt execute() const;
+
+      Modular_Exponentiator* copy() const
+         { return new Fixed_Window_Exponentiator(*this); }
+
+      Fixed_Window_Exponentiator(const BigInt&, Power_Mod::Usage_Hints);
+   private:
+      Modular_Reducer reducer;
+      BigInt exp;
+      u32bit window_bits;
+      std::vector<BigInt> g;
+      Power_Mod::Usage_Hints hints;
+   };
+
+/*************************************************
+* Montgomery Exponentiator                       *
+*************************************************/
+class BOTAN_DLL Montgomery_Exponentiator : public Modular_Exponentiator
+   {
+   public:
+      void set_exponent(const BigInt&);
+      void set_base(const BigInt&);
+      BigInt execute() const;
+
+      Modular_Exponentiator* copy() const
+         { return new Montgomery_Exponentiator(*this); }
+
+      Montgomery_Exponentiator(const BigInt&, Power_Mod::Usage_Hints);
+   private:
+      BigInt exp, modulus;
+      BigInt R2, R_mod;
+      std::vector<BigInt> g;
+      word mod_prime;
+      u32bit mod_words, exp_bits, window_bits;
+      Power_Mod::Usage_Hints hints;
+   };
+
+}
+
+#endif
diff --git a/src/bigint/divide.cpp b/src/bigint/divide.cpp
new file mode 100644
index 000000000..003a06d83
--- /dev/null
+++ b/src/bigint/divide.cpp
@@ -0,0 +1,104 @@
+/*************************************************
+* Division Algorithm Source File                 *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/numthry.h>
+#include <botan/mp_core.h>
+
+namespace Botan {
+
+namespace {
+
+/*************************************************
+* Handle signed operands, if necessary           *
+*************************************************/
+void sign_fixup(const BigInt& x, const BigInt& y, BigInt& q, BigInt& r)
+   {
+   if(x.sign() == BigInt::Negative)
+      {
+      q.flip_sign();
+      if(r.is_nonzero()) { --q; r = y.abs() - r; }
+      }
+   if(y.sign() == BigInt::Negative)
+      q.flip_sign();
+   }
+
+}
+
+/*************************************************
+* Solve x = q * y + r                            *
+*************************************************/
+void divide(const BigInt& x, const BigInt& y_arg, BigInt& q, BigInt& r)
+   {
+   if(y_arg.is_zero())
+      throw BigInt::DivideByZero();
+
+   BigInt y = y_arg;
+   const u32bit y_words = y.sig_words();
+   r = x;
+
+   r.set_sign(BigInt::Positive);
+   y.set_sign(BigInt::Positive);
+
+   s32bit compare = r.cmp(y);
+
+   if(compare < 0)
+      q = 0;
+   else if(compare ==  0)
+      {
+      q = 1;
+      r = 0;
+      }
+   else
+      {
+      u32bit shifts = 0;
+      word y_top = y[y.sig_words()-1];
+      while(y_top < MP_WORD_TOP_BIT) { y_top <<= 1; ++shifts; }
+      y <<= shifts;
+      r <<= shifts;
+
+      const u32bit n = r.sig_words() - 1, t = y_words - 1;
+
+      q.get_reg().create(n - t + 1);
+      if(n <= t)
+         {
+         while(r > y) { r -= y; ++q; }
+         r >>= shifts;
+         sign_fixup(x, y_arg, q, r);
+         return;
+         }
+
+      BigInt temp = y << (MP_WORD_BITS * (n-t));
+
+      while(r >= temp) { r -= temp; ++q[n-t]; }
+
+      for(u32bit j = n; j != t; --j)
+         {
+         const word x_j0  = r.word_at(j);
+         const word x_j1 = r.word_at(j-1);
+         const word y_t  = y.word_at(t);
+
+         if(x_j0 == y_t)
+            q[j-t-1] = MP_WORD_MAX;
+         else
+            q[j-t-1] = bigint_divop(x_j0, x_j1, y_t);
+
+         while(bigint_divcore(q[j-t-1], y_t, y.word_at(t-1),
+                              x_j0, x_j1, r.word_at(j-2)))
+            --q[j-t-1];
+
+         r -= (q[j-t-1] * y) << (MP_WORD_BITS * (j-t-1));
+         if(r.is_negative())
+            {
+            r += y << (MP_WORD_BITS * (j-t-1));
+            --q[j-t-1];
+            }
+         }
+      r >>= shifts;
+      }
+
+   sign_fixup(x, y_arg, q, r);
+   }
+
+}
diff --git a/src/bigint/dsa_gen.cpp b/src/bigint/dsa_gen.cpp
new file mode 100644
index 000000000..baaba66ee
--- /dev/null
+++ b/src/bigint/dsa_gen.cpp
@@ -0,0 +1,134 @@
+/*************************************************
+* DSA Parameter Generation Source File           *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/dl_group.h>
+#include <botan/numthry.h>
+#include <botan/lookup.h>
+#include <botan/parsing.h>
+#include <algorithm>
+#include <memory>
+
+namespace Botan {
+
+namespace {
+
+/*************************************************
+* Check if this size is allowed by FIPS 186-3    *
+*************************************************/
+bool fips186_3_valid_size(u32bit pbits, u32bit qbits)
+   {
+   if(qbits == 160)
+      return (pbits == 512 || pbits == 768 || pbits == 1024);
+
+   if(qbits == 224)
+      return (pbits == 2048);
+
+   if(qbits == 256)
+      return (pbits == 2048 || pbits == 3072);
+
+   return false;
+   }
+
+}
+
+/*************************************************
+* Attempt DSA prime generation with given seed   *
+*************************************************/
+bool DL_Group::generate_dsa_primes(RandomNumberGenerator& rng,
+                                   BigInt& p, BigInt& q,
+                                   u32bit pbits, u32bit qbits,
+                                   const MemoryRegion<byte>& seed_c)
+   {
+   if(!fips186_3_valid_size(pbits, qbits))
+      throw Invalid_Argument(
+         "FIPS 186-3 does not allow DSA domain parameters of " +
+         to_string(pbits) + "/" + to_string(qbits) + " bits long");
+
+   if(qbits == 224)
+      throw Invalid_Argument(
+         "DSA parameter generation with a q of 224 bits not supported");
+
+   if(seed_c.size() * 8 < qbits)
+      throw Invalid_Argument(
+         "Generating a DSA parameter set with a " + to_string(qbits) +
+         "long q requires a seed at least as many bits long");
+
+   std::auto_ptr<HashFunction> hash(get_hash("SHA-" + to_string(qbits)));
+
+   const u32bit HASH_SIZE = hash->OUTPUT_LENGTH;
+
+   class Seed
+      {
+      public:
+         Seed(const MemoryRegion<byte>& s) : seed(s) {}
+
+         operator MemoryRegion<byte>& () { return seed; }
+
+         Seed& operator++()
+            {
+            for(u32bit j = seed.size(); j > 0; --j)
+               if(++seed[j-1])
+                  break;
+            return (*this);
+            }
+      private:
+         SecureVector<byte> seed;
+      };
+
+   Seed seed(seed_c);
+
+   q.binary_decode(hash->process(seed));
+   q.set_bit(qbits-1);
+   q.set_bit(0);
+
+   if(!is_prime(q, rng))
+      return false;
+
+   const u32bit n = (pbits-1) / (HASH_SIZE * 8),
+                b = (pbits-1) % (HASH_SIZE * 8);
+
+   BigInt X;
+   SecureVector<byte> V(HASH_SIZE * (n+1));
+
+   for(u32bit j = 0; j != 4096; ++j)
+      {
+      for(u32bit k = 0; k <= n; ++k)
+         {
+         ++seed;
+         hash->update(seed);
+         hash->final(V + HASH_SIZE * (n-k));
+         }
+
+      X.binary_decode(V + (HASH_SIZE - 1 - b/8),
+                      V.size() - (HASH_SIZE - 1 - b/8));
+      X.set_bit(pbits-1);
+
+      p = X - (X % (2*q) - 1);
+
+      if(p.bits() == pbits && is_prime(p, rng))
+         return true;
+      }
+   return false;
+   }
+
+/*************************************************
+* Generate DSA Primes                            *
+*************************************************/
+SecureVector<byte> DL_Group::generate_dsa_primes(RandomNumberGenerator& rng,
+                                                 BigInt& p, BigInt& q,
+                                                 u32bit pbits, u32bit qbits)
+   {
+   SecureVector<byte> seed(qbits/8);
+
+   while(true)
+      {
+      rng.randomize(seed, seed.size());
+
+      if(generate_dsa_primes(rng, p, q, pbits, qbits, seed))
+         return seed;
+      }
+   }
+
+}
diff --git a/src/bigint/jacobi.cpp b/src/bigint/jacobi.cpp
new file mode 100644
index 000000000..57c78508a
--- /dev/null
+++ b/src/bigint/jacobi.cpp
@@ -0,0 +1,51 @@
+/*************************************************
+* Jacobi Function Source File                    *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/numthry.h>
+
+namespace Botan {
+
+/*************************************************
+* Calculate the Jacobi symbol                    *
+*************************************************/
+s32bit jacobi(const BigInt& a, const BigInt& n)
+   {
+   if(a.is_negative())
+      throw Invalid_Argument("jacobi: first argument must be non-negative");
+   if(n.is_even() || n < 2)
+      throw Invalid_Argument("jacobi: second argument must be odd and > 1");
+
+   BigInt x = a, y = n;
+   s32bit J = 1;
+
+   while(y > 1)
+      {
+      x %= y;
+      if(x > y / 2)
+         {
+         x = y - x;
+         if(y % 4 == 3)
+            J = -J;
+         }
+      if(x.is_zero())
+         return 0;
+
+      u32bit shifts = low_zero_bits(x);
+      x >>= shifts;
+      if(shifts % 2)
+         {
+         word y_mod_8 = y % 8;
+         if(y_mod_8 == 3 || y_mod_8 == 5)
+            J = -J;
+         }
+
+      if(x % 4 == 3 && y % 4 == 3)
+         J = -J;
+      std::swap(x, y);
+      }
+   return J;
+   }
+
+}
diff --git a/src/bigint/make_prm.cpp b/src/bigint/make_prm.cpp
new file mode 100644
index 000000000..dc26a0400
--- /dev/null
+++ b/src/bigint/make_prm.cpp
@@ -0,0 +1,79 @@
+/*************************************************
+* Prime Generation Source File                   *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/numthry.h>
+#include <botan/parsing.h>
+#include <algorithm>
+
+namespace Botan {
+
+/*************************************************
+* Generate a random prime                        *
+*************************************************/
+BigInt random_prime(RandomNumberGenerator& rng,
+                    u32bit bits, const BigInt& coprime,
+                    u32bit equiv, u32bit modulo)
+   {
+   if(bits <= 1)
+      throw Invalid_Argument("random_prime: Can't make a prime of " +
+                             to_string(bits) + " bits");
+   else if(bits == 2)
+      return ((rng.next_byte() % 1) ? 2 : 3);
+   else if(bits == 3)
+      return ((rng.next_byte() % 1) ? 5 : 7);
+   else if(bits == 4)
+      return ((rng.next_byte() % 1) ? 11 : 13);
+
+   if(coprime <= 0)
+      throw Invalid_Argument("random_prime: coprime must be > 0");
+   if(modulo % 2 == 1 || modulo == 0)
+      throw Invalid_Argument("random_prime: Invalid modulo value");
+   if(equiv >= modulo || equiv % 2 == 0)
+      throw Invalid_Argument("random_prime: equiv must be < modulo, and odd");
+
+   while(true)
+      {
+      BigInt p(rng, bits);
+      p.set_bit(bits - 2);
+      p.set_bit(0);
+
+      if(p % modulo != equiv)
+         p += (modulo - p % modulo) + equiv;
+
+      const u32bit sieve_size = std::min(bits / 2, PRIME_TABLE_SIZE);
+      SecureVector<u32bit> sieve(sieve_size);
+
+      for(u32bit j = 0; j != sieve.size(); ++j)
+         sieve[j] = p % PRIMES[j];
+
+      u32bit counter = 0;
+      while(true)
+         {
+         if(counter == 4096 || p.bits() > bits)
+            break;
+
+         bool passes_sieve = true;
+         ++counter;
+         p += modulo;
+
+         if(p.bits() > bits)
+            break;
+
+         for(u32bit j = 0; j != sieve.size(); ++j)
+            {
+            sieve[j] = (sieve[j] + modulo) % PRIMES[j];
+            if(sieve[j] == 0)
+               passes_sieve = false;
+            }
+
+         if(!passes_sieve || gcd(p - 1, coprime) != 1)
+            continue;
+         if(passes_mr_tests(rng, p))
+            return p;
+         }
+      }
+   }
+
+}
diff --git a/src/bigint/modinfo.txt b/src/bigint/modinfo.txt
new file mode 100644
index 000000000..571c3550f
--- /dev/null
+++ b/src/bigint/modinfo.txt
@@ -0,0 +1,43 @@
+realname "BigInt"
+
+load_on auto
+
+define BIGINT
+
+<add>
+big_code.cpp
+big_io.cpp
+big_ops2.cpp
+big_ops3.cpp
+big_rand.cpp
+bigint.cpp
+bigint.h
+blinding.cpp
+blinding.h
+def_powm.h
+divide.cpp
+dsa_gen.cpp
+jacobi.cpp
+make_prm.cpp
+mp_asm.cpp
+mp_comba.cpp
+mp_core.h
+mp_karat.cpp
+mp_misc.cpp
+mp_monty.cpp
+mp_mulop.cpp
+mp_numth.cpp
+mp_shift.cpp
+numthry.cpp
+numthry.h
+pow_mod.cpp
+pow_mod.h
+powm_fw.cpp
+powm_mnt.cpp
+reducer.cpp
+reducer.h
+ressol.cpp
+mp_types.h
+mp_asm.h
+mp_asmi.h
+</add>
diff --git a/src/bigint/mp_amd64/modinfo.txt b/src/bigint/mp_amd64/xxxinfo.txt
index a042a3976..a042a3976 100644
--- a/src/bigint/mp_amd64/modinfo.txt
+++ b/src/bigint/mp_amd64/xxxinfo.txt
diff --git a/src/bigint/mp_asm.cpp b/src/bigint/mp_asm.cpp
new file mode 100644
index 000000000..e5d1fe0d6
--- /dev/null
+++ b/src/bigint/mp_asm.cpp
@@ -0,0 +1,177 @@
+/*************************************************
+* Lowest Level MPI Algorithms Source File        *
+* (C) 1999-2008 Jack Lloyd                       *
+*     2006 Luca Piccarreta                       *
+*************************************************/
+
+#include <botan/mp_asm.h>
+#include <botan/mp_asmi.h>
+#include <botan/mp_core.h>
+#include <botan/mem_ops.h>
+
+namespace Botan {
+
+extern "C" {
+
+/*************************************************
+* Two Operand Addition, No Carry                 *
+*************************************************/
+word bigint_add2_nc(word x[], u32bit x_size, const word y[], u32bit y_size)
+   {
+   word carry = 0;
+
+   const u32bit blocks = y_size - (y_size % 8);
+
+   for(u32bit j = 0; j != blocks; j += 8)
+      carry = word8_add2(x + j, y + j, carry);
+
+   for(u32bit j = blocks; j != y_size; ++j)
+      x[j] = word_add(x[j], y[j], &carry);
+
+   if(!carry)
+      return 0;
+
+   for(u32bit j = y_size; j != x_size; ++j)
+      if(++x[j])
+         return 0;
+
+   return 1;
+   }
+
+/*************************************************
+* Three Operand Addition, No Carry               *
+*************************************************/
+word bigint_add3_nc(word z[], const word x[], u32bit x_size,
+                              const word y[], u32bit y_size)
+   {
+   if(x_size < y_size)
+      { return bigint_add3_nc(z, y, y_size, x, x_size); }
+
+   word carry = 0;
+
+   const u32bit blocks = y_size - (y_size % 8);
+
+   for(u32bit j = 0; j != blocks; j += 8)
+      carry = word8_add3(z + j, x + j, y + j, carry);
+
+   for(u32bit j = blocks; j != y_size; ++j)
+      z[j] = word_add(x[j], y[j], &carry);
+
+   for(u32bit j = y_size; j != x_size; ++j)
+      {
+      word x_j = x[j] + carry;
+      if(carry && x_j)
+         carry = 0;
+      z[j] = x_j;
+      }
+
+   return carry;
+   }
+
+/*************************************************
+* Two Operand Addition                           *
+*************************************************/
+void bigint_add2(word x[], u32bit x_size, const word y[], u32bit y_size)
+   {
+   if(bigint_add2_nc(x, x_size, y, y_size))
+      ++x[x_size];
+   }
+
+/*************************************************
+* Three Operand Addition                         *
+*************************************************/
+void bigint_add3(word z[], const word x[], u32bit x_size,
+                           const word y[], u32bit y_size)
+   {
+   if(bigint_add3_nc(z, x, x_size, y, y_size))
+      ++z[(x_size > y_size ? x_size : y_size)];
+   }
+
+/*************************************************
+* Two Operand Subtraction                        *
+*************************************************/
+void bigint_sub2(word x[], u32bit x_size, const word y[], u32bit y_size)
+   {
+   word carry = 0;
+
+   const u32bit blocks = y_size - (y_size % 8);
+
+   for(u32bit j = 0; j != blocks; j += 8)
+      carry = word8_sub2(x + j, y + j, carry);
+
+   for(u32bit j = blocks; j != y_size; ++j)
+      x[j] = word_sub(x[j], y[j], &carry);
+
+   if(!carry) return;
+
+   for(u32bit j = y_size; j != x_size; ++j)
+      {
+      --x[j];
+      if(x[j] != MP_WORD_MAX) return;
+      }
+   }
+
+/*************************************************
+* Three Operand Subtraction                      *
+*************************************************/
+void bigint_sub3(word z[], const word x[], u32bit x_size,
+                           const word y[], u32bit y_size)
+   {
+   word carry = 0;
+
+   const u32bit blocks = y_size - (y_size % 8);
+
+   for(u32bit j = 0; j != blocks; j += 8)
+      carry = word8_sub3(z + j, x + j, y + j, carry);
+
+   for(u32bit j = blocks; j != y_size; ++j)
+      z[j] = word_sub(x[j], y[j], &carry);
+
+   for(u32bit j = y_size; j != x_size; ++j)
+      {
+      word x_j = x[j] - carry;
+      if(carry && x_j != MP_WORD_MAX)
+         carry = 0;
+      z[j] = x_j;
+      }
+   }
+
+/*************************************************
+* Two Operand Linear Multiply                    *
+*************************************************/
+void bigint_linmul2(word x[], u32bit x_size, word y)
+   {
+   const u32bit blocks = x_size - (x_size % 8);
+
+   word carry = 0;
+
+   for(u32bit j = 0; j != blocks; j += 8)
+      carry = word8_linmul2(x + j, y, carry);
+
+   for(u32bit j = blocks; j != x_size; ++j)
+      x[j] = word_madd2(x[j], y, &carry);
+
+   x[x_size] = carry;
+   }
+
+/*************************************************
+* Three Operand Linear Multiply                  *
+*************************************************/
+void bigint_linmul3(word z[], const word x[], u32bit x_size, word y)
+   {
+   const u32bit blocks = x_size - (x_size % 8);
+
+   word carry = 0;
+
+   for(u32bit j = 0; j != blocks; j += 8)
+      carry = word8_linmul3(z + j, x + j, y, carry);
+
+   for(u32bit j = blocks; j != x_size; ++j)
+      z[j] = word_madd2(x[j], y, &carry);
+
+   z[x_size] = carry;
+   }
+
+}
+
+}
diff --git a/src/bigint/mp_asm.h b/src/bigint/mp_asm.h
new file mode 100644
index 000000000..e62a57110
--- /dev/null
+++ b/src/bigint/mp_asm.h
@@ -0,0 +1,52 @@
+/*************************************************
+* Lowest Level MPI Algorithms Header File        *
+* (C) 1999-2008 Jack Lloyd                       *
+*     2006 Luca Piccarreta                       *
+*************************************************/
+
+#ifndef BOTAN_MP_ASM_H__
+#define BOTAN_MP_ASM_H__
+
+#include <botan/mp_types.h>
+
+#if (BOTAN_MP_WORD_BITS == 8)
+  typedef Botan::u16bit dword;
+#elif (BOTAN_MP_WORD_BITS == 16)
+  typedef Botan::u32bit dword;
+#elif (BOTAN_MP_WORD_BITS == 32)
+  typedef Botan::u64bit dword;
+#elif (BOTAN_MP_WORD_BITS == 64)
+  #error BOTAN_MP_WORD_BITS can be 64 only with assembly support
+#else
+  #error BOTAN_MP_WORD_BITS must be 8, 16, 32, or 64
+#endif
+
+namespace Botan {
+
+extern "C" {
+
+/*************************************************
+* Word Multiply/Add                              *
+*************************************************/
+inline word word_madd2(word a, word b, word* c)
+   {
+   dword z = (dword)a * b + *c;
+   *c = (word)(z >> BOTAN_MP_WORD_BITS);
+   return (word)z;
+   }
+
+/*************************************************
+* Word Multiply/Add                              *
+*************************************************/
+inline word word_madd3(word a, word b, word c, word* d)
+   {
+   dword z = (dword)a * b + c + *d;
+   *d = (word)(z >> BOTAN_MP_WORD_BITS);
+   return (word)z;
+   }
+
+}
+
+}
+
+#endif
diff --git a/src/bigint/mp_asm64/modinfo.txt b/src/bigint/mp_asm64/xxxinfo.txt
index a9e5d53da..a9e5d53da 100644
--- a/src/bigint/mp_asm64/modinfo.txt
+++ b/src/bigint/mp_asm64/xxxinfo.txt
diff --git a/src/bigint/mp_asmi.h b/src/bigint/mp_asmi.h
new file mode 100644
index 000000000..d15295154
--- /dev/null
+++ b/src/bigint/mp_asmi.h
@@ -0,0 +1,189 @@
+/*************************************************
+* Lowest Level MPI Algorithms Header File        *
+* (C) 1999-2008 Jack Lloyd                       *
+*     2006 Luca Piccarreta                       *
+*************************************************/
+
+#ifndef BOTAN_MP_ASM_INTERNAL_H__
+#define BOTAN_MP_ASM_INTERNAL_H__
+
+#include <botan/mp_asm.h>
+
+namespace Botan {
+
+extern "C" {
+
+/*************************************************
+* Word Addition                                  *
+*************************************************/
+inline word word_add(word x, word y, word* carry)
+   {
+   word z = x + y;
+   word c1 = (z < x);
+   z += *carry;
+   *carry = c1 | (z < *carry);
+   return z;
+   }
+
+/*************************************************
+* Eight Word Block Addition, Two Argument        *
+*************************************************/
+inline word word8_add2(word x[8], const word y[8], word carry)
+   {
+   x[0] = word_add(x[0], y[0], &carry);
+   x[1] = word_add(x[1], y[1], &carry);
+   x[2] = word_add(x[2], y[2], &carry);
+   x[3] = word_add(x[3], y[3], &carry);
+   x[4] = word_add(x[4], y[4], &carry);
+   x[5] = word_add(x[5], y[5], &carry);
+   x[6] = word_add(x[6], y[6], &carry);
+   x[7] = word_add(x[7], y[7], &carry);
+   return carry;
+   }
+
+/*************************************************
+* Eight Word Block Addition, Three Argument      *
+*************************************************/
+inline word word8_add3(word z[8], const word x[8],
+                       const word y[8], word carry)
+   {
+   z[0] = word_add(x[0], y[0], &carry);
+   z[1] = word_add(x[1], y[1], &carry);
+   z[2] = word_add(x[2], y[2], &carry);
+   z[3] = word_add(x[3], y[3], &carry);
+   z[4] = word_add(x[4], y[4], &carry);
+   z[5] = word_add(x[5], y[5], &carry);
+   z[6] = word_add(x[6], y[6], &carry);
+   z[7] = word_add(x[7], y[7], &carry);
+   return carry;
+   }
+
+/*************************************************
+* Word Subtraction                               *
+*************************************************/
+inline word word_sub(word x, word y, word* carry)
+   {
+   word t0 = x - y;
+   word c1 = (t0 > x);
+   word z = t0 - *carry;
+   *carry = c1 | (z > t0);
+   return z;
+   }
+
+/*************************************************
+* Eight Word Block Subtraction, Two Argument     *
+*************************************************/
+inline word word8_sub2(word x[4], const word y[4], word carry)
+   {
+   x[0] = word_sub(x[0], y[0], &carry);
+   x[1] = word_sub(x[1], y[1], &carry);
+   x[2] = word_sub(x[2], y[2], &carry);
+   x[3] = word_sub(x[3], y[3], &carry);
+   x[4] = word_sub(x[4], y[4], &carry);
+   x[5] = word_sub(x[5], y[5], &carry);
+   x[6] = word_sub(x[6], y[6], &carry);
+   x[7] = word_sub(x[7], y[7], &carry);
+   return carry;
+   }
+
+/*************************************************
+* Eight Word Block Subtraction, Three Argument   *
+*************************************************/
+inline word word8_sub3(word z[8], const word x[8],
+                       const word y[8], word carry)
+   {
+   z[0] = word_sub(x[0], y[0], &carry);
+   z[1] = word_sub(x[1], y[1], &carry);
+   z[2] = word_sub(x[2], y[2], &carry);
+   z[3] = word_sub(x[3], y[3], &carry);
+   z[4] = word_sub(x[4], y[4], &carry);
+   z[5] = word_sub(x[5], y[5], &carry);
+   z[6] = word_sub(x[6], y[6], &carry);
+   z[7] = word_sub(x[7], y[7], &carry);
+   return carry;
+   }
+
+/*************************************************
+* Eight Word Block Linear Multiplication         *
+*************************************************/
+inline word word8_linmul2(word x[4], word y, word carry)
+   {
+   x[0] = word_madd2(x[0], y, &carry);
+   x[1] = word_madd2(x[1], y, &carry);
+   x[2] = word_madd2(x[2], y, &carry);
+   x[3] = word_madd2(x[3], y, &carry);
+   x[4] = word_madd2(x[4], y, &carry);
+   x[5] = word_madd2(x[5], y, &carry);
+   x[6] = word_madd2(x[6], y, &carry);
+   x[7] = word_madd2(x[7], y, &carry);
+   return carry;
+   }
+
+/*************************************************
+* Eight Word Block Linear Multiplication         *
+*************************************************/
+inline word word8_linmul3(word z[8], const word x[8], word y, word carry)
+   {
+   z[0] = word_madd2(x[0], y, &carry);
+   z[1] = word_madd2(x[1], y, &carry);
+   z[2] = word_madd2(x[2], y, &carry);
+   z[3] = word_madd2(x[3], y, &carry);
+   z[4] = word_madd2(x[4], y, &carry);
+   z[5] = word_madd2(x[5], y, &carry);
+   z[6] = word_madd2(x[6], y, &carry);
+   z[7] = word_madd2(x[7], y, &carry);
+   return carry;
+   }
+
+/*************************************************
+* Eight Word Block Multiply/Add                  *
+*************************************************/
+inline word word8_madd3(word z[8], const word x[8], word y, word carry)
+   {
+   z[0] = word_madd3(x[0], y, z[0], &carry);
+   z[1] = word_madd3(x[1], y, z[1], &carry);
+   z[2] = word_madd3(x[2], y, z[2], &carry);
+   z[3] = word_madd3(x[3], y, z[3], &carry);
+   z[4] = word_madd3(x[4], y, z[4], &carry);
+   z[5] = word_madd3(x[5], y, z[5], &carry);
+   z[6] = word_madd3(x[6], y, z[6], &carry);
+   z[7] = word_madd3(x[7], y, z[7], &carry);
+   return carry;
+   }
+
+/*************************************************
+* Multiply-Add Accumulator                       *
+*************************************************/
+inline void word3_muladd(word* w2, word* w1, word* w0, word a, word b)
+   {
+   word carry = *w0;
+   *w0 = word_madd2(a, b, &carry);
+   *w1 += carry;
+   *w2 += (*w1 < carry) ? 1 : 0;
+   }
+
+/*************************************************
+* Multiply-Add Accumulator                       *
+*************************************************/
+inline void word3_muladd_2(word* w2, word* w1, word* w0, word a, word b)
+   {
+   word carry = 0;
+   a = word_madd2(a, b, &carry);
+   b = carry;
+
+   word top = (b >> (BOTAN_MP_WORD_BITS-1));
+   b <<= 1;
+   b |= (a >> (BOTAN_MP_WORD_BITS-1));
+   a <<= 1;
+
+   carry = 0;
+   *w0 = word_add(*w0, a, &carry);
+   *w1 = word_add(*w1, b, &carry);
+   *w2 = word_add(*w2, top, &carry);
+   }
+
+}
+
+}
+
+#endif
diff --git a/src/bigint/mp_comba.cpp b/src/bigint/mp_comba.cpp
new file mode 100644
index 000000000..c7a9c964c
--- /dev/null
+++ b/src/bigint/mp_comba.cpp
@@ -0,0 +1,918 @@
+/*************************************************
+* Comba Multiplication and Squaring Source File  *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/mp_core.h>
+#include <botan/mp_asmi.h>
+
+namespace Botan {
+
+extern "C" {
+
+/*************************************************
+* Comba 4x4 Squaring                             *
+*************************************************/
+void bigint_comba_sqr4(word z[8], const word x[4])
+   {
+   word w2 = 0, w1 = 0, w0 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], x[0]);
+   z[0] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[1]);
+   z[1] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[2]);
+   word3_muladd(&w2, &w1, &w0, x[1], x[1]);
+   z[2] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[3]);
+   word3_muladd_2(&w2, &w1, &w0, x[1], x[2]);
+   z[3] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[1], x[3]);
+   word3_muladd(&w2, &w1, &w0, x[2], x[2]);
+   z[4] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[2], x[3]);
+   z[5] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[3], x[3]);
+   z[6] = w0;
+   z[7] = w1;
+   }
+
+/*************************************************
+* Comba 4x4 Multiplication                       *
+*************************************************/
+void bigint_comba_mul4(word z[8], const word x[4], const word y[4])
+   {
+   word w2 = 0, w1 = 0, w0 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[0]);
+   z[0] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[0]);
+   z[1] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[0]);
+   z[2] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[0]);
+   z[3] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[1], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[1]);
+   z[4] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[2], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[2]);
+   z[5] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[3], y[3]);
+   z[6] = w0;
+   z[7] = w1;
+   }
+
+/*************************************************
+* Comba 6x6 Squaring                             *
+*************************************************/
+void bigint_comba_sqr6(word z[12], const word x[6])
+   {
+   word w2 = 0, w1 = 0, w0 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], x[0]);
+   z[0] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[1]);
+   z[1] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[2]);
+   word3_muladd(&w2, &w1, &w0, x[1], x[1]);
+   z[2] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[3]);
+   word3_muladd_2(&w2, &w1, &w0, x[1], x[2]);
+   z[3] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[4]);
+   word3_muladd_2(&w2, &w1, &w0, x[1], x[3]);
+   word3_muladd(&w2, &w1, &w0, x[2], x[2]);
+   z[4] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[5]);
+   word3_muladd_2(&w2, &w1, &w0, x[1], x[4]);
+   word3_muladd_2(&w2, &w1, &w0, x[2], x[3]);
+   z[5] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[1], x[5]);
+   word3_muladd_2(&w2, &w1, &w0, x[2], x[4]);
+   word3_muladd(&w2, &w1, &w0, x[3], x[3]);
+   z[6] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[2], x[5]);
+   word3_muladd_2(&w2, &w1, &w0, x[3], x[4]);
+   z[7] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[3], x[5]);
+   word3_muladd(&w2, &w1, &w0, x[4], x[4]);
+   z[8] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[4], x[5]);
+   z[9] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[5], x[5]);
+   z[10] = w0;
+   z[11] = w1;
+   }
+
+/*************************************************
+* Comba 6x6 Multiplication                       *
+*************************************************/
+void bigint_comba_mul6(word z[12], const word x[6], const word y[6])
+   {
+   word w2 = 0, w1 = 0, w0 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[0]);
+   z[0] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[0]);
+   z[1] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[0]);
+   z[2] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[0]);
+   z[3] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[4]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[4], y[0]);
+   z[4] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[5]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[4]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[4], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[5], y[0]);
+   z[5] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[1], y[5]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[4]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[4], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[5], y[1]);
+   z[6] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[2], y[5]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[4]);
+   word3_muladd(&w2, &w1, &w0, x[4], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[5], y[2]);
+   z[7] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[3], y[5]);
+   word3_muladd(&w2, &w1, &w0, x[4], y[4]);
+   word3_muladd(&w2, &w1, &w0, x[5], y[3]);
+   z[8] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[4], y[5]);
+   word3_muladd(&w2, &w1, &w0, x[5], y[4]);
+   z[9] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[5], y[5]);
+   z[10] = w0;
+   z[11] = w1;
+   }
+
+/*************************************************
+* Comba 8x8 Squaring                             *
+*************************************************/
+void bigint_comba_sqr8(word z[16], const word x[8])
+   {
+   word w2 = 0, w1 = 0, w0 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], x[0]);
+   z[0] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[1]);
+   z[1] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[2]);
+   word3_muladd(&w2, &w1, &w0, x[1], x[1]);
+   z[2] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[3]);
+   word3_muladd_2(&w2, &w1, &w0, x[1], x[2]);
+   z[3] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[4]);
+   word3_muladd_2(&w2, &w1, &w0, x[1], x[3]);
+   word3_muladd(&w2, &w1, &w0, x[2], x[2]);
+   z[4] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[5]);
+   word3_muladd_2(&w2, &w1, &w0, x[1], x[4]);
+   word3_muladd_2(&w2, &w1, &w0, x[2], x[3]);
+   z[5] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[6]);
+   word3_muladd_2(&w2, &w1, &w0, x[1], x[5]);
+   word3_muladd_2(&w2, &w1, &w0, x[2], x[4]);
+   word3_muladd(&w2, &w1, &w0, x[3], x[3]);
+   z[6] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[0], x[7]);
+   word3_muladd_2(&w2, &w1, &w0, x[1], x[6]);
+   word3_muladd_2(&w2, &w1, &w0, x[2], x[5]);
+   word3_muladd_2(&w2, &w1, &w0, x[3], x[4]);
+   z[7] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[1], x[7]);
+   word3_muladd_2(&w2, &w1, &w0, x[2], x[6]);
+   word3_muladd_2(&w2, &w1, &w0, x[3], x[5]);
+   word3_muladd(&w2, &w1, &w0, x[4], x[4]);
+   z[8] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[2], x[7]);
+   word3_muladd_2(&w2, &w1, &w0, x[3], x[6]);
+   word3_muladd_2(&w2, &w1, &w0, x[4], x[5]);
+   z[9] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[3], x[7]);
+   word3_muladd_2(&w2, &w1, &w0, x[4], x[6]);
+   word3_muladd(&w2, &w1, &w0, x[5], x[5]);
+   z[10] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[4], x[7]);
+   word3_muladd_2(&w2, &w1, &w0, x[5], x[6]);
+   z[11] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[5], x[7]);
+   word3_muladd(&w2, &w1, &w0, x[6], x[6]);
+   z[12] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[6], x[7]);
+   z[13] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[7], x[7]);
+   z[14] = w0;
+   z[15] = w1;
+   }
+
+/*************************************************
+* Comba 8x8 Multiplication                       *
+*************************************************/
+void bigint_comba_mul8(word z[16], const word x[8], const word y[8])
+   {
+   word w2 = 0, w1 = 0, w0 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[0]);
+   z[0] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[0]);
+   z[1] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[0]);
+   z[2] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[0]);
+   z[3] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[4]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[4], y[0]);
+   z[4] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[5]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[4]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[4], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[5], y[0]);
+   z[5] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[6]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[5]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[4]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[4], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[5], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[6], y[0]);
+   z[6] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[0], y[7]);
+   word3_muladd(&w2, &w1, &w0, x[1], y[6]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[5]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[4]);
+   word3_muladd(&w2, &w1, &w0, x[4], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[5], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[6], y[1]);
+   word3_muladd(&w2, &w1, &w0, x[7], y[0]);
+   z[7] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[1], y[7]);
+   word3_muladd(&w2, &w1, &w0, x[2], y[6]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[5]);
+   word3_muladd(&w2, &w1, &w0, x[4], y[4]);
+   word3_muladd(&w2, &w1, &w0, x[5], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[6], y[2]);
+   word3_muladd(&w2, &w1, &w0, x[7], y[1]);
+   z[8] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[2], y[7]);
+   word3_muladd(&w2, &w1, &w0, x[3], y[6]);
+   word3_muladd(&w2, &w1, &w0, x[4], y[5]);
+   word3_muladd(&w2, &w1, &w0, x[5], y[4]);
+   word3_muladd(&w2, &w1, &w0, x[6], y[3]);
+   word3_muladd(&w2, &w1, &w0, x[7], y[2]);
+   z[9] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[3], y[7]);
+   word3_muladd(&w2, &w1, &w0, x[4], y[6]);
+   word3_muladd(&w2, &w1, &w0, x[5], y[5]);
+   word3_muladd(&w2, &w1, &w0, x[6], y[4]);
+   word3_muladd(&w2, &w1, &w0, x[7], y[3]);
+   z[10] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[4], y[7]);
+   word3_muladd(&w2, &w1, &w0, x[5], y[6]);
+   word3_muladd(&w2, &w1, &w0, x[6], y[5]);
+   word3_muladd(&w2, &w1, &w0, x[7], y[4]);
+   z[11] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[5], y[7]);
+   word3_muladd(&w2, &w1, &w0, x[6], y[6]);
+   word3_muladd(&w2, &w1, &w0, x[7], y[5]);
+   z[12] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[6], y[7]);
+   word3_muladd(&w2, &w1, &w0, x[7], y[6]);
+   z[13] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[7], y[7]);
+   z[14] = w0;
+   z[15] = w1;
+   }
+
+/*************************************************
+* Comba 16x16 Squaring                           *
+*************************************************/
+void bigint_comba_sqr16(word z[32], const word x[16])
+   {
+   word w2 = 0, w1 = 0, w0 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], x[ 0]);
+   z[ 0] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[ 1]);
+   z[ 1] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], x[ 1]);
+   z[ 2] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[ 3]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[ 2]);
+   z[ 3] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[ 4]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], x[ 2]);
+   z[ 4] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[ 5]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[ 4]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[ 3]);
+   z[ 5] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[ 6]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[ 5]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], x[ 3]);
+   z[ 6] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[ 7]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[ 6]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[ 5]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 3], x[ 4]);
+   z[ 7] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[ 8]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[ 7]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[ 6]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 3], x[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], x[ 4]);
+   z[ 8] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[ 9]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[ 8]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[ 7]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 3], x[ 6]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 4], x[ 5]);
+   z[ 9] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[10]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[ 9]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[ 8]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 3], x[ 7]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 4], x[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], x[ 5]);
+   z[10] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[11]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[10]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[ 9]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 3], x[ 8]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 4], x[ 7]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 5], x[ 6]);
+   z[11] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[12]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[11]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[10]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 3], x[ 9]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 4], x[ 8]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 5], x[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], x[ 6]);
+   z[12] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[13]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[12]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[11]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 3], x[10]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 4], x[ 9]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 5], x[ 8]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 6], x[ 7]);
+   z[13] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[14]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[13]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[12]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 3], x[11]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 4], x[10]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 5], x[ 9]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 6], x[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], x[ 7]);
+   z[14] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 0], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[14]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[13]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 3], x[12]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 4], x[11]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 5], x[10]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 6], x[ 9]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 7], x[ 8]);
+   z[15] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 1], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[14]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 3], x[13]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 4], x[12]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 5], x[11]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 6], x[10]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 7], x[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], x[ 8]);
+   z[16] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 2], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 3], x[14]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 4], x[13]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 5], x[12]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 6], x[11]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 7], x[10]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 8], x[ 9]);
+   z[17] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 3], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 4], x[14]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 5], x[13]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 6], x[12]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 7], x[11]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 8], x[10]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], x[ 9]);
+   z[18] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 4], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 5], x[14]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 6], x[13]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 7], x[12]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 8], x[11]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 9], x[10]);
+   z[19] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 5], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 6], x[14]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 7], x[13]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 8], x[12]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 9], x[11]);
+   word3_muladd(&w2, &w1, &w0, x[10], x[10]);
+   z[20] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 6], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 7], x[14]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 8], x[13]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 9], x[12]);
+   word3_muladd_2(&w2, &w1, &w0, x[10], x[11]);
+   z[21] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 7], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 8], x[14]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 9], x[13]);
+   word3_muladd_2(&w2, &w1, &w0, x[10], x[12]);
+   word3_muladd(&w2, &w1, &w0, x[11], x[11]);
+   z[22] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 8], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[ 9], x[14]);
+   word3_muladd_2(&w2, &w1, &w0, x[10], x[13]);
+   word3_muladd_2(&w2, &w1, &w0, x[11], x[12]);
+   z[23] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[ 9], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[10], x[14]);
+   word3_muladd_2(&w2, &w1, &w0, x[11], x[13]);
+   word3_muladd(&w2, &w1, &w0, x[12], x[12]);
+   z[24] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[10], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[11], x[14]);
+   word3_muladd_2(&w2, &w1, &w0, x[12], x[13]);
+   z[25] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[11], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[12], x[14]);
+   word3_muladd(&w2, &w1, &w0, x[13], x[13]);
+   z[26] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[12], x[15]);
+   word3_muladd_2(&w2, &w1, &w0, x[13], x[14]);
+   z[27] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[13], x[15]);
+   word3_muladd(&w2, &w1, &w0, x[14], x[14]);
+   z[28] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd_2(&w2, &w1, &w0, x[14], x[15]);
+   z[29] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[15], x[15]);
+   z[30] = w0;
+   z[31] = w1;
+   }
+
+/*************************************************
+* Comba 16x16 Multiplication                     *
+*************************************************/
+void bigint_comba_mul16(word z[32], const word x[16], const word y[16])
+   {
+   word w2 = 0, w1 = 0, w0 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[ 0]);
+   z[0] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[ 0]);
+   z[1] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[ 0]);
+   z[2] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[ 0]);
+   z[3] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[ 0]);
+   z[4] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[ 0]);
+   z[5] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[ 0]);
+   z[6] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[ 0]);
+   z[7] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[ 0]);
+   z[8] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[ 0]);
+   z[9] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[ 0]);
+   z[10] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[ 0]);
+   z[11] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[ 0]);
+   z[12] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[ 0]);
+   z[13] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[ 0]);
+   z[14] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 0], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[ 1]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[ 0]);
+   z[15] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 1], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[ 2]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[ 1]);
+   z[16] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 2], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[ 3]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[ 2]);
+   z[17] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 3], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[ 4]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[ 3]);
+   z[18] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 4], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[ 5]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[ 4]);
+   z[19] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 5], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[ 6]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[ 5]);
+   z[20] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 6], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[ 7]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[ 6]);
+   z[21] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 7], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[ 8]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[ 7]);
+   z[22] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 8], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[ 9]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[ 8]);
+   z[23] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[ 9], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[10], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[10]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[ 9]);
+   z[24] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[10], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[11], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[11]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[10]);
+   z[25] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[11], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[12], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[12]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[11]);
+   z[26] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[12], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[13], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[13]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[12]);
+   z[27] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[13], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[14], y[14]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[13]);
+   z[28] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[14], y[15]);
+   word3_muladd(&w2, &w1, &w0, x[15], y[14]);
+   z[29] = w0; w0 = w1; w1 = w2; w2 = 0;
+
+   word3_muladd(&w2, &w1, &w0, x[15], y[15]);
+   z[30] = w0;
+   z[31] = w1;
+   }
+
+}
+
+}
diff --git a/src/bigint/mp_core.h b/src/bigint/mp_core.h
new file mode 100644
index 000000000..92949cd83
--- /dev/null
+++ b/src/bigint/mp_core.h
@@ -0,0 +1,96 @@
+/*************************************************
+* MPI Algorithms Header File                     *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#ifndef BOTAN_MP_CORE_H__
+#define BOTAN_MP_CORE_H__
+
+#include <botan/mp_types.h>
+
+namespace Botan {
+
+/*************************************************
+* The size of the word type, in bits             *
+*************************************************/
+const u32bit MP_WORD_BITS = BOTAN_MP_WORD_BITS;
+
+extern "C" {
+
+/*************************************************
+* Addition/Subtraction Operations                *
+*************************************************/
+void bigint_add2(word[], u32bit, const word[], u32bit);
+void bigint_add3(word[], const word[], u32bit, const word[], u32bit);
+
+word bigint_add2_nc(word[], u32bit, const word[], u32bit);
+word bigint_add3_nc(word[], const word[], u32bit, const word[], u32bit);
+
+void bigint_sub2(word[], u32bit, const word[], u32bit);
+void bigint_sub3(word[], const word[], u32bit, const word[], u32bit);
+
+/*************************************************
+* Shift Operations                               *
+*************************************************/
+void bigint_shl1(word[], u32bit, u32bit, u32bit);
+void bigint_shl2(word[], const word[], u32bit, u32bit, u32bit);
+void bigint_shr1(word[], u32bit, u32bit, u32bit);
+void bigint_shr2(word[], const word[], u32bit, u32bit, u32bit);
+
+/*************************************************
+* Simple O(N^2) Multiplication and Squaring      *
+*************************************************/
+void bigint_simple_mul(word z[], const word x[], u32bit x_size,
+                       const word y[], u32bit y_size);
+void bigint_simple_sqr(word z[], const word x[], u32bit x_size);
+
+/*************************************************
+* Linear Multiply                                *
+*************************************************/
+void bigint_linmul2(word[], u32bit, word);
+void bigint_linmul3(word[], const word[], u32bit, word);
+void bigint_linmul_add(word[], u32bit, const word[], u32bit, word);
+
+/*************************************************
+* Montgomery Reduction                           *
+*************************************************/
+void bigint_monty_redc(word[], u32bit, const word[], u32bit, word);
+
+/*************************************************
+* Misc Utility Operations                        *
+*************************************************/
+u32bit bigint_divcore(word, word, word, word, word, word);
+s32bit bigint_cmp(const word[], u32bit, const word[], u32bit);
+word bigint_divop(word, word, word);
+word bigint_modop(word, word, word);
+void bigint_wordmul(word, word, word*, word*);
+
+/*************************************************
+* Comba Multiplication / Squaring                *
+*************************************************/
+void bigint_comba_mul4(word[8], const word[4], const word[4]);
+void bigint_comba_mul6(word[12], const word[6], const word[6]);
+void bigint_comba_mul8(word[16], const word[8], const word[8]);
+void bigint_comba_mul16(word[32], const word[16], const word[16]);
+
+void bigint_comba_sqr4(word[8], const word[4]);
+void bigint_comba_sqr6(word[12], const word[6]);
+void bigint_comba_sqr8(word[16], const word[8]);
+void bigint_comba_sqr8(word[32], const word[16]);
+void bigint_comba_sqr16(word[64], const word[32]);
+
+}
+
+/*************************************************
+* High Level Multiplication/Squaring Interfaces  *
+*************************************************/
+void bigint_mul(word[], u32bit, word[],
+                const word[], u32bit, u32bit,
+                const word[], u32bit, u32bit);
+
+void bigint_sqr(word[], u32bit, word[],
+                const word[], u32bit, u32bit);
+
+}
+
+#endif
diff --git a/src/bigint/mp_ia32/modinfo.txt b/src/bigint/mp_ia32/xxxinfo.txt
index cf4959250..cf4959250 100644
--- a/src/bigint/mp_ia32/modinfo.txt
+++ b/src/bigint/mp_ia32/xxxinfo.txt
diff --git a/src/bigint/mp_ia32_msvc/modinfo.txt b/src/bigint/mp_ia32_msvc/xxxinfo.txt
index 36d9d0290..36d9d0290 100644
--- a/src/bigint/mp_ia32_msvc/modinfo.txt
+++ b/src/bigint/mp_ia32_msvc/xxxinfo.txt
diff --git a/src/bigint/mp_karat.cpp b/src/bigint/mp_karat.cpp
new file mode 100644
index 000000000..15b0551fd
--- /dev/null
+++ b/src/bigint/mp_karat.cpp
@@ -0,0 +1,334 @@
+/*************************************************
+* Karatsuba Multiplication/Squaring Source File  *
+* (C) 1999-2008 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/mp_core.h>
+#include <botan/mem_ops.h>
+#include <botan/mp_asmi.h>
+
+namespace Botan {
+
+namespace {
+
+/*************************************************
+* Karatsuba Multiplication Operation             *
+*************************************************/
+void karatsuba_mul(word z[], const word x[], const word y[], u32bit N,
+                   word workspace[])
+   {
+   if(N == 6)
+      bigint_comba_mul6(z, x, y);
+   else if(N == 8)
+      bigint_comba_mul8(z, x, y);
+   else if(N == 16)
+      bigint_comba_mul16(z, x, y);
+   else if(N < BOTAN_KARAT_MUL_THRESHOLD || N % 2)
+      bigint_simple_mul(z, x, N, y, N);
+   else
+      {
+      const u32bit N2 = N / 2;
+
+      const word* x0 = x;
+      const word* x1 = x + N2;
+      const word* y0 = y;
+      const word* y1 = y + N2;
+      word* z0 = z;
+      word* z1 = z + N;
+
+      const s32bit cmp0 = bigint_cmp(x0, N2, x1, N2);
+      const s32bit cmp1 = bigint_cmp(y1, N2, y0, N2);
+
+      clear_mem(workspace, 2*N);
+
+      if(cmp0 && cmp1)
+         {
+         if(cmp0 > 0)
+            bigint_sub3(z0, x0, N2, x1, N2);
+         else
+            bigint_sub3(z0, x1, N2, x0, N2);
+
+         if(cmp1 > 0)
+            bigint_sub3(z1, y1, N2, y0, N2);
+         else
+            bigint_sub3(z1, y0, N2, y1, N2);
+
+         karatsuba_mul(workspace, z0, z1, N2, workspace+N);
+         }
+
+      karatsuba_mul(z0, x0, y0, N2, workspace+N);
+      karatsuba_mul(z1, x1, y1, N2, workspace+N);
+
+      const u32bit blocks_of_8 = N - (N % 8);
+
+      word carry = 0;
+
+      for(u32bit j = 0; j != blocks_of_8; j += 8)
+         carry = word8_add3(workspace + N + j, z0 + j, z1 + j, carry);
+
+      for(u32bit j = blocks_of_8; j != N; ++j)
+         workspace[N + j] = word_add(z0[j], z1[j], &carry);
+
+      word carry2 = 0;
+
+      for(u32bit j = 0; j != blocks_of_8; j += 8)
+         carry2 = word8_add2(z + N2 + j, workspace + N + j, carry2);
+
+      for(u32bit j = blocks_of_8; j != N; ++j)
+         z[N2 + j] = word_add(z[N2 + j], workspace[N + j], &carry2);
+
+      z[N + N2] = word_add(z[N + N2], carry2, &carry);
+
+      if(carry)
+         for(u32bit j = 1; j != N2; ++j)
+            if(++z[N + N2 + j])
+               break;
+
+      if((cmp0 == cmp1) || (cmp0 == 0) || (cmp1 == 0))
+         bigint_add2(z + N2, 2*N-N2, workspace, N);
+      else
+         bigint_sub2(z + N2, 2*N-N2, workspace, N);
+      }
+   }
+
+/*************************************************
+* Karatsuba Squaring Operation                   *
+*************************************************/
+void karatsuba_sqr(word z[], const word x[], u32bit N, word workspace[])
+   {
+   if(N == 6)
+      bigint_comba_sqr6(z, x);
+   else if(N == 8)
+      bigint_comba_sqr8(z, x);
+   else if(N == 16)
+      bigint_comba_sqr16(z, x);
+   else if(N < BOTAN_KARAT_SQR_THRESHOLD || N % 2)
+      bigint_simple_sqr(z, x, N);
+   else
+      {
+      const u32bit N2 = N / 2;
+
+      const word* x0 = x;
+      const word* x1 = x + N2;
+      word* z0 = z;
+      word* z1 = z + N;
+
+      const s32bit cmp = bigint_cmp(x0, N2, x1, N2);
+
+      clear_mem(workspace, 2*N);
+
+      if(cmp)
+         {
+         if(cmp > 0)
+            bigint_sub3(z0, x0, N2, x1, N2);
+         else
+            bigint_sub3(z0, x1, N2, x0, N2);
+
+         karatsuba_sqr(workspace, z0, N2, workspace+N);
+         }
+
+      karatsuba_sqr(z0, x0, N2, workspace+N);
+      karatsuba_sqr(z1, x1, N2, workspace+N);
+
+      const u32bit blocks_of_8 = N - (N % 8);
+
+      word carry = 0;
+
+      for(u32bit j = 0; j != blocks_of_8; j += 8)
+         carry = word8_add3(workspace + N + j, z0 + j, z1 + j, carry);
+
+      for(u32bit j = blocks_of_8; j != N; ++j)
+         workspace[N + j] = word_add(z0[j], z1[j], &carry);
+
+      word carry2 = 0;
+
+      for(u32bit j = 0; j != blocks_of_8; j += 8)
+         carry2 = word8_add2(z + N2 + j, workspace + N + j, carry2);
+
+      for(u32bit j = blocks_of_8; j != N; ++j)
+         z[N2 + j] = word_add(z[N2 + j], workspace[N + j], &carry2);
+
+      z[N + N2] = word_add(z[N + N2], carry2, &carry);
+
+      if(carry)
+         for(u32bit j = 1; j != N2; ++j)
+            if(++z[N + N2 + j])
+               break;
+
+      if(cmp == 0)
+         bigint_add2(z + N2, 2*N-N2, workspace, N);
+      else
+         bigint_sub2(z + N2, 2*N-N2, workspace, N);
+      }
+   }
+
+/*************************************************
+* Pick a good size for the Karatsuba multiply    *
+*************************************************/
+u32bit karatsuba_size(u32bit z_size,
+                      u32bit x_size, u32bit x_sw,
+                      u32bit y_size, u32bit y_sw)
+   {
+   if(x_sw > x_size || x_sw > y_size || y_sw > x_size || y_sw > y_size)
+      return 0;
+
+   if(((x_size == x_sw) && (x_size % 2)) ||
+      ((y_size == y_sw) && (y_size % 2)))
+      return 0;
+
+   const u32bit start = (x_sw > y_sw) ? x_sw : y_sw;
+   const u32bit end = (x_size < y_size) ? x_size : y_size;
+
+   if(start == end)
+      {
+      if(start % 2)
+         return 0;
+      return start;
+      }
+
+   for(u32bit j = start; j <= end; ++j)
+      {
+      if(j % 2)
+         continue;
+
+      if(2*j > z_size)
+         return 0;
+
+      if(x_sw <= j && j <= x_size && y_sw <= j && j <= y_size)
+         {
+         if(j % 4 == 2 &&
+            (j+2) <= x_size && (j+2) <= y_size && 2*(j+2) <= z_size)
+            return j+2;
+         return j;
+         }
+      }
+
+   return 0;
+   }
+
+/*************************************************
+* Pick a good size for the Karatsuba squaring    *
+*************************************************/
+u32bit karatsuba_size(u32bit z_size, u32bit x_size, u32bit x_sw)
+   {
+   if(x_sw == x_size)
+      {
+      if(x_sw % 2)
+         return 0;
+      return x_sw;
+      }
+
+   for(u32bit j = x_sw; j <= x_size; ++j)
+      {
+      if(j % 2)
+         continue;
+
+      if(2*j > z_size)
+         return 0;
+
+      if(j % 4 == 2 && (j+2) <= x_size && 2*(j+2) <= z_size)
+         return j+2;
+      return j;
+      }
+
+   return 0;
+   }
+
+}
+
+/*************************************************
+* Multiplication Algorithm Dispatcher            *
+*************************************************/
+void bigint_mul(word z[], u32bit z_size, word workspace[],
+                const word x[], u32bit x_size, u32bit x_sw,
+                const word y[], u32bit y_size, u32bit y_sw)
+   {
+   if(x_sw == 1)
+      {
+      bigint_linmul3(z, y, y_sw, x[0]);
+      }
+   else if(y_sw == 1)
+      {
+      bigint_linmul3(z, x, x_sw, y[0]);
+      }
+   else if(x_sw <= 4 && x_size >= 4 &&
+           y_sw <= 4 && y_size >= 4 && z_size >= 8)
+      {
+      bigint_comba_mul4(z, x, y);
+      }
+   else if(x_sw <= 6 && x_size >= 6 &&
+           y_sw <= 6 && y_size >= 6 && z_size >= 12)
+      {
+      bigint_comba_mul6(z, x, y);
+      }
+   else if(x_sw <= 8 && x_size >= 8 &&
+           y_sw <= 8 && y_size >= 8 && z_size >= 16)
+      {
+      bigint_comba_mul8(z, x, y);
+      }
+   else if(x_sw <= 16 && x_size >= 16 &&
+           y_sw <= 16 && y_size >= 16 && z_size >= 32)
+      {
+      bigint_comba_mul16(z, x, y);
+      }
+   else if(x_sw < BOTAN_KARAT_MUL_THRESHOLD || y_sw < BOTAN_KARAT_MUL_THRESHOLD)
+      bigint_simple_mul(z, x, x_sw, y, y_sw);
+   else
+      {
+      const u32bit N = karatsuba_size(z_size, x_size, x_sw, y_size, y_sw);
+
+      if(N)
+         {
+         clear_mem(workspace, 2*N);
+         karatsuba_mul(z, x, y, N, workspace);
+         }
+      else
+         bigint_simple_mul(z, x, x_sw, y, y_sw);
+      }
+   }
+
+/*************************************************
+* Squaring Algorithm Dispatcher                  *
+*************************************************/
+void bigint_sqr(word z[], u32bit z_size, word workspace[],
+                const word x[], u32bit x_size, u32bit x_sw)
+   {
+   if(x_sw == 1)
+      {
+      bigint_linmul3(z, x, x_sw, x[0]);
+      }
+   else if(x_sw <= 4 && x_size >= 4 && z_size >= 8)
+      {
+      bigint_comba_sqr4(z, x);
+      }
+   else if(x_sw <= 6 && x_size >= 6 && z_size >= 12)
+      {
+      bigint_comba_sqr6(z, x);
+      }
+   else if(x_sw <= 8 && x_size >= 8 && z_size >= 16)
+      {
+      bigint_comba_sqr8(z, x);
+      }
+   else if(x_sw <= 16 && x_size >= 16 && z_size >= 32)
+      {
+      bigint_comba_sqr16(z, x);
+      }
+   else if(x_size < BOTAN_KARAT_SQR_THRESHOLD)
+      {
+      bigint_simple_sqr(z, x, x_sw);
+      }
+   else
+      {
+      const u32bit N = karatsuba_size(z_size, x_size, x_sw);
+
+      if(N)
+         {
+         clear_mem(workspace, 2*N);
+         karatsuba_sqr(z, x, N, workspace);
+         }
+      else
+         bigint_simple_sqr(z, x, x_sw);
+      }
+   }
+
+}
diff --git a/src/bigint/mp_misc.cpp b/src/bigint/mp_misc.cpp
new file mode 100644
index 000000000..db9c8cda0
--- /dev/null
+++ b/src/bigint/mp_misc.cpp
@@ -0,0 +1,92 @@
+/*************************************************
+* MP Misc Functions Source File                  *
+* (C) 1999-2008 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/mp_core.h>
+#include <botan/mp_asm.h>
+
+namespace Botan {
+
+extern "C" {
+
+/*************************************************
+* Core Division Operation                        *
+*************************************************/
+u32bit bigint_divcore(word q, word y1, word y2,
+                      word x1, word x2, word x3)
+   {
+   word y0 = 0;
+   y2 = word_madd2(q, y2, &y0);
+   y1 = word_madd2(q, y1, &y0);
+
+   if(y0 > x1) return 1;
+   if(y0 < x1) return 0;
+   if(y1 > x2)  return 1;
+   if(y1 < x2)  return 0;
+   if(y2 > x3)  return 1;
+   if(y2 < x3)  return 0;
+   return 0;
+   }
+
+/*************************************************
+* Compare two MP integers                        *
+*************************************************/
+s32bit bigint_cmp(const word x[], u32bit x_size,
+                  const word y[], u32bit y_size)
+   {
+   if(x_size < y_size) { return (-bigint_cmp(y, y_size, x, x_size)); }
+
+   while(x_size > y_size)
+      {
+      if(x[x_size-1])
+         return 1;
+      x_size--;
+      }
+   for(u32bit j = x_size; j > 0; --j)
+      {
+      if(x[j-1] > y[j-1]) return 1;
+      if(x[j-1] < y[j-1]) return -1;
+      }
+   return 0;
+   }
+
+/*************************************************
+* Do a 2-word/1-word Division                    *
+*************************************************/
+word bigint_divop(word n1, word n0, word d)
+   {
+   word high = n1 % d, quotient = 0;
+
+   for(u32bit j = 0; j != MP_WORD_BITS; ++j)
+      {
+      word high_top_bit = (high & MP_WORD_TOP_BIT);
+
+      high <<= 1;
+      high |= (n0 >> (MP_WORD_BITS-1-j)) & 1;
+      quotient <<= 1;
+
+      if(high_top_bit || high >= d)
+         {
+         high -= d;
+         quotient |= 1;
+         }
+      }
+
+   return quotient;
+   }
+
+/*************************************************
+* Do a 2-word/1-word Modulo                      *
+*************************************************/
+word bigint_modop(word n1, word n0, word d)
+   {
+   word z = bigint_divop(n1, n0, d);
+   word dummy = 0;
+   z = word_madd2(z, d, &dummy);
+   return (n0-z);
+   }
+
+}
+
+}
diff --git a/src/bigint/mp_monty.cpp b/src/bigint/mp_monty.cpp
new file mode 100644
index 000000000..c162bfd4f
--- /dev/null
+++ b/src/bigint/mp_monty.cpp
@@ -0,0 +1,76 @@
+/*************************************************
+* Montgomery Reduction Source File               *
+* (C) 1999-2008 Jack Lloyd                       *
+*     2006 Luca Piccarreta                       *
+*************************************************/
+
+#include <botan/mp_core.h>
+#include <botan/mp_asm.h>
+#include <botan/mp_asmi.h>
+
+namespace Botan {
+
+extern "C" {
+
+/*************************************************
+* Montgomery Reduction Algorithm                 *
+*************************************************/
+void bigint_monty_redc(word z[], u32bit z_size,
+                       const word x[], u32bit x_size, word u)
+   {
+   const u32bit blocks_of_8 = x_size - (x_size % 8);
+
+   for(u32bit i = 0; i != x_size; ++i)
+      {
+      word* z_i = z + i;
+
+      const word y = z_i[0] * u;
+
+      word carry = 0;
+
+      for(u32bit j = 0; j != blocks_of_8; j += 8)
+         carry = word8_madd3(z_i + j, x + j, y, carry);
+
+      for(u32bit j = blocks_of_8; j != x_size; ++j)
+         z_i[j] = word_madd3(x[j], y, z_i[j], &carry);
+
+      word z_sum = z_i[x_size] + carry;
+      carry = (z_sum < z_i[x_size]);
+      z_i[x_size] = z_sum;
+
+      for(u32bit j = x_size + 1; carry && j != z_size - i; ++j)
+         {
+         ++z_i[j];
+         carry = !z_i[j];
+         }
+      }
+
+   // Check if z[x_size...x_size+1] >= x[0...x_size] using bigint_cmp (inlined)
+   if(!z[x_size + x_size])
+      {
+      for(u32bit i = x_size; i > 0; --i)
+         {
+         if(z[x_size + i - 1] > x[i-1])
+            break;
+
+         if(z[x_size + i - 1] < x[i-1])
+            return;
+         }
+      }
+
+   // If the compare above is true, subtract using bigint_sub2 (inlined)
+   word carry = 0;
+
+   for(u32bit i = 0; i != blocks_of_8; i += 8)
+      carry = word8_sub2(z + x_size + i, x + i, carry);
+
+   for(u32bit i = blocks_of_8; i != x_size; ++i)
+      z[x_size + i] = word_sub(z[x_size + i], x[i], &carry);
+
+   if(carry)
+      --z[x_size+x_size];
+   }
+
+}
+
+}
diff --git a/src/bigint/mp_mulop.cpp b/src/bigint/mp_mulop.cpp
new file mode 100644
index 000000000..3ab28d306
--- /dev/null
+++ b/src/bigint/mp_mulop.cpp
@@ -0,0 +1,72 @@
+/*************************************************
+* Simple O(N^2) Multiplication and Squaring      *
+* (C) 1999-2008 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/mp_asm.h>
+#include <botan/mp_asmi.h>
+#include <botan/mp_core.h>
+#include <botan/mem_ops.h>
+
+namespace Botan {
+
+extern "C" {
+
+/*************************************************
+* Simple O(N^2) Multiplication                   *
+*************************************************/
+void bigint_simple_mul(word z[], const word x[], u32bit x_size,
+                                 const word y[], u32bit y_size)
+   {
+   const u32bit x_size_8 = x_size - (x_size % 8);
+
+   clear_mem(z, x_size + y_size);
+
+   for(u32bit i = 0; i != y_size; ++i)
+      {
+      const word y_i = y[i];
+
+      word carry = 0;
+
+      for(u32bit j = 0; j != x_size_8; j += 8)
+         carry = word8_madd3(z + i + j, x + j, y_i, carry);
+
+      for(u32bit j = x_size_8; j != x_size; ++j)
+         z[i+j] = word_madd3(x[j], y_i, z[i+j], &carry);
+
+      z[x_size+i] = carry;
+      }
+   }
+
+/*************************************************
+* Simple O(N^2) Squaring
+
+This is exactly the same algorithm as bigint_simple_mul,
+however because C/C++ compilers suck at alias analysis it
+is good to have the version where the compiler knows
+that x == y
+*************************************************/
+void bigint_simple_sqr(word z[], const word x[], u32bit x_size)
+   {
+   const u32bit x_size_8 = x_size - (x_size % 8);
+
+   clear_mem(z, 2*x_size);
+
+   for(u32bit i = 0; i != x_size; ++i)
+      {
+      const word x_i = x[i];
+      word carry = 0;
+
+      for(u32bit j = 0; j != x_size_8; j += 8)
+         carry = word8_madd3(z + i + j, x + j, x_i, carry);
+
+      for(u32bit j = x_size_8; j != x_size; ++j)
+         z[i+j] = word_madd3(x[j], x_i, z[i+j], &carry);
+
+      z[x_size+i] = carry;
+      }
+   }
+
+}
+
+}
diff --git a/src/bigint/mp_numth.cpp b/src/bigint/mp_numth.cpp
new file mode 100644
index 000000000..b024d2e2d
--- /dev/null
+++ b/src/bigint/mp_numth.cpp
@@ -0,0 +1,69 @@
+/*************************************************
+* Fused and Important MP Algorithms Source File  *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/numthry.h>
+#include <botan/mp_core.h>
+#include <botan/util.h>
+#include <algorithm>
+
+namespace Botan {
+
+/*************************************************
+* Square a BigInt                                *
+*************************************************/
+BigInt square(const BigInt& x)
+   {
+   const u32bit x_sw = x.sig_words();
+
+   BigInt z(BigInt::Positive, round_up(2*x_sw, 16));
+   SecureVector<word> workspace(z.size());
+
+   bigint_sqr(z.get_reg(), z.size(), workspace,
+              x.data(), x.size(), x_sw);
+   return z;
+   }
+
+/*************************************************
+* Multiply-Add Operation                         *
+*************************************************/
+BigInt mul_add(const BigInt& a, const BigInt& b, const BigInt& c)
+   {
+   if(c.is_negative() || c.is_zero())
+      throw Invalid_Argument("mul_add: Third argument must be > 0");
+
+   BigInt::Sign sign = BigInt::Positive;
+   if(a.sign() != b.sign())
+      sign = BigInt::Negative;
+
+   const u32bit a_sw = a.sig_words();
+   const u32bit b_sw = b.sig_words();
+   const u32bit c_sw = c.sig_words();
+
+   BigInt r(sign, std::max(a.size() + b.size(), c_sw) + 1);
+   SecureVector<word> workspace(r.size());
+
+   bigint_mul(r.get_reg(), r.size(), workspace,
+              a.data(), a.size(), a_sw,
+              b.data(), b.size(), b_sw);
+   const u32bit r_size = std::max(r.sig_words(), c_sw);
+   bigint_add2(r.get_reg(), r_size, c.data(), c_sw);
+   return r;
+   }
+
+/*************************************************
+* Subtract-Multiply Operation                    *
+*************************************************/
+BigInt sub_mul(const BigInt& a, const BigInt& b, const BigInt& c)
+   {
+   if(a.is_negative() || b.is_negative())
+      throw Invalid_Argument("sub_mul: First two arguments must be >= 0");
+
+   BigInt r = a;
+   r -= b;
+   r *= c;
+   return r;
+   }
+
+}
diff --git a/src/bigint/mp_shift.cpp b/src/bigint/mp_shift.cpp
new file mode 100644
index 000000000..033774e46
--- /dev/null
+++ b/src/bigint/mp_shift.cpp
@@ -0,0 +1,136 @@
+/*************************************************
+* MP Shift Algorithms Source File                *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/mp_core.h>
+#include <botan/mem_ops.h>
+
+namespace Botan {
+
+extern "C" {
+
+/*************************************************
+* Single Operand Left Shift                      *
+*************************************************/
+void bigint_shl1(word x[], u32bit x_size, u32bit word_shift, u32bit bit_shift)
+   {
+   if(word_shift)
+      {
+      for(u32bit j = 1; j != x_size + 1; ++j)
+         x[(x_size - j) + word_shift] = x[x_size - j];
+      clear_mem(x, word_shift);
+      }
+
+   if(bit_shift)
+      {
+      word carry = 0;
+      for(u32bit j = word_shift; j != x_size + word_shift + 1; ++j)
+         {
+         word temp = x[j];
+         x[j] = (temp << bit_shift) | carry;
+         carry = (temp >> (MP_WORD_BITS - bit_shift));
+         }
+      }
+   }
+
+/*************************************************
+* Single Operand Right Shift                     *
+*************************************************/
+void bigint_shr1(word x[], u32bit x_size, u32bit word_shift, u32bit bit_shift)
+   {
+   if(x_size < word_shift)
+      {
+      clear_mem(x, x_size);
+      return;
+      }
+
+   if(word_shift)
+      {
+      copy_mem(x, x + word_shift, x_size - word_shift);
+      clear_mem(x + x_size - word_shift, word_shift);
+      }
+
+   if(bit_shift)
+      {
+      word carry = 0;
+
+      u32bit top = x_size - word_shift;
+
+      while(top >= 4)
+         {
+         word w = x[top-1];
+         x[top-1] = (w >> bit_shift) | carry;
+         carry = (w << (MP_WORD_BITS - bit_shift));
+
+         w = x[top-2];
+         x[top-2] = (w >> bit_shift) | carry;
+         carry = (w << (MP_WORD_BITS - bit_shift));
+
+         w = x[top-3];
+         x[top-3] = (w >> bit_shift) | carry;
+         carry = (w << (MP_WORD_BITS - bit_shift));
+
+         w = x[top-4];
+         x[top-4] = (w >> bit_shift) | carry;
+         carry = (w << (MP_WORD_BITS - bit_shift));
+
+         top -= 4;
+         }
+
+      while(top)
+         {
+         word w = x[top-1];
+         x[top-1] = (w >> bit_shift) | carry;
+         carry = (w << (MP_WORD_BITS - bit_shift));
+
+         top--;
+         }
+      }
+   }
+
+/*************************************************
+* Two Operand Left Shift                         *
+*************************************************/
+void bigint_shl2(word y[], const word x[], u32bit x_size,
+                 u32bit word_shift, u32bit bit_shift)
+   {
+   for(u32bit j = 0; j != x_size; ++j)
+      y[j + word_shift] = x[j];
+   if(bit_shift)
+      {
+      word carry = 0;
+      for(u32bit j = word_shift; j != x_size + word_shift + 1; ++j)
+         {
+         word w = y[j];
+         y[j] = (w << bit_shift) | carry;
+         carry = (w >> (MP_WORD_BITS - bit_shift));
+         }
+      }
+   }
+
+/*************************************************
+* Two Operand Right Shift                        *
+*************************************************/
+void bigint_shr2(word y[], const word x[], u32bit x_size,
+                 u32bit word_shift, u32bit bit_shift)
+   {
+   if(x_size < word_shift) return;
+
+   for(u32bit j = 0; j != x_size - word_shift; ++j)
+      y[j] = x[j + word_shift];
+   if(bit_shift)
+      {
+      word carry = 0;
+      for(u32bit j = x_size - word_shift; j > 0; --j)
+         {
+         word w = y[j-1];
+         y[j-1] = (w >> bit_shift) | carry;
+         carry = (w << (MP_WORD_BITS - bit_shift));
+         }
+      }
+   }
+
+}
+
+}
diff --git a/src/bigint/mp_types.h b/src/bigint/mp_types.h
new file mode 100644
index 000000000..81b6d7395
--- /dev/null
+++ b/src/bigint/mp_types.h
@@ -0,0 +1,31 @@
+/*************************************************
+* Low Level MPI Types Header File                *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#ifndef BOTAN_MPI_TYPES_H__
+#define BOTAN_MPI_TYPES_H__
+
+#include <botan/types.h>
+
+namespace Botan {
+
+#if (BOTAN_MP_WORD_BITS == 8)
+  typedef byte word;
+#elif (BOTAN_MP_WORD_BITS == 16)
+  typedef u16bit word;
+#elif (BOTAN_MP_WORD_BITS == 32)
+  typedef u32bit word;
+#elif (BOTAN_MP_WORD_BITS == 64)
+  typedef u64bit word;
+#else
+  #error BOTAN_MP_WORD_BITS must be 8, 16, 32, or 64
+#endif
+
+const word MP_WORD_MASK = ~static_cast<word>(0);
+const word MP_WORD_TOP_BIT = static_cast<word>(1) << (8*sizeof(word) - 1);
+const word MP_WORD_MAX = MP_WORD_MASK;
+
+}
+
+#endif
diff --git a/src/bigint/numthry.cpp b/src/bigint/numthry.cpp
new file mode 100644
index 000000000..ffd523e82
--- /dev/null
+++ b/src/bigint/numthry.cpp
@@ -0,0 +1,339 @@
+/*************************************************
+* Number Theory Source File                      *
+* (C) 1999-2008 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/numthry.h>
+#include <botan/bit_ops.h>
+#include <algorithm>
+
+namespace Botan {
+
+namespace {
+
+/*************************************************
+* Miller-Rabin Iterations                        *
+*************************************************/
+u32bit miller_rabin_test_iterations(u32bit bits, bool verify)
+   {
+   struct mapping { u32bit bits; u32bit verify_iter; u32bit check_iter; };
+
+   static const mapping tests[] = {
+      {   50, 55, 25 },
+      {  100, 38, 22 },
+      {  160, 32, 18 },
+      {  163, 31, 17 },
+      {  168, 30, 16 },
+      {  177, 29, 16 },
+      {  181, 28, 15 },
+      {  185, 27, 15 },
+      {  190, 26, 15 },
+      {  195, 25, 14 },
+      {  201, 24, 14 },
+      {  208, 23, 14 },
+      {  215, 22, 13 },
+      {  222, 21, 13 },
+      {  231, 20, 13 },
+      {  241, 19, 12 },
+      {  252, 18, 12 },
+      {  264, 17, 12 },
+      {  278, 16, 11 },
+      {  294, 15, 10 },
+      {  313, 14,  9 },
+      {  334, 13,  8 },
+      {  360, 12,  8 },
+      {  392, 11,  7 },
+      {  430, 10,  7 },
+      {  479,  9,  6 },
+      {  542,  8,  6 },
+      {  626,  7,  5 },
+      {  746,  6,  4 },
+      {  926,  5,  3 },
+      { 1232,  4,  2 },
+      { 1853,  3,  2 },
+      {    0,  0,  0 }
+   };
+
+   for(u32bit j = 0; tests[j].bits; ++j)
+      {
+      if(bits <= tests[j].bits)
+         {
+         if(verify)
+            return tests[j].verify_iter;
+         else
+            return tests[j].check_iter;
+         }
+      }
+   return 2;
+   }
+
+}
+
+/*************************************************
+* Return the number of 0 bits at the end of n    *
+*************************************************/
+u32bit low_zero_bits(const BigInt& n)
+   {
+   if(n.is_negative() || n.is_zero()) return 0;
+
+   u32bit low_zero = 0;
+
+   if(n.is_positive() && n.is_nonzero())
+      {
+      for(u32bit i = 0; i != n.size(); ++i)
+         {
+         word x = n[i];
+
+         if(x)
+            {
+            low_zero += ctz(x);
+            break;
+            }
+         else
+            low_zero += BOTAN_MP_WORD_BITS;
+         }
+      }
+
+   return low_zero;
+   }
+
+/*************************************************
+* Calculate the GCD                              *
+*************************************************/
+BigInt gcd(const BigInt& a, const BigInt& b)
+   {
+   if(a.is_zero() || b.is_zero()) return 0;
+   if(a == 1 || b == 1)           return 1;
+
+   BigInt x = a, y = b;
+   x.set_sign(BigInt::Positive);
+   y.set_sign(BigInt::Positive);
+   u32bit shift = std::min(low_zero_bits(x), low_zero_bits(y));
+
+   x >>= shift;
+   y >>= shift;
+
+   while(x.is_nonzero())
+      {
+      x >>= low_zero_bits(x);
+      y >>= low_zero_bits(y);
+      if(x >= y) { x -= y; x >>= 1; }
+      else       { y -= x; y >>= 1; }
+      }
+
+   return (y << shift);
+   }
+
+/*************************************************
+* Calculate the LCM                              *
+*************************************************/
+BigInt lcm(const BigInt& a, const BigInt& b)
+   {
+   return ((a * b) / gcd(a, b));
+   }
+
+/*************************************************
+* Find the Modular Inverse                       *
+*************************************************/
+BigInt inverse_mod(const BigInt& n, const BigInt& mod)
+   {
+   if(mod.is_zero())
+      throw BigInt::DivideByZero();
+   if(mod.is_negative() || n.is_negative())
+      throw Invalid_Argument("inverse_mod: arguments must be non-negative");
+
+   if(n.is_zero() || (n.is_even() && mod.is_even()))
+      return 0;
+
+   BigInt x = mod, y = n, u = mod, v = n;
+   BigInt A = 1, B = 0, C = 0, D = 1;
+
+   while(u.is_nonzero())
+      {
+      u32bit zero_bits = low_zero_bits(u);
+      u >>= zero_bits;
+      for(u32bit j = 0; j != zero_bits; ++j)
+         {
+         if(A.is_odd() || B.is_odd())
+            { A += y; B -= x; }
+         A >>= 1; B >>= 1;
+         }
+
+      zero_bits = low_zero_bits(v);
+      v >>= zero_bits;
+      for(u32bit j = 0; j != zero_bits; ++j)
+         {
+         if(C.is_odd() || D.is_odd())
+            { C += y; D -= x; }
+         C >>= 1; D >>= 1;
+         }
+
+      if(u >= v) { u -= v; A -= C; B -= D; }
+      else       { v -= u; C -= A; D -= B; }
+      }
+
+   if(v != 1)
+      return 0;
+
+   while(D.is_negative()) D += mod;
+   while(D >= mod) D -= mod;
+
+   return D;
+   }
+
+/*************************************************
+* Modular Exponentiation                         *
+*************************************************/
+BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod)
+   {
+   Power_Mod pow_mod(mod);
+   pow_mod.set_base(base);
+   pow_mod.set_exponent(exp);
+   return pow_mod.execute();
+   }
+
+/*************************************************
+* Do simple tests of primality                   *
+*************************************************/
+s32bit simple_primality_tests(const BigInt& n)
+   {
+   const s32bit NOT_PRIME = -1, UNKNOWN = 0, PRIME = 1;
+
+   if(n == 2)
+      return PRIME;
+   if(n <= 1 || n.is_even())
+      return NOT_PRIME;
+
+   if(n <= PRIMES[PRIME_TABLE_SIZE-1])
+      {
+      const word num = n.word_at(0);
+      for(u32bit j = 0; PRIMES[j]; ++j)
+         {
+         if(num == PRIMES[j]) return PRIME;
+         if(num <  PRIMES[j]) return NOT_PRIME;
+         }
+      return NOT_PRIME;
+      }
+
+   u32bit check_first = std::min(n.bits() / 32, PRIME_PRODUCTS_TABLE_SIZE);
+   for(u32bit j = 0; j != check_first; ++j)
+      if(gcd(n, PRIME_PRODUCTS[j]) != 1)
+         return NOT_PRIME;
+
+   return UNKNOWN;
+   }
+
+/*************************************************
+* Fast check of primality                        *
+*************************************************/
+bool check_prime(const BigInt& n, RandomNumberGenerator& rng)
+   {
+   return run_primality_tests(rng, n, 0);
+   }
+
+/*************************************************
+* Test for primality                             *
+*************************************************/
+bool is_prime(const BigInt& n, RandomNumberGenerator& rng)
+   {
+   return run_primality_tests(rng, n, 1);
+   }
+
+/*************************************************
+* Verify primality                               *
+*************************************************/
+bool verify_prime(const BigInt& n, RandomNumberGenerator& rng)
+   {
+   return run_primality_tests(rng, n, 2);
+   }
+
+/*************************************************
+* Verify primality                               *
+*************************************************/
+bool run_primality_tests(RandomNumberGenerator& rng,
+                         const BigInt& n, u32bit level)
+   {
+   s32bit simple_tests = simple_primality_tests(n);
+   if(simple_tests) return (simple_tests == 1) ? true : false;
+   return passes_mr_tests(rng, n, level);
+   }
+
+/*************************************************
+* Test for primaility using Miller-Rabin         *
+*************************************************/
+bool passes_mr_tests(RandomNumberGenerator& rng,
+                     const BigInt& n, u32bit level)
+   {
+   const u32bit PREF_NONCE_BITS = 40;
+
+   if(level > 2)
+      level = 2;
+
+   MillerRabin_Test mr(n);
+
+   if(!mr.passes_test(2))
+      return false;
+
+   if(level == 0)
+      return true;
+
+   const u32bit NONCE_BITS = std::min(n.bits() - 1, PREF_NONCE_BITS);
+
+   const bool verify = (level == 2);
+
+   u32bit tests = miller_rabin_test_iterations(n.bits(), verify);
+
+   BigInt nonce;
+   for(u32bit j = 0; j != tests; ++j)
+      {
+      if(verify) nonce.randomize(rng, NONCE_BITS);
+      else       nonce = PRIMES[j];
+
+      if(!mr.passes_test(nonce))
+         return false;
+      }
+   return true;
+   }
+
+/*************************************************
+* Miller-Rabin Test                              *
+*************************************************/
+bool MillerRabin_Test::passes_test(const BigInt& a)
+   {
+   if(a < 2 || a >= n_minus_1)
+      throw Invalid_Argument("Bad size for nonce in Miller-Rabin test");
+
+   BigInt y = pow_mod(a);
+   if(y == 1 || y == n_minus_1)
+      return true;
+
+   for(u32bit j = 1; j != s; ++j)
+      {
+      y = reducer.square(y);
+
+      if(y == 1)
+         return false;
+      if(y == n_minus_1)
+         return true;
+      }
+   return false;
+   }
+
+/*************************************************
+* Miller-Rabin Constructor                       *
+*************************************************/
+MillerRabin_Test::MillerRabin_Test(const BigInt& num)
+   {
+   if(num.is_even() || num < 3)
+      throw Invalid_Argument("MillerRabin_Test: Invalid number for testing");
+
+   n = num;
+   n_minus_1 = n - 1;
+   s = low_zero_bits(n_minus_1);
+   r = n_minus_1 >> s;
+
+   pow_mod = Fixed_Exponent_Power_Mod(r, n);
+   reducer = Modular_Reducer(n);
+   }
+
+}
diff --git a/src/bigint/numthry.h b/src/bigint/numthry.h
new file mode 100644
index 000000000..371621c2d
--- /dev/null
+++ b/src/bigint/numthry.h
@@ -0,0 +1,103 @@
+/*************************************************
+* Number Theory Header File                      *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#ifndef BOTAN_NUMBTHRY_H__
+#define BOTAN_NUMBTHRY_H__
+
+#include <botan/base.h>
+#include <botan/bigint.h>
+#include <botan/reducer.h>
+#include <botan/pow_mod.h>
+
+namespace Botan {
+
+/*************************************************
+* Fused Arithmetic Operations                    *
+*************************************************/
+BigInt BOTAN_DLL mul_add(const BigInt&, const BigInt&, const BigInt&);
+BigInt BOTAN_DLL sub_mul(const BigInt&, const BigInt&, const BigInt&);
+
+/*************************************************
+* Number Theory Functions                        *
+*************************************************/
+inline BigInt abs(const BigInt& n) { return n.abs(); }
+
+void BOTAN_DLL divide(const BigInt&, const BigInt&, BigInt&, BigInt&);
+
+BigInt BOTAN_DLL gcd(const BigInt&, const BigInt&);
+BigInt BOTAN_DLL lcm(const BigInt&, const BigInt&);
+
+BigInt BOTAN_DLL square(const BigInt&);
+BigInt BOTAN_DLL inverse_mod(const BigInt&, const BigInt&);
+s32bit BOTAN_DLL jacobi(const BigInt&, const BigInt&);
+
+BigInt BOTAN_DLL power_mod(const BigInt&, const BigInt&, const BigInt&);
+
+/*************************************************
+* Compute the square root of x modulo a prime    *
+* using the Shanks-Tonnelli algorithm            *
+*************************************************/
+BigInt ressol(const BigInt& x, const BigInt& p);
+
+/*************************************************
+* Utility Functions                              *
+*************************************************/
+u32bit BOTAN_DLL low_zero_bits(const BigInt&);
+
+/*************************************************
+* Primality Testing                              *
+*************************************************/
+bool BOTAN_DLL check_prime(const BigInt&, RandomNumberGenerator&);
+bool BOTAN_DLL is_prime(const BigInt&, RandomNumberGenerator&);
+bool BOTAN_DLL verify_prime(const BigInt&, RandomNumberGenerator&);
+
+s32bit BOTAN_DLL simple_primality_tests(const BigInt&);
+
+bool BOTAN_DLL passes_mr_tests(RandomNumberGenerator&,
+                               const BigInt&, u32bit = 1);
+
+bool BOTAN_DLL run_primality_tests(RandomNumberGenerator&,
+                                   const BigInt&, u32bit = 1);
+
+/*************************************************
+* Random Number Generation                       *
+*************************************************/
+BigInt BOTAN_DLL random_integer(RandomNumberGenerator&,
+                                const BigInt&, const BigInt&);
+
+BigInt BOTAN_DLL random_prime(RandomNumberGenerator&,
+                              u32bit bits, const BigInt& coprime = 1,
+                              u32bit equiv = 1, u32bit equiv_mod = 2);
+
+BigInt BOTAN_DLL random_safe_prime(RandomNumberGenerator&,
+                                   u32bit);
+
+/*************************************************
+* Prime Numbers                                  *
+*************************************************/
+const u32bit PRIME_TABLE_SIZE = 6541;
+const u32bit PRIME_PRODUCTS_TABLE_SIZE = 256;
+
+extern const u16bit BOTAN_DLL PRIMES[];
+extern const u64bit PRIME_PRODUCTS[];
+
+/*************************************************
+* Miller-Rabin Primality Tester                  *
+*************************************************/
+class BOTAN_DLL MillerRabin_Test
+   {
+   public:
+      bool passes_test(const BigInt&);
+      MillerRabin_Test(const BigInt&);
+   private:
+      BigInt n, r, n_minus_1;
+      u32bit s;
+      Fixed_Exponent_Power_Mod pow_mod;
+      Modular_Reducer reducer;
+   };
+
+}
+
+#endif
diff --git a/src/bigint/pow_mod.cpp b/src/bigint/pow_mod.cpp
new file mode 100644
index 000000000..17ca7b796
--- /dev/null
+++ b/src/bigint/pow_mod.cpp
@@ -0,0 +1,155 @@
+/*************************************************
+* Modular Exponentiation Proxy Source File       *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/pow_mod.h>
+#include <botan/engine.h>
+
+namespace Botan {
+
+/*************************************************
+* Power_Mod Constructor                          *
+*************************************************/
+Power_Mod::Power_Mod(const BigInt& n, Usage_Hints hints)
+   {
+   core = 0;
+   set_modulus(n, hints);
+   }
+
+/*************************************************
+* Power_Mod Copy Constructor                     *
+*************************************************/
+Power_Mod::Power_Mod(const Power_Mod& other)
+   {
+   core = 0;
+   if(other.core)
+      core = other.core->copy();
+   }
+
+/*************************************************
+* Power_Mod Assignment Operator                  *
+*************************************************/
+Power_Mod& Power_Mod::operator=(const Power_Mod& other)
+   {
+   delete core;
+   core = 0;
+   if(other.core)
+      core = other.core->copy();
+   return (*this);
+   }
+
+/*************************************************
+* Power_Mod Destructor                           *
+*************************************************/
+Power_Mod::~Power_Mod()
+   {
+   delete core;
+   }
+
+/*************************************************
+* Set the modulus                                *
+*************************************************/
+void Power_Mod::set_modulus(const BigInt& n, Usage_Hints hints) const
+   {
+   delete core;
+   core = ((n == 0) ? 0 : Engine_Core::mod_exp(n, hints));
+   }
+
+/*************************************************
+* Set the base                                   *
+*************************************************/
+void Power_Mod::set_base(const BigInt& b) const
+   {
+   if(b.is_zero() || b.is_negative())
+      throw Invalid_Argument("Power_Mod::set_base: arg must be > 0");
+
+   if(!core)
+      throw Internal_Error("Power_Mod::set_base: core was NULL");
+   core->set_base(b);
+   }
+
+/*************************************************
+* Set the exponent                               *
+*************************************************/
+void Power_Mod::set_exponent(const BigInt& e) const
+   {
+   if(e.is_negative())
+      throw Invalid_Argument("Power_Mod::set_exponent: arg must be > 0");
+
+   if(!core)
+      throw Internal_Error("Power_Mod::set_exponent: core was NULL");
+   core->set_exponent(e);
+   }
+
+/*************************************************
+* Compute the result                             *
+*************************************************/
+BigInt Power_Mod::execute() const
+   {
+   if(!core)
+      throw Internal_Error("Power_Mod::execute: core was NULL");
+   return core->execute();
+   }
+
+namespace {
+
+/*************************************************
+* Choose potentially useful hints                *
+*************************************************/
+Power_Mod::Usage_Hints choose_base_hints(const BigInt& b, const BigInt& n)
+   {
+   if(b == 2)
+      return Power_Mod::Usage_Hints(Power_Mod::BASE_IS_2 |
+                                    Power_Mod::BASE_IS_SMALL);
+
+   const u32bit b_bits = b.bits();
+   const u32bit n_bits = n.bits();
+
+   if(b_bits < n_bits / 32)
+      return Power_Mod::BASE_IS_SMALL;
+   if(b_bits > n_bits / 4)
+      return Power_Mod::BASE_IS_LARGE;
+
+   return Power_Mod::NO_HINTS;
+   }
+
+/*************************************************
+* Choose potentially useful hints                *
+*************************************************/
+Power_Mod::Usage_Hints choose_exp_hints(const BigInt& e, const BigInt& n)
+   {
+   const u32bit e_bits = e.bits();
+   const u32bit n_bits = n.bits();
+
+   if(e_bits < n_bits / 32)
+      return Power_Mod::BASE_IS_SMALL;
+   if(e_bits > n_bits / 4)
+      return Power_Mod::BASE_IS_LARGE;
+   return Power_Mod::NO_HINTS;
+   }
+
+}
+
+/*************************************************
+* Fixed_Exponent_Power_Mod Constructor           *
+*************************************************/
+Fixed_Exponent_Power_Mod::Fixed_Exponent_Power_Mod(const BigInt& e,
+                                                   const BigInt& n,
+                                                   Usage_Hints hints) :
+   Power_Mod(n, Usage_Hints(hints | EXP_IS_FIXED | choose_exp_hints(e, n)))
+   {
+   set_exponent(e);
+   }
+
+/*************************************************
+* Fixed_Base_Power_Mod Constructor               *
+*************************************************/
+Fixed_Base_Power_Mod::Fixed_Base_Power_Mod(const BigInt& b, const BigInt& n,
+                                           Usage_Hints hints) :
+   Power_Mod(n, Usage_Hints(hints | BASE_IS_FIXED | choose_base_hints(b, n)))
+   {
+   set_base(b);
+   }
+
+}
diff --git a/src/bigint/pow_mod.h b/src/bigint/pow_mod.h
new file mode 100644
index 000000000..37e0871da
--- /dev/null
+++ b/src/bigint/pow_mod.h
@@ -0,0 +1,91 @@
+/*************************************************
+* Modular Exponentiator Header File              *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#ifndef BOTAN_POWER_MOD_H__
+#define BOTAN_POWER_MOD_H__
+
+#include <botan/bigint.h>
+
+namespace Botan {
+
+/*************************************************
+* Modular Exponentiator Interface                *
+*************************************************/
+class BOTAN_DLL Modular_Exponentiator
+   {
+   public:
+      virtual void set_base(const BigInt&) = 0;
+      virtual void set_exponent(const BigInt&) = 0;
+      virtual BigInt execute() const = 0;
+      virtual Modular_Exponentiator* copy() const = 0;
+      virtual ~Modular_Exponentiator() {}
+   };
+
+/*************************************************
+* Modular Exponentiator Proxy                    *
+*************************************************/
+class BOTAN_DLL Power_Mod
+   {
+   public:
+      enum Usage_Hints {
+         NO_HINTS        = 0x0000,
+
+         BASE_IS_FIXED   = 0x0001,
+         BASE_IS_SMALL   = 0x0002,
+         BASE_IS_LARGE   = 0x0004,
+         BASE_IS_2       = 0x0008,
+
+         EXP_IS_FIXED    = 0x0100,
+         EXP_IS_SMALL    = 0x0200,
+         EXP_IS_LARGE    = 0x0400
+      };
+
+      void set_modulus(const BigInt&, Usage_Hints = NO_HINTS) const;
+      void set_base(const BigInt&) const;
+      void set_exponent(const BigInt&) const;
+
+      BigInt execute() const;
+
+      Power_Mod& operator=(const Power_Mod&);
+
+      Power_Mod(const BigInt& = 0, Usage_Hints = NO_HINTS);
+      Power_Mod(const Power_Mod&);
+      ~Power_Mod();
+   private:
+      mutable Modular_Exponentiator* core;
+      Usage_Hints hints;
+   };
+
+/*************************************************
+* Fixed Exponent Modular Exponentiator Proxy     *
+*************************************************/
+class BOTAN_DLL Fixed_Exponent_Power_Mod : public Power_Mod
+   {
+   public:
+      BigInt operator()(const BigInt& b) const
+         { set_base(b); return execute(); }
+
+      Fixed_Exponent_Power_Mod() {}
+      Fixed_Exponent_Power_Mod(const BigInt&, const BigInt&,
+                               Usage_Hints = NO_HINTS);
+   };
+
+/*************************************************
+* Fixed Base Modular Exponentiator Proxy         *
+*************************************************/
+class BOTAN_DLL Fixed_Base_Power_Mod : public Power_Mod
+   {
+   public:
+      BigInt operator()(const BigInt& e) const
+         { set_exponent(e); return execute(); }
+
+      Fixed_Base_Power_Mod() {}
+      Fixed_Base_Power_Mod(const BigInt&, const BigInt&,
+                           Usage_Hints = NO_HINTS);
+   };
+
+}
+
+#endif
diff --git a/src/bigint/powm_fw.cpp b/src/bigint/powm_fw.cpp
new file mode 100644
index 000000000..c29b9f311
--- /dev/null
+++ b/src/bigint/powm_fw.cpp
@@ -0,0 +1,102 @@
+/*************************************************
+* Fixed Window Exponentiation Source File        *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/def_powm.h>
+#include <botan/numthry.h>
+#include <vector>
+
+namespace Botan {
+
+namespace {
+
+/*************************************************
+* Try to choose a good window size               *
+*************************************************/
+u32bit choose_window_bits(u32bit exp_bits, u32bit,
+                          Power_Mod::Usage_Hints hints)
+   {
+   static const u32bit wsize[][2] = {
+      { 2048, 7 }, { 1024, 6 }, { 256, 5 }, { 128, 4 }, { 64, 3 }, { 0, 0 }
+   };
+
+   u32bit window_bits = 3;
+
+   if(exp_bits)
+      {
+      for(u32bit j = 0; wsize[j][0]; ++j)
+         {
+         if(exp_bits >= wsize[j][0])
+            {
+            window_bits += wsize[j][1];
+            break;
+            }
+         }
+      }
+
+   if(hints & Power_Mod::EXP_IS_FIXED)
+      window_bits += 2;
+   if(hints & Power_Mod::EXP_IS_LARGE)
+      window_bits += 2;
+   if(hints & Power_Mod::BASE_IS_FIXED)
+      ++window_bits;
+
+   return window_bits;
+   }
+
+}
+
+/*************************************************
+* Set the exponent                               *
+*************************************************/
+void Fixed_Window_Exponentiator::set_exponent(const BigInt& e)
+   {
+   exp = e;
+   }
+
+/*************************************************
+* Set the base                                   *
+*************************************************/
+void Fixed_Window_Exponentiator::set_base(const BigInt& base)
+   {
+   window_bits = choose_window_bits(exp.bits(), base.bits(), hints);
+
+   g.resize((1 << window_bits) - 1);
+   g[0] = base;
+   for(u32bit j = 1; j != g.size(); ++j)
+      g[j] = reducer.multiply(g[j-1], g[0]);
+   }
+
+/*************************************************
+* Compute the result                             *
+*************************************************/
+BigInt Fixed_Window_Exponentiator::execute() const
+   {
+   const u32bit exp_nibbles = (exp.bits() + window_bits - 1) / window_bits;
+
+   BigInt x = 1;
+   for(u32bit j = exp_nibbles; j > 0; --j)
+      {
+      for(u32bit k = 0; k != window_bits; ++k)
+         x = reducer.square(x);
+
+      u32bit nibble = exp.get_substring(window_bits*(j-1), window_bits);
+      if(nibble)
+         x = reducer.multiply(x, g[nibble-1]);
+      }
+   return x;
+   }
+
+/*************************************************
+* Fixed_Window_Exponentiator Constructor         *
+*************************************************/
+Fixed_Window_Exponentiator::Fixed_Window_Exponentiator(const BigInt& n,
+   Power_Mod::Usage_Hints hints)
+   {
+   reducer = Modular_Reducer(n);
+   this->hints = hints;
+   window_bits = 0;
+   }
+
+}
diff --git a/src/bigint/powm_mnt.cpp b/src/bigint/powm_mnt.cpp
new file mode 100644
index 000000000..6091d467a
--- /dev/null
+++ b/src/bigint/powm_mnt.cpp
@@ -0,0 +1,178 @@
+/*************************************************
+* Montgomery Exponentiation Source File          *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/def_powm.h>
+#include <botan/numthry.h>
+#include <botan/mp_core.h>
+
+namespace Botan {
+
+namespace {
+
+/*************************************************
+* Try to choose a good window size               *
+*************************************************/
+u32bit choose_window_bits(u32bit exp_bits, u32bit,
+                          Power_Mod::Usage_Hints hints)
+   {
+   static const u32bit wsize[][2] = {
+      { 2048, 4 }, { 1024, 3 }, { 256, 2 }, { 128, 1 }, { 0, 0 }
+   };
+
+   u32bit window_bits = 1;
+
+   if(exp_bits)
+      {
+      for(u32bit j = 0; wsize[j][0]; ++j)
+         {
+         if(exp_bits >= wsize[j][0])
+            {
+            window_bits += wsize[j][1];
+            break;
+            }
+         }
+      }
+
+   if(hints & Power_Mod::BASE_IS_FIXED)
+      window_bits += 2;
+   if(hints & Power_Mod::EXP_IS_LARGE)
+      ++window_bits;
+
+   return window_bits;
+   }
+
+/*************************************************
+* Montgomery Reduction                           *
+*************************************************/
+inline void montgomery_reduce(BigInt& out, MemoryRegion<word>& z_buf,
+                              const BigInt& x_bn, u32bit x_size, word u)
+   {
+   const word* x = x_bn.data();
+   word* z = z_buf.begin();
+   u32bit z_size = z_buf.size();
+
+   bigint_monty_redc(z, z_size, x, x_size, u);
+
+   out.get_reg().set(z + x_size, x_size + 1);
+   }
+
+}
+
+/*************************************************
+* Set the exponent                               *
+*************************************************/
+void Montgomery_Exponentiator::set_exponent(const BigInt& exp)
+   {
+   this->exp = exp;
+   exp_bits = exp.bits();
+   }
+
+/*************************************************
+* Set the base                                   *
+*************************************************/
+void Montgomery_Exponentiator::set_base(const BigInt& base)
+   {
+   window_bits = choose_window_bits(exp.bits(), base.bits(), hints);
+
+   g.resize((1 << window_bits) - 1);
+
+   SecureVector<word> z(2 * (mod_words + 1));
+   SecureVector<word> workspace(z.size());
+
+   g[0] = (base >= modulus) ? (base % modulus) : base;
+   bigint_mul(z.begin(), z.size(), workspace,
+              g[0].data(), g[0].size(), g[0].sig_words(),
+              R2.data(), R2.size(), R2.sig_words());
+
+   montgomery_reduce(g[0], z, modulus, mod_words, mod_prime);
+
+   const BigInt& x = g[0];
+   const u32bit x_sig = x.sig_words();
+
+   for(u32bit j = 1; j != g.size(); ++j)
+      {
+      const BigInt& y = g[j-1];
+      const u32bit y_sig = y.sig_words();
+
+      z.clear();
+      bigint_mul(z.begin(), z.size(), workspace,
+                 x.data(), x.size(), x_sig,
+                 y.data(), y.size(), y_sig);
+
+      montgomery_reduce(g[j], z, modulus, mod_words, mod_prime);
+      }
+   }
+
+/*************************************************
+* Compute the result                             *
+*************************************************/
+BigInt Montgomery_Exponentiator::execute() const
+   {
+   const u32bit exp_nibbles = (exp_bits + window_bits - 1) / window_bits;
+
+   BigInt x = R_mod;
+   SecureVector<word> z(2 * (mod_words + 1));
+   SecureVector<word> workspace(2 * (mod_words + 1));
+
+   for(u32bit j = exp_nibbles; j > 0; --j)
+      {
+      for(u32bit k = 0; k != window_bits; ++k)
+         {
+         z.clear();
+         bigint_sqr(z.begin(), z.size(), workspace,
+                    x.data(), x.size(), x.sig_words());
+
+         montgomery_reduce(x, z, modulus, mod_words, mod_prime);
+         }
+
+      u32bit nibble = exp.get_substring(window_bits*(j-1), window_bits);
+      if(nibble)
+         {
+         const BigInt& y = g[nibble-1];
+
+         z.clear();
+         bigint_mul(z.begin(), z.size(), workspace,
+                    x.data(), x.size(), x.sig_words(),
+                    y.data(), y.size(), y.sig_words());
+
+         montgomery_reduce(x, z, modulus, mod_words, mod_prime);
+         }
+      }
+
+   z.clear();
+   z.copy(x.data(), x.size());
+
+   montgomery_reduce(x, z, modulus, mod_words, mod_prime);
+   return x;
+   }
+
+/*************************************************
+* Montgomery_Exponentiator Constructor           *
+*************************************************/
+Montgomery_Exponentiator::Montgomery_Exponentiator(const BigInt& mod,
+   Power_Mod::Usage_Hints hints)
+   {
+   if(!mod.is_positive())
+      throw Exception("Montgomery_Exponentiator: modulus must be positive");
+   if(mod.is_even())
+      throw Exception("Montgomery_Exponentiator: modulus must be odd");
+
+   window_bits = 0;
+   this->hints = hints;
+   modulus = mod;
+
+   mod_words = modulus.sig_words();
+
+   BigInt mod_prime_bn(BigInt::Power2, MP_WORD_BITS);
+   mod_prime = (mod_prime_bn - inverse_mod(modulus, mod_prime_bn)).word_at(0);
+
+   R_mod = BigInt(BigInt::Power2, MP_WORD_BITS * mod_words);
+   R_mod %= modulus;
+
+   R2 = BigInt(BigInt::Power2, 2 * MP_WORD_BITS * mod_words);
+   R2 %= modulus;
+   }
+
+}
diff --git a/src/bigint/reducer.cpp b/src/bigint/reducer.cpp
new file mode 100644
index 000000000..47c5c20fc
--- /dev/null
+++ b/src/bigint/reducer.cpp
@@ -0,0 +1,95 @@
+/*************************************************
+* Modular Reducer Source File                    *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#include <botan/reducer.h>
+#include <botan/numthry.h>
+#include <botan/mp_core.h>
+
+namespace Botan {
+
+/*************************************************
+* Modular_Reducer Constructor                    *
+*************************************************/
+Modular_Reducer::Modular_Reducer(const BigInt& mod)
+   {
+   if(mod <= 0)
+      throw Invalid_Argument("Modular_Reducer: modulus must be positive");
+
+   modulus = mod;
+   mod_words = modulus.sig_words();
+
+   modulus_2 = Botan::square(modulus);
+   mod2_words = modulus_2.sig_words();
+
+   mu = BigInt(BigInt::Power2, 2 * MP_WORD_BITS * mod_words) / modulus;
+   mu_words = mu.sig_words();
+   }
+
+/*************************************************
+* Barrett Reduction                              *
+*************************************************/
+BigInt Modular_Reducer::reduce(const BigInt& x) const
+   {
+   if(mod_words == 0)
+      throw Invalid_State("Modular_Reducer: Never initalized");
+
+   BigInt t1 = x;
+   t1.set_sign(BigInt::Positive);
+
+   if(t1 < modulus)
+      {
+      if(x.is_negative() && t1.is_nonzero())
+         return modulus - t1;
+      return x;
+      }
+
+   if(t1 >= modulus_2)
+      return (x % modulus);
+
+   t1 >>= (MP_WORD_BITS * (mod_words - 1));
+   t1 *= mu;
+   t1 >>= (MP_WORD_BITS * (mod_words + 1));
+
+   t1 *= modulus;
+   t1.mask_bits(MP_WORD_BITS * (mod_words+1));
+
+   BigInt t2 = x;
+   t2.set_sign(BigInt::Positive);
+   t2.mask_bits(MP_WORD_BITS * (mod_words+1));
+
+   t1 = t2 - t1;
+
+   if(t1.is_negative())
+      {
+      BigInt b_to_k1(BigInt::Power2, MP_WORD_BITS * (mod_words+1));
+      t1 += b_to_k1;
+      }
+
+   while(t1 >= modulus)
+      t1 -= modulus;
+
+   if(x.is_negative() && t1.is_nonzero())
+      t1 = modulus - t1;
+
+   return t1;
+   }
+
+/*************************************************
+* Multiply, followed by a reduction              *
+*************************************************/
+BigInt Modular_Reducer::multiply(const BigInt& x, const BigInt& y) const
+   {
+   return reduce(x * y);
+   }
+
+/*************************************************
+* Square, followed by a reduction                *
+*************************************************/
+BigInt Modular_Reducer::square(const BigInt& x) const
+   {
+   return reduce(Botan::square(x));
+   }
+
+}
diff --git a/src/bigint/reducer.h b/src/bigint/reducer.h
new file mode 100644
index 000000000..48008e73b
--- /dev/null
+++ b/src/bigint/reducer.h
@@ -0,0 +1,34 @@
+/*************************************************
+* Modular Reducer Header File                    *
+* (C) 1999-2007 Jack Lloyd                       *
+*************************************************/
+
+#ifndef BOTAN_MODARITH_H__
+#define BOTAN_MODARITH_H__
+
+#include <botan/bigint.h>
+
+namespace Botan {
+
+/*************************************************
+* Modular Reducer                                *
+*************************************************/
+class BOTAN_DLL Modular_Reducer
+   {
+   public:
+      BigInt multiply(const BigInt&, const BigInt&) const;
+      BigInt square(const BigInt&) const;
+      BigInt reduce(const BigInt&) const;
+
+      bool initialized() const { return (mod_words != 0); }
+
+      Modular_Reducer() { mod_words = 0; }
+      Modular_Reducer(const BigInt&);
+   private:
+      BigInt modulus, modulus_2, mu;
+      u32bit mod_words, mod2_words, mu_words;
+   };
+
+}
+
+#endif
diff --git a/src/bigint/ressol.cpp b/src/bigint/ressol.cpp
new file mode 100644
index 000000000..0cd2b988a
--- /dev/null
+++ b/src/bigint/ressol.cpp
@@ -0,0 +1,82 @@
+/*************************************************
+* Shanks-Tonnelli (RESSOL) Source File           *
+* (C) 2007-2008 Falko Strenzke, FlexSecure GmbH  *
+* (C) 2008 Jack Lloyd                            *
+*************************************************/
+
+#include <botan/numthry.h>
+#include <botan/reducer.h>
+
+#include <iostream>
+
+namespace Botan {
+
+/*************************************************
+* Shanks-Tonnelli algorithm                      *
+*************************************************/
+BigInt ressol(const BigInt& a, const BigInt& p)
+   {
+   if(a < 0)
+      throw Invalid_Argument("ressol(): a to solve for must be positive");
+   if(p <= 1)
+      throw Invalid_Argument("ressol(): prime must be > 1");
+
+   if(a == 0)
+      return 0;
+   if(p == 2)
+      return a;
+
+   if(jacobi(a, p) != 1) // not a quadratic residue
+      return -BigInt(1);
+
+   if(p % 4 == 3)
+      return power_mod(a, ((p+1) >> 2), p);
+
+   u32bit s = low_zero_bits(p - 1);
+   BigInt q = p >> s;
+
+   q -= 1;
+   q >>= 1;
+
+   Modular_Reducer mod_p(p);
+
+   BigInt r = power_mod(a, q, p);
+   BigInt n = mod_p.multiply(a, mod_p.square(r));
+   r = mod_p.multiply(r, a);
+
+   if(n == 1)
+      return r;
+
+   // find random non quadratic residue z
+   BigInt z = 2;
+   while(jacobi(z, p) == 1) // while z quadratic residue
+      ++z;
+
+   BigInt c = power_mod(z, (q << 1) + 1, p);
+
+   while(n > 1)
+      {
+      q = n;
+
+      u32bit i = 0;
+      while(q != 1)
+         {
+         q = mod_p.square(q);
+         ++i;
+         }
+      u32bit t = s;
+
+      if(t <= i)
+         return -BigInt(1);
+
+      c = power_mod(c, BigInt(BigInt::Power2, t-i-1), p);
+      r = mod_p.multiply(r, c);
+      c = mod_p.square(c);
+      n = mod_p.multiply(n, c);
+      s = i;
+      }
+
+   return r;
+   }
+
+}