Move GF(p) math code from pk/ecdsa to math/gfpmath

7 files changed, 2892 insertions, 0 deletions

diff --git a/src/math/gfpmath/curve_gfp.cpp b/src/math/gfpmath/curve_gfp.cpp
new file mode 100644
index 000000000..89fa74c49
--- /dev/null
+++ b/src/math/gfpmath/curve_gfp.cpp
@@ -0,0 +1,157 @@
+/******************************************************
+ * Elliptic curves over GF(p) (source file)           *
+ *                                                    *
+ * (C) 2007 Martin Doering                            *
+ *          [email protected]    *
+ *          Christoph Ludwig                          *
+ *          [email protected]                        *
+ *          Falko Strenzke                            *
+ *          [email protected]                    *
+ ******************************************************/
+
+#include <botan/curve_gfp.h>
+
+namespace Botan {
+
+CurveGFp::CurveGFp(GFpElement const& a, GFpElement const& b,
+                   const BigInt& p) :
+   mA(a), mB(b)
+   {
+   if(p != mA.get_p() || p != mB.get_p())
+      throw Invalid_Argument("could not construct curve: moduli of arguments differ");
+
+   std::tr1::shared_ptr<GFpModulus> p_mod = std::tr1::shared_ptr<GFpModulus>(new GFpModulus(p));
+   // the above is the creation of the GFpModuls object which will be shared point-wide
+   // (in the context of a point of course)
+   set_shrd_mod(p_mod);
+   }
+
+// copy constructor
+CurveGFp::CurveGFp(CurveGFp const& other)
+   :	mA(other.get_a()),
+        mB(other.get_b())
+   {
+   mp_mod = std::tr1::shared_ptr<GFpModulus>(new GFpModulus(*other.mp_mod));
+   //assert(mp_mod->p_equal_to(mA.get_p()));
+   //assert(mp_mod->p_equal_to(mB.get_p()));
+   set_shrd_mod(mp_mod);
+   if(other.mp_mres_a.get())
+      {
+      mp_mres_a = std::tr1::shared_ptr<GFpElement>(new GFpElement(*other.mp_mres_a));
+      }
+   if(other.mp_mres_b.get())
+      {
+      mp_mres_b = std::tr1::shared_ptr<GFpElement>(new GFpElement(*other.mp_mres_b));
+      }
+   if(other.mp_mres_one.get())
+      {
+      mp_mres_one = std::tr1::shared_ptr<GFpElement>(new GFpElement(*other.mp_mres_one));
+      }
+
+   }
+
+// assignment operator
+CurveGFp const& CurveGFp::operator=(CurveGFp const& other)
+   {
+   // for exception safety...
+   GFpElement a_tmp = other.mA;
+   GFpElement b_tmp = other.mB;
+   mA.swap(a_tmp);
+   mB.swap(b_tmp);
+
+   std::tr1::shared_ptr<GFpModulus> p_mod = std::tr1::shared_ptr<GFpModulus>(new GFpModulus(*other.mp_mod));
+   set_shrd_mod(p_mod);
+
+   // exception safety note: no problem if we have a throw from here on...
+   if(other.mp_mres_a.get())
+      {
+      mp_mres_a = std::tr1::shared_ptr<GFpElement>(new GFpElement(*other.mp_mres_a));
+      }
+   if(other.mp_mres_b.get())
+      {
+      mp_mres_b = std::tr1::shared_ptr<GFpElement>(new GFpElement(*other.mp_mres_b));
+      }
+   if(other.mp_mres_one.get())
+      {
+      mp_mres_one = std::tr1::shared_ptr<GFpElement>(new GFpElement(*other.mp_mres_one));
+      }
+   return *this;
+   }
+
+// getters
+GFpElement const CurveGFp::get_a() const
+   {
+   return mA;
+   }
+GFpElement const CurveGFp::get_b() const
+   {
+   return mB;
+   }
+
+BigInt const CurveGFp::get_p() const
+   {
+   //assert(mp_mod.get() != 0);
+   return mp_mod->get_p();
+   }
+
+// swaps the states of *this and other, does not throw
+void CurveGFp::swap(CurveGFp& other)
+   {
+   mA.swap(other.mA);
+   mB.swap(other.mB);
+   mp_mod.swap(other.mp_mod);
+   std::swap(mp_mres_a, other.mp_mres_a);
+   std::swap(mp_mres_b, other.mp_mres_b);
+   std::swap(mp_mres_one, other.mp_mres_one);
+   }
+
+void CurveGFp::set_shrd_mod(std::tr1::shared_ptr<GFpModulus> const mod)
+   {
+   mp_mod = mod;
+   mA.turn_off_sp_red_mul();// m.m. is not needed, must be trf. back
+   mB.turn_off_sp_red_mul();// m.m. is not needed, must be trf. back
+   //ok, above we destroy any evantually computated montg. mult. values,
+   // but that won't influence performance in usual applications
+   mA.set_shrd_mod(mod);
+   mB.set_shrd_mod(mod);
+   }
+
+GFpElement const CurveGFp::get_mres_a() const
+   {
+   if(mp_mres_a.get() == 0)
+      {
+      mp_mres_a = std::tr1::shared_ptr<GFpElement>(new GFpElement(mA));
+      mp_mres_a->turn_on_sp_red_mul();
+      mp_mres_a->get_mres();
+      }
+   return GFpElement(*mp_mres_a);
+   }
+
+GFpElement const CurveGFp::get_mres_b() const
+   {
+   if(mp_mres_b.get() == 0)
+      {
+      mp_mres_b = std::tr1::shared_ptr<GFpElement>(new GFpElement(mB));
+      mp_mres_b->turn_on_sp_red_mul();
+      mp_mres_b->get_mres();
+      }
+   return GFpElement(*mp_mres_b);
+   }
+
+std::tr1::shared_ptr<GFpElement const> const CurveGFp::get_mres_one() const
+   {
+   if(mp_mres_one.get() == 0)
+      {
+      mp_mres_one = std::tr1::shared_ptr<GFpElement>(new GFpElement(mp_mod->get_p(), 1));
+      mp_mres_one->turn_on_sp_red_mul();
+      mp_mres_one->get_mres();
+      }
+   return mp_mres_one;
+   }
+
+bool operator==(CurveGFp const& lhs, CurveGFp const& rhs)
+   {
+   return (lhs.get_p() == rhs.get_p() && lhs.get_a() == rhs.get_a() && lhs.get_b() == rhs.get_b());
+   }
+
+}
diff --git a/src/math/gfpmath/curve_gfp.h b/src/math/gfpmath/curve_gfp.h
new file mode 100644
index 000000000..49688b5dc
--- /dev/null
+++ b/src/math/gfpmath/curve_gfp.h
@@ -0,0 +1,165 @@
+/******************************************************
+ * Elliptic curves over GF(p) (header file)           *
+ *                                                    *
+ * (C) 2007 Martin Doering                             *
+ *          [email protected]    *
+ *          Christoph Ludwig                          *
+ *          [email protected]                        *
+ *          Falko Strenzke                            *
+ *          [email protected]                    *
+ ******************************************************/
+
+#ifndef BOTAN_EC_CURVE_GFP_H__
+#define BOTAN_EC_CURVE_GFP_H__
+
+#include <botan/gfp_element.h>
+#include <botan/bigint.h>
+#include <tr1/memory>
+
+namespace Botan {
+
+/**
+* This class represents an elliptic curve over GF(p)
+*/
+class CurveGFp
+   {
+   public:
+      /**
+      * Construct the elliptic curve E: y^2 = x^3 + ax + b over GF(p)
+      * @param a first coefficient
+      * @param b second coefficient
+      * @param p prime number of the field
+      */
+      CurveGFp(GFpElement const& a, GFpElement const& b,
+               const BigInt& p);
+
+      /**
+      * Copy constructor
+      * @param other The curve to clone
+      */
+      CurveGFp(CurveGFp const& other);
+
+      /**
+      * Assignment operator
+      * @param other The curve to use as source for the assignment
+      */
+      CurveGFp const& operator=(CurveGFp const& other);
+
+      /**
+      * Set the shared GFpModulus object.
+      * Warning: do not use this function unless you know in detail how
+      * the sharing of values
+      * in the various EC related objects works.
+      * Do NOT spread pointers to a GFpModulus over different threads!
+      * @param mod a shared pointer to a GFpModulus object suitable for
+      * *this.
+      */
+      void set_shrd_mod(std::tr1::shared_ptr<Botan::GFpModulus> const mod);
+
+      // getters
+
+      /**
+      * Get coefficient a
+      * @result coefficient a
+      */
+      GFpElement const get_a() const;
+
+      /**
+      * Get coefficient b
+      * @result coefficient b
+      */
+      GFpElement const get_b() const;
+
+      /**
+      * Get the GFpElement coefficient a  transformed
+      * to its m-residue. This can be used for efficency reasons: the curve
+      * stores the transformed version after the first invocation of this
+      * function.
+      * @result the coefficient a, transformed to its m-residue
+      */
+      GFpElement const get_mres_a() const;
+
+      /**
+      * Get the GFpElement coefficient b transformed
+      * to it´s m-residue. This can be used for efficency reasons: the curve
+      * stores the transformed version after the first invocation of this
+      * function.
+      * @result the coefficient b, transformed to it´s m-residue
+      */
+      GFpElement const get_mres_b() const;
+
+
+      /**
+      * Get the GFpElement 1  transformed
+      * to it´s m-residue. This can be used for efficency reasons: the curve
+      * stores the transformed version after the first invocation of this
+      * function.
+      * @result the GFpElement 1, transformed to it´s m-residue
+      */
+      std::tr1::shared_ptr<GFpElement const> const get_mres_one() const;
+
+      /**
+      * Get prime modulus of the field of the curve
+      * @result prime modulus of the field of the curve
+      */
+      BigInt const get_p() const;
+      /*inline std::tr1::shared_ptr<BigInt> const get_ptr_p() const
+      {
+      return mp_p;
+      }*/
+
+      /**
+      * Retrieve a shared pointer to the curves GFpModulus object for efficient storage
+      * and computation of montgomery multiplication related data members and functions.
+      * Warning: do not use this function unless you know in detail how the sharing of values
+      * in the various EC related objects works.
+      * Do NOT spread pointers to a GFpModulus over different threads!
+      * @result a shared pointer to a GFpModulus object
+      */
+      inline std::tr1::shared_ptr<Botan::GFpModulus> const get_ptr_mod() const
+         {
+         return mp_mod;
+         }
+
+      /**
+      * swaps the states of *this and other, does not throw
+      * @param other The curve to swap values with
+      */
+      void swap(CurveGFp& other);
+
+   private:
+      std::tr1::shared_ptr<Botan::GFpModulus> mp_mod;
+      GFpElement mA;
+      GFpElement mB;
+      mutable std::tr1::shared_ptr<GFpElement> mp_mres_a;
+      mutable std::tr1::shared_ptr<GFpElement> mp_mres_b;
+      mutable std::tr1::shared_ptr<GFpElement> mp_mres_one;
+   };
+
+// relational operators
+bool operator==(CurveGFp const& lhs, CurveGFp const& rhs);
+inline bool operator!=(CurveGFp const& lhs, CurveGFp const& rhs) {
+return !operator==(lhs, rhs);
+}
+
+// swaps the states of curve1 and curve2, does not throw!
+// cf. Meyers, Item 25
+inline
+void swap(CurveGFp& curve1, CurveGFp& curve2) {
+curve1.swap(curve2);
+}
+
+}
+
+namespace std {
+
+template<> inline void swap<Botan::CurveGFp>(
+   Botan::CurveGFp& curve1,
+   Botan::CurveGFp& curve2)
+   {
+   curve1.swap(curve2);
+   }
+
+}
+
+#endif
diff --git a/src/math/gfpmath/gfp_element.cpp b/src/math/gfpmath/gfp_element.cpp
new file mode 100644
index 000000000..686cc338b
--- /dev/null
+++ b/src/math/gfpmath/gfp_element.cpp
@@ -0,0 +1,674 @@
+/******************************************************
+ * Arithmetic for prime fields GF(p) (source file)    *
+ *                                                    *
+ * (C) 2007 Martin Doering                            *
+ *          [email protected]    *
+ *          Christoph Ludwig                          *
+ *          [email protected]                        *
+ *          Falko Strenzke                            *
+ *          [email protected]                    *
+ ******************************************************/
+
+#include <botan/gfp_element.h>
+#include <botan/numthry.h>
+#include <botan/def_powm.h>
+#include <botan/mp_asm.h>
+#include <botan/mp_asmi.h>
+#include <botan/mp_core.h>
+
+namespace Botan {
+
+namespace {
+
+/**
+*calculates R=b^n (here b=2) with R>m (and R beeing as small as possible) for an odd modulus m.
+* no check for oddity is performed!
+*/
+BigInt montgm_calc_r_oddmod(const BigInt& prime)
+   {
+   u32bit n = prime.sig_words();
+   BigInt result(1);
+   result <<= n*BOTAN_MP_WORD_BITS;
+   return result;
+
+   }
+
+/**
+*calculates m' with r*r^-1 - m*m' = 1
+* where r^-1 is the multiplicative inverse of r to the modulus m
+*/
+BigInt montgm_calc_m_dash(const BigInt& r, const BigInt& m, const BigInt& r_inv )
+   {
+   BigInt result = ((r * r_inv) - BigInt(1))/m;
+   return result;
+   }
+BigInt montg_trf_to_mres(const BigInt& ord_res, const BigInt& r, const BigInt& m)
+   {
+   BigInt result(ord_res);
+   result *= r;
+   result %= m;
+   return result;
+   }
+BigInt montg_trf_to_ordres(const BigInt& m_res, const BigInt& m, const BigInt& r_inv)
+   {
+   BigInt result(m_res);
+   result *= r_inv;
+   result %= m;
+   return result;
+   }
+
+void inner_montg_mult_sos(word result[], const word* a_bar,
+                          const word* b_bar, const word* n,
+                          const word* n_dash, u32bit s)
+   {
+   SecureVector<word> t;
+   t.grow_to(2*s+1);
+
+   for (u32bit i=0; i<s; i++)
+      {
+      word C = 0;
+      word S = 0;
+      for (u32bit j=0; j<s; j++)
+         {
+         S = word_madd3(a_bar[j], b_bar[i], t[i+j], &C);
+         t[i+j] = S;
+         }
+      t[i+s] = C;
+      }
+
+
+   for (u32bit i=0; i<s; i++)
+      {
+      // word word_madd2(word a, word b, word c, word* carry)
+      // returns a * b + c, resets the carry
+
+      word C = 0;
+      word m = word_madd2(t[i], n_dash[0], &C);
+      // C = 0; ???
+
+      for (u32bit j=0; j<s; j++)
+         {
+         word S = word_madd3(m, n[j], t[i+j], &C);
+         t[i+j] = S;
+         }
+
+      //// mp_mulop.cpp:
+      ////word bigint_mul_add_words(word z[], const word x[], u32bit x_size, word y)
+      u32bit cnt = 0;
+      while (C > 0)
+         {
+         // we need not worry here about C > 1, because the other operand is zero
+         word tmp = word_add(t[i+s+cnt], 0, &C);
+         t[i+s+cnt] = tmp;
+         cnt++;
+         }
+      }
+   SecureVector<word> u;
+   u.grow_to(s+1);
+   for (u32bit j=0; j<s+1; j++)
+      {
+      u[j] = t[j+s];
+      }
+
+   word B = 0;
+   word D = 0;
+   for (u32bit i=0; i<s; i++)
+      {
+      D = word_sub(u[i], n[i], &B);
+      t[i] = D;
+      }
+   D = word_sub(u[s], 0, &B);
+   t[s] = D;
+   if (B == 0)
+      {
+      for (u32bit i=0; i<s; i++)
+         {
+         result[i] = t[i];
+         }
+      }
+   else
+      {
+      for (u32bit i=0; i<s; i++)
+         {
+         result[i] = u[i];
+         }
+      }
+   }
+
+void montg_mult(BigInt& result, BigInt& a_bar, BigInt& b_bar, const BigInt& m, const BigInt& m_dash, const BigInt )
+   {
+   if (m.is_zero() || m_dash.is_zero())
+      {
+      throw Invalid_Argument("montg_mult(): neither modulus nor m_dash may be zero (and one of them was)");
+      }
+   if (a_bar.is_zero() || b_bar.is_zero())
+      {
+      result = 0;
+      }
+   u32bit s = m.sig_words();
+   a_bar.grow_to(s);
+   b_bar.grow_to(s);
+
+   result.grow_to(s);
+
+   inner_montg_mult_sos(result.get_reg(), a_bar.data(), b_bar.data(), m.data(), m_dash.data(), s);
+
+   }
+
+
+}
+
+GFpElement::GFpElement ( const BigInt& p, const BigInt& value, bool use_montgm)
+   : mp_mod(),
+     m_value(value %p),
+     m_use_montgm(use_montgm),
+     m_is_trf(false)
+   {
+   //assert(mp_mod.get() == 0);
+   mp_mod = std::tr1::shared_ptr<GFpModulus>(new GFpModulus(p));
+   //assert(mp_mod->m_p_dash == 0);
+   if (m_use_montgm)
+      {
+      ensure_montgm_precomp();
+
+      }
+
+   }
+
+GFpElement::GFpElement(std::tr1::shared_ptr<GFpModulus> const mod, const BigInt& value, bool use_montgm)
+   : mp_mod(),
+     m_value(value % mod->m_p),
+     m_use_montgm(use_montgm),
+     m_is_trf(false)
+   {
+   //assert(mp_mod.get() == 0);
+   mp_mod = mod;
+   }
+
+GFpElement::GFpElement ( GFpElement const& other )
+   : m_value ( other.m_value ),
+     m_use_montgm(other.m_use_montgm),
+     m_is_trf(other.m_is_trf)
+
+   {
+   //creates an independent copy
+   //assert((other.m_is_trf && other.m_use_montgm) || !other.m_is_trf);
+   mp_mod.reset(new GFpModulus(*other.mp_mod)); // copy-ctor of GFpModulus
+   }
+
+void GFpElement::turn_on_sp_red_mul() const
+   {
+   ensure_montgm_precomp();
+   m_use_montgm = true;
+   }
+void GFpElement::turn_off_sp_red_mul() const
+   {
+   if (m_is_trf)
+      {
+      trf_to_ordres();
+      // will happen soon anyway, so we can do it here already
+      // (this is not lazy but way more secure concerning our internal logic here)
+      }
+   m_use_montgm = false;
+   }
+void GFpElement::ensure_montgm_precomp() const
+   {
+   if ((!mp_mod->m_r.is_zero()) && (!mp_mod->m_r_inv.is_zero()) && (!mp_mod->m_p_dash.is_zero()))
+      {
+      // values are already set, nothing more to do
+      }
+   else
+      {
+      BigInt tmp_r(montgm_calc_r_oddmod(mp_mod->m_p));
+
+      BigInt tmp_r_inv(inverse_mod(tmp_r, mp_mod->m_p));
+
+      BigInt tmp_p_dash(montgm_calc_m_dash(tmp_r, mp_mod->m_p, tmp_r_inv));
+
+      mp_mod->m_r.grow_reg(tmp_r.size());
+      mp_mod->m_r_inv.grow_reg(tmp_r_inv.size());
+      mp_mod->m_p_dash.grow_reg(tmp_p_dash.size());
+
+      mp_mod->m_r = tmp_r;
+      mp_mod->m_r_inv = tmp_r_inv;
+      mp_mod->m_p_dash = tmp_p_dash;
+
+      //assert(!mp_mod->m_r.is_zero());
+      //assert(!mp_mod->m_r_inv.is_zero());
+      //assert(!mp_mod->m_p_dash.is_zero());
+      }
+
+   }
+
+void GFpElement::set_shrd_mod(std::tr1::shared_ptr<GFpModulus> const p_mod)
+   {
+   mp_mod = p_mod;
+   }
+void GFpElement::trf_to_mres() const
+   {
+   if (!m_use_montgm)
+      {
+      throw Illegal_Transformation("GFpElement is not allowed to be transformed to m-residue");
+      }
+   //assert(m_is_trf == false);
+   //assert(!mp_mod->m_r_inv.is_zero());
+   //assert(!mp_mod->m_p_dash.is_zero());
+   m_value = montg_trf_to_mres(m_value, mp_mod->m_r, mp_mod->m_p);
+   m_is_trf = true;
+   }
+void GFpElement::trf_to_ordres() const
+   {
+   //assert(m_is_trf == true);
+   m_value = montg_trf_to_ordres(m_value, mp_mod->m_p, mp_mod->m_r_inv);
+   m_is_trf = false;
+   }
+bool GFpElement::align_operands_res(GFpElement const& lhs, GFpElement const& rhs) //static
+   {
+   //assert(lhs.mp_mod->m_p == rhs.mp_mod->m_p);
+   if (lhs.m_use_montgm && rhs.m_use_montgm)
+      {
+
+      //assert(rhs.mp_mod->m_p_dash == lhs.mp_mod->m_p_dash);
+      //assert(rhs.mp_mod->m_r == lhs.mp_mod->m_r);
+      //assert(rhs.mp_mod->m_r_inv == lhs.mp_mod->m_r_inv);
+      if (!lhs.m_is_trf && !rhs.m_is_trf)
+         {
+         return false;
+         }
+      else if (lhs.m_is_trf && rhs.m_is_trf)
+         {
+         return true;
+         }
+      else // one is transf., the other not
+         {
+         if (!lhs.m_is_trf)
+            {
+            lhs.trf_to_mres();
+            //assert(rhs.m_is_trf==true);
+            return true;
+            }
+         //assert(rhs.m_is_trf==false);
+         //assert(lhs.m_is_trf==true);
+         rhs.trf_to_mres(); // the only possibility left...
+         return true;
+         }
+      }
+   else // at least one of them does not use mm
+      // (so it is impossible that both use it)
+      {
+      if (lhs.m_is_trf)
+         {
+         lhs.trf_to_ordres();
+         //assert(rhs.m_is_trf == false);
+         return false;
+         }
+      if (rhs.m_is_trf)
+         {
+         rhs.trf_to_ordres();
+         //assert(lhs.m_is_trf == false);
+         return false;
+         }
+      return false;
+      }
+   //assert(false);
+
+   }
+bool GFpElement::is_trf_to_mres() const
+   {
+   return m_is_trf;
+
+   }
+BigInt const GFpElement::get_p() const
+   {
+   return (mp_mod->m_p);
+   }
+BigInt const GFpElement::get_value() const
+   {
+
+   if (m_is_trf)
+      {
+      //assert(m_use_montgm);
+      trf_to_ordres();
+      }
+   return m_value;
+   }
+BigInt const GFpElement::get_mres() const
+   {
+   if (!m_use_montgm)
+      {
+      // does the following exception really make sense?
+      // wouldn´t it be better to simply turn on montg.mult. when
+      // this explicit request is made?
+      throw Illegal_Transformation("GFpElement is not allowed to be transformed to m-residue");
+      }
+   if (!m_is_trf)
+      {
+
+      trf_to_mres();
+
+      }
+   return m_value;
+   }
+GFpElement const& GFpElement::operator= ( GFpElement const& other )
+   {
+   m_value.grow_reg(other.m_value.size()); // grow first for exception safety
+
+   //m_value = other.m_value;
+
+   //              m_use_montgm = other.m_use_montgm;
+   //              m_is_trf = other.m_is_trf;
+   // we want to keep the member pointers, which might be part of a "sharing group"
+   // but we may not simply overwrite the BigInt values with those of the argument!!
+   // if ours already contains precomputations, it would be hazardous to
+   // set them back to zero.
+   // thus we first check for equality of the moduli,
+   // then whether either of the two objects already contains
+   // precomputed values.
+
+   // we also deal with the case were the pointers themsevles are equal:
+   if (mp_mod.get() == other.mp_mod.get())
+      {
+      // everything ok, we are in the same sharing group anyway, nothing to do
+      m_value = other.m_value; // cannot throw
+      m_use_montgm = other.m_use_montgm;
+      m_is_trf = other.m_is_trf;
+      return *this;
+      }
+   if (mp_mod->m_p != other.mp_mod->m_p)
+      {
+      // the moduli are different, this is a special case
+      // which will not occur in usual applications,
+      // so we don´t hesitate to simply create new objects
+      // (we do want to create an independent copy)
+      mp_mod.reset(new GFpModulus(*other.mp_mod)); // this could throw,
+      // and because of this
+      // we haven't modified
+      // anything so far
+      m_value = other.m_value; // can't throw
+      m_use_montgm = other.m_use_montgm;
+      m_is_trf = other.m_is_trf;
+      return *this;
+      }
+   // exception safety note: from now on we are on the safe
+   // side with respect to the modulus,
+   // so we can assign the value now:
+   m_value = other.m_value;
+   m_use_montgm = other.m_use_montgm;
+   m_is_trf = other.m_is_trf;
+   // the moduli are equal, but we deal with different sharing groups.
+   // we will NOT fuse the sharing goups
+   // and we will NOT reset already precomputed values
+   if (mp_mod->has_precomputations())
+      {
+      // our own sharing group already has precomputed values,
+      // so nothing to do.
+      return *this;
+      }
+   else
+      {
+      // let´s see whether the argument has something for us...
+      if (other.mp_mod->has_precomputations())
+         {
+         // fetch them for our sharing group
+         // exc. safety note: grow first
+         mp_mod->m_p_dash.grow_reg(other.mp_mod->m_p_dash.size());
+         mp_mod->m_r.grow_reg(other.mp_mod->m_r.size());
+         mp_mod->m_r_inv.grow_reg(other.mp_mod->m_r_inv.size());
+
+         mp_mod->m_p_dash = other.mp_mod->m_p_dash;
+         mp_mod->m_r = other.mp_mod->m_r;
+         mp_mod->m_r_inv = other.mp_mod->m_r_inv;
+         return *this;
+         }
+      }
+   // our precomputations aren´t set, the arguments neither,
+   // so we let them alone
+   return *this;
+
+
+   }
+void GFpElement::share_assign(GFpElement const& other)
+   {
+   //assert((other.m_is_trf && other.m_use_montgm) || !other.m_is_trf);
+
+   // use grow_to to make it exc safe
+   m_value.grow_reg(other.m_value.size());
+   m_value = other.m_value;
+
+   m_use_montgm = other.m_use_montgm;
+   m_is_trf = other.m_is_trf;
+   mp_mod = other.mp_mod; // cannot throw
+
+   }
+GFpElement& GFpElement::operator+= ( GFpElement const& rhs )
+   {
+   GFpElement::align_operands_res(*this, rhs);
+
+   workspace = m_value;
+   workspace += rhs.m_value;
+   if (workspace >= mp_mod->m_p)
+      {
+      workspace -= mp_mod->m_p;
+      }
+
+   m_value = workspace;
+   //assert(m_value < mp_mod->m_p);
+   //assert(m_value >= 0);
+
+   return *this;
+   }
+GFpElement& GFpElement::operator-= ( GFpElement const& rhs )
+   {
+   GFpElement::align_operands_res(*this, rhs);
+
+   workspace = m_value;
+
+   workspace -= rhs.m_value;
+
+   if (workspace.is_negative())
+      {
+      workspace += mp_mod->m_p;
+      }
+
+   m_value = workspace;
+   //assert(m_value < mp_mod->m_p);
+   //assert(m_value >= 0);
+   return *this;
+
+   }
+GFpElement& GFpElement::operator*= (u32bit rhs)
+   {
+   workspace = m_value;
+   workspace *= rhs;
+   workspace %= mp_mod->m_p;
+   m_value = workspace;
+   return *this;
+   }
+GFpElement& GFpElement::operator*= ( GFpElement const& rhs )
+   {
+   //assert(rhs.mp_mod->m_p == mp_mod->m_p);
+   // here, we do not use align_operands_res() for one simple reason:
+   // we want to enforce the transformation to an m-residue, otherwise it would
+   // never happen
+   if (m_use_montgm && rhs.m_use_montgm)
+      {
+      //assert(rhs.mp_mod->m_p == mp_mod->m_p); // is montgm. mult is on, then precomps must be there
+      //assert(rhs.mp_mod->m_p_dash == mp_mod->m_p_dash);
+      //assert(rhs.mp_mod->m_r == mp_mod->m_r);
+      if (!m_is_trf)
+         {
+         trf_to_mres();
+         }
+      if (!rhs.m_is_trf)
+         {
+         rhs.trf_to_mres();
+         }
+      workspace = m_value;
+      montg_mult(m_value, workspace, rhs.m_value, mp_mod->m_p, mp_mod->m_p_dash, mp_mod->m_r);
+      }
+   else // ordinary multiplication
+      {
+      if (m_is_trf)
+         {
+         //assert(m_use_montgm);
+         trf_to_ordres();
+         }
+      if (rhs.m_is_trf)
+         {
+         //assert(rhs.m_use_montgm);
+         rhs.trf_to_ordres();
+         }
+      workspace = m_value;
+      workspace *= rhs.m_value;
+      workspace %= mp_mod->m_p;
+      m_value = workspace;
+      }
+   return *this;
+   }
+GFpElement& GFpElement::operator/= ( GFpElement const& rhs )
+   {
+   bool use_mres = GFpElement::align_operands_res(*this, rhs);
+   //assert((this->m_is_trf && rhs.m_is_trf) || !(this->m_is_trf && rhs.m_is_trf) );
+   // (internal note: see C86)
+   if (use_mres)
+      {
+      //assert(m_use_montgm && rhs.m_use_montgm);
+      GFpElement rhs_ordres(rhs);
+      rhs_ordres.trf_to_ordres();
+      rhs_ordres.inverse_in_place();
+      workspace = m_value;
+      workspace *=  rhs_ordres.get_value();
+      workspace %= mp_mod->m_p;
+      m_value = workspace;
+
+      }
+   else
+      {
+      GFpElement inv_rhs(rhs);
+      inv_rhs.inverse_in_place();
+      *this *= inv_rhs;
+      }
+   return *this;
+   }
+bool GFpElement::is_zero()
+   {
+   return (m_value.is_zero());
+   // this is correct because x_bar = x * r = x = 0 for x = 0
+   }
+GFpElement& GFpElement::inverse_in_place()
+   {
+   m_value = inverse_mod(m_value, mp_mod->m_p);
+   if (m_is_trf)
+      {
+      //assert(m_use_montgm);
+
+      m_value *= mp_mod->m_r;
+      m_value *= mp_mod->m_r;
+      m_value %= mp_mod->m_p;
+      }
+   //assert(m_value <= mp_mod->m_p);
+   return *this;
+
+   }
+GFpElement& GFpElement::negate()
+   {
+   m_value = mp_mod->m_p - m_value;
+   //assert(m_value <= mp_mod->m_p);
+   return *this;
+   }
+void GFpElement::swap ( GFpElement& other )
+   {
+   m_value.swap ( other.m_value );
+   mp_mod.swap(other.mp_mod);
+   std::swap<bool>(m_use_montgm,other.m_use_montgm);
+   std::swap<bool>(m_is_trf,other.m_is_trf);
+   }
+
+bool operator== ( GFpElement const& lhs, GFpElement const& rhs )
+   {
+   // for effeciency reasons we firstly check whether
+   //the modulus pointers are different in the first place:
+   if (lhs.get_ptr_mod() != rhs.get_ptr_mod())
+      {
+      if (lhs.get_p() != rhs.get_p())
+         {
+         return false;
+         }
+      }
+   // so the modulus is equal, now check the values
+   bool use_mres = GFpElement::align_operands_res(lhs, rhs);
+
+   if (use_mres)
+      {
+      return (lhs.get_mres() == rhs.get_mres());
+      }
+   else
+      {
+      return(lhs.get_value() == rhs.get_value());
+      }
+   }
+
+GFpElement operator+ ( GFpElement const& lhs, GFpElement const& rhs )
+   {
+   // consider the case that lhs and rhs both use montgm:
+   // then += returns an element which uses montgm.
+   // thus the return value of op+ here will be an element
+   // using montgm in this case
+   // NOTE: the rhs might be transformed when using op+, the lhs never
+   GFpElement result ( lhs );
+   result += rhs;
+   return result;
+   }
+GFpElement operator- ( GFpElement const& lhs, GFpElement const& rhs )
+   {
+   GFpElement result ( lhs );
+   result -= rhs;
+   return result;
+   // NOTE: the rhs might be transformed when using op-, the lhs never
+   }
+GFpElement operator- ( GFpElement const& lhs )
+   {
+   return ( GFpElement ( lhs ) ).negate();
+   }
+GFpElement operator* ( GFpElement const& lhs, GFpElement const& rhs )
+   {
+   // consider the case that lhs and rhs both use montgm:
+   // then *= returns an element which uses montgm.
+   // thus the return value of op* here will be an element
+   // using montgm in this case
+   GFpElement result ( lhs );
+   result *= rhs;
+   return result;
+   }
+GFpElement operator*(GFpElement const& lhs, u32bit rhs)
+   {
+   GFpElement result(lhs);
+   result *= rhs;
+   return result;
+
+   }
+GFpElement operator*(u32bit lhs, GFpElement const& rhs )
+   {
+   return rhs*lhs;
+   }
+GFpElement operator/ ( GFpElement const& lhs, GFpElement const& rhs )
+   {
+   GFpElement result (lhs);
+   result /= rhs;
+   return result;
+   }
+SecureVector<byte> FE2OSP ( GFpElement const& elem )
+   {
+   return BigInt::encode_1363 ( elem.get_value(), elem.get_p().bytes() );
+   }
+GFpElement OS2FEP ( MemoryRegion<byte> const& os, BigInt p )
+   {
+
+   return GFpElement ( p, BigInt::decode ( os.begin(), os.size() ) );
+   }
+GFpElement inverse ( GFpElement const& elem )
+   {
+   return GFpElement ( elem ).inverse_in_place();
+   }
+
+} // namespace Botan
diff --git a/src/math/gfpmath/gfp_element.h b/src/math/gfpmath/gfp_element.h
new file mode 100644
index 000000000..e9850df30
--- /dev/null
+++ b/src/math/gfpmath/gfp_element.h
@@ -0,0 +1,308 @@
+/******************************************************
+ * Arithmetic for prime fields GF(p) (header file)    *
+ *                                                    *
+ * (C) 2007 Martin Döring                             *
+*          [email protected]    *
+*          Christoph Ludwig                          *
+*          [email protected]                        *
+*          Falko Strenzke                            *
+*          [email protected]                    *
+******************************************************/
+
+#ifndef BOTAN_MATH_GF_GFP_ELEMENT_H_GUARD_
+#define BOTAN_MATH_GF_GFP_ELEMENT_H_GUARD_
+
+#include <botan/gfp_modulus.h>
+#include <botan/bigint.h>
+#include <tr1/memory>
+
+namespace Botan
+{
+
+struct Illegal_Transformation : public Exception
+{
+    Illegal_Transformation(const std::string& err = "Requested transformation is not possible")
+    : Exception(err) {}
+};
+
+/**
+* This class represents one element in GF(p). Enables the convenient, transparent use
+* of the montgomery multiplication.
+*/
+class GFpElement
+   {
+
+   private:
+      std::tr1::shared_ptr<GFpModulus> mp_mod;
+      mutable BigInt m_value; // ordinary residue or m-residue respectively
+      mutable BigInt workspace;
+      // *****************************************
+      // data members for montgomery multiplication
+      mutable bool m_use_montgm;
+      //mutable BigInt m_mres;
+      // this bool tells use whether the m_mres carries
+      // the actual value (in this case mValue doesn´t)
+      mutable bool m_is_trf;
+
+
+      void ensure_montgm_precomp() const;
+      void trf_to_mres() const;
+      void trf_to_ordres() const;
+
+   public:
+
+
+      /** construct an element of GF(p) with the given value.
+      * use_montg defaults to false and determines wether Montgomery multiplications
+      * will be use when applying operators '*' , '*='.
+      * @param p the prime number of the field
+      * @param value the element value
+      * @param use_montgm whether this object will use Montgomery multiplication
+      */
+      explicit GFpElement ( const BigInt& p, const BigInt& value, bool use_montgm = false );
+
+
+      /** construct an element of GF(p) with the given value (defaults to 0).
+      * use_montg defaults to false and determines wether montgomery multiplications
+      * will be use when applying operators '*' , '*='.
+      * Use this constructor for efficient use of Montgomery multiplication in a context with a
+      * fixed a modulus.
+      * Warning: do not use this function unless you know in detail about
+      * the implications of using
+      * the shared GFpModulus objects!
+      * @param mod shared pointer to the GFpModulus to be shared
+      * @param value the element value
+      * @param use_montgm whether this object will use Montgomery multiplication
+      */
+      explicit GFpElement(std::tr1::shared_ptr<GFpModulus> const mod, const BigInt& value, bool use_mongm = false);
+
+      /**
+      * Copy constructor
+      * @param other The element to clone
+      */
+      GFpElement ( GFpElement const& other );
+
+      /**
+      * Assignment operator.
+      * makes *this a totally independent object
+      * (gives *this independent modulus specific values).
+      *
+      * @param other The element to assign to our object
+      */
+      GFpElement const& operator= ( GFpElement const& other );
+
+      /**
+      * Works like the assignment operator, but lets
+      * *this share the modulus dependend value with other.
+      * Warning: do not use this function unless you know in detail about
+      * the implications of using
+      * the shared GFpModulus objects!
+      * @param other The element to assign to our object
+      */
+      void share_assign(GFpElement const& other);
+
+      /**
+      * Switch Montgomery multiplcation optimizations ON
+      */
+      void turn_on_sp_red_mul() const;
+
+      /**
+      * Switch Montgomery multiplcation optimizations OFF
+      */
+      void turn_off_sp_red_mul() const;
+
+      /**
+      * += Operator
+      * @param rhs the GFpElement to add to the local value
+      * @result *this
+      */
+      GFpElement& operator+= ( GFpElement const& rhs );
+
+      /**
+      * -= Operator
+      * @param rhs the GFpElement to subtract from the local value
+      * @result *this
+      */
+      GFpElement& operator-= ( GFpElement const& rhs );
+
+      /**
+      * *= Operator
+      * @param rhs the GFpElement to multiply with the local value
+      * @result *this
+      */
+      GFpElement& operator*= ( GFpElement const& rhs );
+      /**
+      * /= Operator
+      * @param rhs the GFpElement to divide the local value by
+      * @result *this
+      */
+      GFpElement& operator/= ( GFpElement const& rhs );
+
+      /**
+      * *= Operator
+      * @param rhs the value to multiply with the local value
+      * @result *this
+      */
+      GFpElement& operator*= (u32bit rhs);
+
+      /**
+      * Negate internal value ( *this *= -1 )
+      * @return *this
+      */
+      GFpElement& negate();
+
+      /**
+      * Assigns the inverse of *this to *this, i.e.
+      * *this = (*this)^(-1)
+      * @result *this
+      */
+      GFpElement& inverse_in_place();
+
+      /**
+      * checks whether the value is zero (without provoking
+      * a backtransformation to the ordinary-residue)
+      * @result true, if the value is zero, false otherwise.
+      */
+      bool is_zero();
+
+      /**
+      * return prime number of GF(p)
+      * @result a prime number
+      */
+      BigInt const get_p() const;
+
+      /**
+      * Return the represented value in GF(p)
+      * @result The value in GF(p)
+      */
+      BigInt const get_value() const;
+
+      /**
+      * Returns the shared pointer to the GFpModulus of *this.
+      * Warning: do not use this function unless you know in detail about
+      * the implications of using
+      * the shared GFpModulus objects!
+      * @result the shared pointer to the GFpModulus of *this
+      */
+      inline std::tr1::shared_ptr<GFpModulus> const get_ptr_mod() const
+         {
+         return mp_mod;
+         }
+
+
+      /**
+      * Sets the shared pointer to the GFpModulus of *this.
+      * Warning: do not use this function unless you know in detail about
+      * the implications of using
+      * the shared GFpModulus objects!
+      * @param mod a shared pointer to a GFpModulus that will be held in *this
+      */
+      void set_shrd_mod(std::tr1::shared_ptr<GFpModulus> const mod);
+
+      /**
+      * Tells whether this GFpElement is currently transformed to it´ m-residue,
+      * i.e. in the form x_bar = x * r mod m.
+      * @result true if it is currently transformed to it´s m-residue.
+      */
+      bool is_trf_to_mres() const;
+
+      /**
+      * Transforms this to x_bar = x * r mod m
+      * @result return the value x_bar.
+      */
+      BigInt const get_mres() const;
+
+      /**
+      * Check, if montgomery multiplication is used.
+      * @result true, if montgomery multiplication is used, false otherwise
+      */
+      bool is_use_montgm() const
+         {
+         return m_use_montgm;
+         }
+
+      /**
+      * Transforms the arguments in such way that either both
+      * are in m-residue representation (returns true) or both are
+      * in ordinary residue representation (returns false).
+      * m-residue is prefered in case of ambiguity.
+      * does not toggle m_use_montgm of the arguments.
+      * Don´t be confused about the constness of the arguments:
+      * the transformation between normal residue and m-residue is
+      * considered as leaving the object const.
+      * @param lhs the first operand to be aligned
+      * @param rhs the second operand to be aligned
+      * @result true if both are transformed to their m-residue,
+      * false it both are transformed to their normal residue.
+      */
+      static bool align_operands_res(GFpElement const& lhs, GFpElement const& rhs);
+
+      //friend declarations for non-member functions
+
+      /**
+      * write a GFpElement to an output stream.
+      * @param output the output stream to write to
+      * @param elem the object to write
+      * @result the output stream
+      */
+      friend class Point_Coords_GFp;
+
+      /**
+      * swaps the states of *this and other, does not throw!
+      * @param other The value to swap with
+      */
+      void swap ( GFpElement& other );
+
+   };
+
+// relational operators
+bool operator== ( GFpElement const& lhs, GFpElement const& rhs );
+inline bool operator!= ( GFpElement const& lhs, GFpElement const& rhs )
+   {
+   return !operator== ( lhs, rhs );
+   }
+
+// arithmetic operators
+GFpElement operator+ ( GFpElement const& lhs, GFpElement const& rhs );
+GFpElement operator- ( GFpElement const& lhs, GFpElement const& rhs );
+GFpElement operator- ( GFpElement const& lhs );
+
+GFpElement operator* ( GFpElement const& lhs, GFpElement const& rhs );
+GFpElement operator/ ( GFpElement const& lhs, GFpElement const& rhs );
+GFpElement operator* (GFpElement const& lhs, u32bit rhs);
+GFpElement operator* (u32bit rhs, GFpElement const& lhs);
+
+// return (*this)^(-1)
+GFpElement inverse ( GFpElement const& elem );
+
+// encoding and decoding
+SecureVector<byte> FE2OSP ( GFpElement const& elem );
+GFpElement OS2FEP ( MemoryRegion<byte> const& os, BigInt p );
+
+
+// swaps the states of elem1 and elem2, does not throw!
+// cf. Meyers, Item 25
+inline
+void swap ( GFpElement& elem1, GFpElement& elem2 )
+   {
+   elem1.swap ( elem2 );
+   }
+
+} // namespace Botan
+
+namespace std
+{
+
+// swaps the states of elem1 and elem2, does not throw!
+// cf. Meyers, Item 25
+template<>
+inline
+void swap< Botan::GFpElement>(Botan::GFpElement& elem1,
+                              Botan::GFpElement& elem2)
+   {
+   elem1.swap(elem2);
+   }
+
+} // namespace std
+
+#endif
diff --git a/src/math/gfpmath/gfp_modulus.h b/src/math/gfpmath/gfp_modulus.h
new file mode 100644
index 000000000..5edf44ba0
--- /dev/null
+++ b/src/math/gfpmath/gfp_modulus.h
@@ -0,0 +1,124 @@
+/******************************************************
+ * Modulus and related data for a specific            *
+ * implementation of  GF(p) (header file)             *
+ *                                                    *
+ * (C) 2008 Martin Döring                             *
+ *          [email protected]    *
+ *          Christoph Ludwig                          *
+ *          [email protected]                        *
+ *          Falko Strenzke                            *
+ *          [email protected]                    *
+ ******************************************************/
+
+#ifndef BOTAN_MATH_GF_GFP_MODULUS_H_GUARD_
+#define BOTAN_MATH_GF_GFP_MODULUS_H_GUARD_
+
+#include <botan/bigint.h>
+
+namespace Botan
+{
+
+class GFpElement;
+/**
+* This class represents a GFpElement modulus including the modulus related
+* values necessary for the montgomery multiplication.
+*/
+class GFpModulus
+   {
+      friend class GFpElement;
+   private:
+      BigInt m_p; // the modulus itself
+      mutable BigInt m_p_dash;
+      mutable BigInt m_r;
+      mutable BigInt m_r_inv;
+   public:
+
+      /**
+      * Construct a GF(P)-Modulus from a BigInt
+      */
+      GFpModulus(BigInt p)
+         : m_p(p),
+           m_p_dash(),
+           m_r(),
+           m_r_inv()
+         {}
+
+      /**
+      * Tells whether the precomputations necessary for the use of the montgomery
+      * multiplication have yet been established.
+      * @result true if the precomputated value are already available.
+      */
+      inline bool has_precomputations() const
+         {
+         return(!m_p_dash.is_zero() && !m_r.is_zero() && !m_r_inv.is_zero());
+         }
+
+      /**
+      * Swaps this with another GFpModulus, does not throw.
+      * @param other the GFpModulus to swap *this with.
+      */
+      inline void swap(GFpModulus& other)
+         {
+         m_p.swap(other.m_p);
+         m_p_dash.swap(other.m_p_dash);
+         m_r.swap(other.m_r);
+         m_r_inv.swap(other.m_r_inv);
+         }
+
+      /**
+      * Tells whether the modulus of *this is equal to the argument.
+      * @param mod the modulus to compare this with
+      * @result true if the modulus of *this and the argument are equal.
+      */
+      inline bool p_equal_to(const BigInt& mod) const
+         {
+         return (m_p == mod);
+         }
+
+      /**
+      * Return the modulus of this GFpModulus.
+      * @result the modulus of *this.
+      */
+      inline const BigInt get_p() const
+         {
+         return m_p;
+         }
+
+      /**
+      * returns the montgomery multiplication related value r.
+      * Warning: will be zero if precomputations have not yet been
+      * performed!
+      * @result r
+      */
+      inline const BigInt get_r() const
+         {
+         return m_r;
+         }
+
+      /**
+      * returns the montgomery multiplication related value r^{-1}.
+      * Warning: will be zero if precomputations have not yet been
+      * performed!
+      * @result r^{-1}
+      */
+      inline const BigInt get_r_inv() const
+         {
+         return m_r_inv;
+         }
+
+      /**
+      * returns the montgomery multiplication related value p'.
+      * Warning: will be zero if precomputations have not yet been
+      * performed!
+      * @result p'
+      */
+      inline const BigInt get_p_dash() const
+         {
+         return m_p_dash;
+         }
+      // default cp-ctor, op= are fine
+   };
+
+}
+
+#endif
diff --git a/src/math/gfpmath/point_gfp.cpp b/src/math/gfpmath/point_gfp.cpp
new file mode 100644
index 000000000..23b6d4518
--- /dev/null
+++ b/src/math/gfpmath/point_gfp.cpp
@@ -0,0 +1,1157 @@
+/******************************************************
+ * Arithmetic for point groups of elliptic curves     *
+ * over GF(p) (source file)                           *
+ *                                                    *
+ * (C) 2007 Martin Döring                             *
+ *          [email protected]    *
+ *          Christoph Ludwig                          *
+ *          [email protected]                        *
+ *          Falko Strenzke                            *
+ *          [email protected]                    *
+ ******************************************************/
+
+#include <botan/point_gfp.h>
+#include <botan/numthry.h>
+
+namespace Botan {
+
+// construct the point at infinity or a random point
+PointGFp::PointGFp(CurveGFp const& curve)
+   :	mC(curve),
+        mX(curve.get_p(), 0),
+        mY(curve.get_p(), 1),
+        mZ(curve.get_p(), 0),
+        mZpow2(curve.get_p(),0),
+        mZpow3(curve.get_p(),0),
+        mAZpow4(curve.get_p(),0),
+        mZpow2_set(false),
+        mZpow3_set(false),
+        mAZpow4_set(false)
+   {
+   // first set the point wide pointer
+
+   set_shrd_mod(mC.get_ptr_mod());
+
+   }
+
+
+// construct a point given its jacobian projective coordinates
+PointGFp::PointGFp(CurveGFp const& curve, GFpElement const& x,
+                   GFpElement const& y, GFpElement const& z)
+   :	mC(curve),
+        mX(x),
+        mY(y),
+        mZ(z),
+        mZpow2(curve.get_p(),0),
+        mZpow3(curve.get_p(),0),
+        mAZpow4(curve.get_p(),0),
+        mZpow2_set(false),
+        mZpow3_set(false),
+        mAZpow4_set(false)
+   {
+   set_shrd_mod(mC.get_ptr_mod());
+   }
+PointGFp::PointGFp ( CurveGFp const& curve, GFpElement const& x,
+                     GFpElement const& y )
+   :mC(curve),
+    mX(x),
+    mY(y),
+    mZ(curve.get_p(),1),
+    mZpow2(curve.get_p(),0),
+    mZpow3(curve.get_p(),0),
+    mAZpow4(curve.get_p(),0),
+    mZpow2_set(false),
+    mZpow3_set(false),
+    mAZpow4_set(false)
+   {
+   set_shrd_mod(mC.get_ptr_mod());
+   }
+
+// copy constructor
+PointGFp::PointGFp(PointGFp const& other)
+   :	mC(other.mC),
+        mX(other.mX),
+        mY(other.mY),
+        mZ(other.mZ),
+        mZpow2(other.mZpow2),
+        mZpow3(other.mZpow3),
+        mAZpow4(other.mAZpow4),
+        mZpow2_set(other.mZpow2_set),
+        mZpow3_set(other.mZpow3_set),
+        mAZpow4_set(other.mAZpow4_set)
+   {
+   set_shrd_mod(mC.get_ptr_mod());
+   }
+
+// assignment operator
+PointGFp const& PointGFp::operator=(PointGFp const& other)
+   {
+   mC = other.get_curve();
+   mX = other.get_jac_proj_x();
+   mY = other.get_jac_proj_y();
+   mZ = other.get_jac_proj_z();
+   mZpow2 = GFpElement(other.mZpow2);
+   mZpow3 = GFpElement(other.mZpow3);
+   mAZpow4 = GFpElement(other.mAZpow4);
+   mZpow2_set = other.mZpow2_set;
+   mZpow3_set = other.mZpow3_set;
+   mAZpow4_set = other.mAZpow4_set;
+   set_shrd_mod(mC.get_ptr_mod());
+   return *this;
+   }
+
+PointGFp const& PointGFp::assign_within_same_curve(PointGFp const& other)
+   {
+   mX = other.get_jac_proj_x();
+   mY = other.get_jac_proj_y();
+   mZ = other.get_jac_proj_z();
+   mZpow2_set = false;
+   mZpow3_set = false;
+   mAZpow4_set = false;
+   // the rest stays!
+   return *this;
+   }
+
+void PointGFp::set_shrd_mod(std::tr1::shared_ptr<GFpModulus> p_mod)
+   {
+   mX.set_shrd_mod(p_mod);
+   mY.set_shrd_mod(p_mod);
+   mZ.set_shrd_mod(p_mod);
+   mZpow2.set_shrd_mod(p_mod);
+   mZpow3.set_shrd_mod(p_mod);
+   mAZpow4.set_shrd_mod(p_mod);
+   }
+
+void PointGFp::ensure_worksp() const
+   {
+   if (mp_worksp_gfp_el.get() != 0)
+      {
+      if ((*mp_worksp_gfp_el).size() == GFPEL_WKSP_SIZE)
+         {
+         return;
+         }
+      else
+         {
+         throw Invalid_State("encountered incorrect size for PointGFp´s GFpElement workspace");
+         }
+      }
+
+   mp_worksp_gfp_el = std::tr1::shared_ptr<std::vector<GFpElement> >(new std::vector<GFpElement>);
+   mp_worksp_gfp_el->reserve(9);
+   for (u32bit i=0; i<GFPEL_WKSP_SIZE; i++)
+      {
+      mp_worksp_gfp_el->push_back(GFpElement(1,0));
+
+      }
+   }
+
+// arithmetic operators
+PointGFp& PointGFp::operator+=(PointGFp const& rhs)
+   {
+   if (is_zero())
+      {
+      *this = rhs;
+      return *this;
+      }
+   if (rhs.is_zero())
+      {
+      return *this;
+      }
+   ensure_worksp();
+
+   if (rhs.mZ == *(mC.get_mres_one()))
+      {
+      //U1 = mX;
+      (*mp_worksp_gfp_el)[0].share_assign(mX);
+
+      //S1 = mY;
+      (*mp_worksp_gfp_el)[2].share_assign(mY);
+      }
+   else
+      {
+      if ((!rhs.mZpow2_set) || (!rhs.mZpow3_set))
+         {
+         rhs.mZpow2 = rhs.mZ;
+         rhs.mZpow2 *= rhs.mZ;
+         rhs.mZpow3 = rhs.mZpow2;
+         rhs.mZpow3 *= rhs.mZ;
+
+         rhs.mZpow2_set = true;
+         rhs.mZpow3_set = true;
+         }
+      //U1 = mX * rhs.mZpow2;
+      (*mp_worksp_gfp_el)[0].share_assign(mX);
+      (*mp_worksp_gfp_el)[0] *= rhs.mZpow2;
+
+      //S1 = mY * rhs.mZpow3;
+      (*mp_worksp_gfp_el)[2].share_assign(mY);
+      (*mp_worksp_gfp_el)[2] *= rhs.mZpow3;
+
+      }
+   if (mZ == *(mC.get_mres_one()))
+      {
+      //U2 = rhs.mX;
+      (*mp_worksp_gfp_el)[1].share_assign(rhs.mX);
+
+      //S2 = rhs.mY;
+      (*mp_worksp_gfp_el)[3].share_assign(rhs.mY);
+      }
+   else
+      {
+      if ((!mZpow2_set) || (!mZpow3_set))
+         {
+         // precomputation can´t be used, because *this changes anyway
+         mZpow2 = mZ;
+         mZpow2 *= mZ;
+
+         mZpow3 = mZpow2;
+         mZpow3 *= mZ;
+         }
+      //U2 = rhs.mX * mZpow2;
+      (*mp_worksp_gfp_el)[1].share_assign(rhs.mX);
+      (*mp_worksp_gfp_el)[1] *= mZpow2;
+
+      //S2 = rhs.mY * mZpow3;
+      (*mp_worksp_gfp_el)[3].share_assign(rhs.mY);
+      (*mp_worksp_gfp_el)[3] *= mZpow3;
+
+      }
+   //GFpElement H(U2 - U1);
+
+   (*mp_worksp_gfp_el)[4].share_assign((*mp_worksp_gfp_el)[1]);
+   (*mp_worksp_gfp_el)[4] -= (*mp_worksp_gfp_el)[0];
+
+   //GFpElement r(S2 - S1);
+   (*mp_worksp_gfp_el)[5].share_assign((*mp_worksp_gfp_el)[3]);
+   (*mp_worksp_gfp_el)[5] -= (*mp_worksp_gfp_el)[2];
+
+   //if(H.is_zero())
+   if ((*mp_worksp_gfp_el)[4].is_zero())
+
+      {
+      if ((*mp_worksp_gfp_el)[5].is_zero())
+
+         {
+         mult2_in_place();
+         return *this;
+         }
+      *this = PointGFp(mC); // setting myself to zero
+      return *this;
+      }
+
+   //U2 = H * H;
+   (*mp_worksp_gfp_el)[1].share_assign((*mp_worksp_gfp_el)[4]);
+   (*mp_worksp_gfp_el)[1] *= (*mp_worksp_gfp_el)[4];
+
+   //S2 = U2 * H;
+   (*mp_worksp_gfp_el)[3].share_assign((*mp_worksp_gfp_el)[1]);
+   (*mp_worksp_gfp_el)[3] *= (*mp_worksp_gfp_el)[4];
+
+   //U2 *= U1;
+   (*mp_worksp_gfp_el)[1] *= (*mp_worksp_gfp_el)[0];
+
+   //GFpElement x(r*r - S2 - (U2+U2));
+   (*mp_worksp_gfp_el)[6].share_assign((*mp_worksp_gfp_el)[5]);
+   (*mp_worksp_gfp_el)[6] *= (*mp_worksp_gfp_el)[5];
+   (*mp_worksp_gfp_el)[6] -= (*mp_worksp_gfp_el)[3];
+   (*mp_worksp_gfp_el)[6] -= (*mp_worksp_gfp_el)[1];
+   (*mp_worksp_gfp_el)[6] -= (*mp_worksp_gfp_el)[1];
+
+   //GFpElement z(S1 * S2);
+   (*mp_worksp_gfp_el)[8].share_assign((*mp_worksp_gfp_el)[2]);
+   (*mp_worksp_gfp_el)[8] *= (*mp_worksp_gfp_el)[3];
+
+   //GFpElement y(r * (U2-x) - z);
+   (*mp_worksp_gfp_el)[7].share_assign((*mp_worksp_gfp_el)[1]);
+   (*mp_worksp_gfp_el)[7] -= (*mp_worksp_gfp_el)[6];
+   (*mp_worksp_gfp_el)[7] *= (*mp_worksp_gfp_el)[5];
+   (*mp_worksp_gfp_el)[7] -= (*mp_worksp_gfp_el)[8];
+
+   if (mZ == *(mC.get_mres_one()))
+      {
+      if (rhs.mZ != *(mC.get_mres_one()))
+         {
+         //z = rhs.mZ * H;
+         (*mp_worksp_gfp_el)[8].share_assign(rhs.mZ);
+         (*mp_worksp_gfp_el)[8] *= (*mp_worksp_gfp_el)[4];
+         }
+      else
+         {
+         //z = H;
+         (*mp_worksp_gfp_el)[8].share_assign((*mp_worksp_gfp_el)[4]);
+         }
+      }
+   else if (rhs.mZ != *(mC.get_mres_one()))
+      {
+      //U1 = mZ * rhs.mZ;
+      (*mp_worksp_gfp_el)[0].share_assign(mZ);
+      (*mp_worksp_gfp_el)[0] *= rhs.mZ;
+
+      //z = U1 * H;
+      (*mp_worksp_gfp_el)[8].share_assign((*mp_worksp_gfp_el)[0]);
+      (*mp_worksp_gfp_el)[8] *= (*mp_worksp_gfp_el)[4];
+
+      }
+   else
+      {
+      //z = mZ * H;
+      (*mp_worksp_gfp_el)[8].share_assign(mZ);
+      (*mp_worksp_gfp_el)[8] *= (*mp_worksp_gfp_el)[4];
+
+      }
+   mZpow2_set = false;
+   mZpow3_set = false;
+   mAZpow4_set = false;
+
+   mX = (*mp_worksp_gfp_el)[6];
+   mY = (*mp_worksp_gfp_el)[7];
+   mZ = (*mp_worksp_gfp_el)[8];
+
+   return *this;
+
+   }
+PointGFp& PointGFp::operator-=(PointGFp const& rhs)
+   {
+   PointGFp minus_rhs = PointGFp(rhs).negate();
+
+   if (is_zero())
+      {
+      *this = minus_rhs;
+      }
+   else
+      {
+      *this += minus_rhs;
+      }
+   return *this;
+   }
+
+PointGFp& PointGFp::mult_this_secure(const BigInt& scalar,
+                                     const BigInt& /*point_order*/,
+                                     const BigInt& /*max_secr*/)
+   {
+   // NOTE: FS: so far this is code duplication of op*=.
+   // we have to see how we deal with this.
+   // fact is that we will probably modify this function
+   // while evaluating the countermeasures
+   // whereas we probably will not start modifying the
+   // function operator*=.
+   // however, in the end both should be merged.
+
+   // use montgomery mult. in this operation
+   this->turn_on_sp_red_mul();
+
+   std::tr1::shared_ptr<PointGFp> H(new PointGFp(this->mC));
+   std::tr1::shared_ptr<PointGFp> tmp; // used for AADA
+
+   PointGFp P(*this);
+   BigInt m(scalar);
+
+   if (m < BigInt(0))
+      {
+      m = -m;
+      P.negate();
+      }
+   if (P.is_zero() || (m == BigInt(0)))
+      {
+      *this = *H;
+      return *this;
+      }
+   if (m == BigInt(1))
+      {
+      return *this;
+      }
+   //
+#ifdef CM_AADA
+#ifndef CM_RAND_EXP
+   int max_secr_bits = max_secr.bits();
+#endif
+#endif
+
+   int mul_bits = m.bits(); // this is used for a determined number of loop runs in
+   // the mult_loop where leading zero´s are padded if necessary.
+   // Here we assign the value that will be used when no countermeasures are specified
+#ifdef CM_RAND_EXP
+   u32bit rand_r_bit_len = 20; // Coron(99) proposes 20 bit for r
+
+#ifdef CM_AADA
+
+   BigInt r_max(1);
+
+#endif // CM_AADA
+
+   // use randomized exponent
+#ifdef TA_COLL_T
+   static BigInt r_randexp;
+   if (new_rand)
+      {
+      r_randexp = random_integer(rand_r_bit_len);
+      }
+   //assert(!r_randexp.is_zero());
+#else
+   BigInt r_randexp(random_integer(rand_r_bit_len));
+#endif
+
+   m += r_randexp * point_order;
+   // determine mul_bits...
+#ifdef CM_AADA
+   // AADA with rand. Exp.
+   //assert(rand_r_bit_len > 0);
+   r_max <<= rand_r_bit_len;
+   r_max -= 1;
+   //assert(r_max.bits() == rand_r_bit_len);
+   mul_bits = (max_secr + point_order * r_max).bits();
+#else
+   // rand. Exp. without AADA
+   mul_bits = m.bits();
+#endif // CM_AADA
+
+
+#endif // CM_RAND_EXP
+
+   // determine mul_bits...
+#if (CM_AADA == 1 && CM_RAND_EXP != 1)
+
+   mul_bits = max_secr_bits;
+#endif // CM_AADA without CM_RAND_EXP
+
+   //assert(mul_bits != 0);
+
+
+   H = mult_loop(mul_bits-1, m, H, tmp, P);
+
+   if (!H->is_zero()) // cannot convert if H == O
+      {
+      *this = H->get_z_to_one();
+      }else
+      {
+      *this = *H;
+      }
+   mX.turn_off_sp_red_mul();
+   mY.turn_off_sp_red_mul();
+   mZ.turn_off_sp_red_mul();
+   return *this;
+   }
+PointGFp& PointGFp::operator*=(const BigInt& scalar)
+   {
+   // use montgomery mult. in this operation
+
+   this->turn_on_sp_red_mul();
+
+   PointGFp H(this->mC); // create as zero
+   H.turn_on_sp_red_mul();
+   PointGFp P(*this);
+   P.turn_on_sp_red_mul();
+   BigInt m(scalar);
+   if (m < BigInt(0))
+      {
+      m = -m;
+      P.negate();
+      }
+   if (P.is_zero() || (m == BigInt(0)))
+      {
+      *this = H;
+      return *this;
+      }
+   if (m == BigInt(1))
+      {
+      //*this == P already
+      return *this;
+      }
+
+   const int l = m.bits() - 1;
+   for (int i=l; i >=0; i--)
+      {
+
+      H.mult2_in_place();
+      if (m.get_bit(i))
+         {
+         H += P;
+         }
+      }
+
+   if (!H.is_zero()) // cannot convert if H == O
+      {
+      *this = H.get_z_to_one();
+      }else
+      {
+      *this = H;
+      }
+   return *this;
+   }
+
+inline std::tr1::shared_ptr<PointGFp> PointGFp::mult_loop(
+   int l,
+   const BigInt& m,
+   std::tr1::shared_ptr<PointGFp> H,
+   std::tr1::shared_ptr<PointGFp> tmp,
+   PointGFp const& P)
+   {
+
+   //assert(l >= (int)m.bits()- 1);
+   tmp = H;
+   std::tr1::shared_ptr<PointGFp> to_add(new PointGFp(P)); // we just need some point
+   // so that we can use op=
+   // inside the loop
+   for (int i=l; i >=0; i--)
+      {
+      H->mult2_in_place();
+
+#ifndef CM_AADA
+
+      if (m.get_bit(i))
+         {
+         *H += P;
+         }
+#else // (CM_AADA is in)
+
+      if (H.get() == to_add.get())
+         {
+         to_add = tmp; // otherwise all pointers might point to the same object
+         // and we always need two objects to be able to switch around
+         }
+      to_add->assign_within_same_curve(*H);
+      tmp = H;
+      *tmp += P; // tmp already points to H
+
+      if (m.get_bit(i))
+         {
+         H = tmp; // NOTE: assign the pointer, not the value!
+         // (so that the operation is fast and thus as difficult
+         // to detect as possible)
+         }
+      else
+         {
+         H = to_add; // NOTE: this is necessary, because the assignment
+         // "*tmp = ..." already changed what H pointed to
+
+
+         }
+#endif // CM_AADA
+
+      }
+   return H;
+   }
+PointGFp& PointGFp::negate()
+   {
+   if (!is_zero())
+      {
+      mY.negate();
+      }
+   return *this;
+   }
+// *this *= 2
+PointGFp& PointGFp::mult2_in_place()
+   {
+   if (is_zero())
+      {
+      return *this;
+      }
+   if (mY.is_zero())
+      {
+
+      *this = PointGFp(mC); // setting myself to zero
+      return *this;
+      }
+   ensure_worksp();
+
+   (*mp_worksp_gfp_el)[0].share_assign(mY);
+   (*mp_worksp_gfp_el)[0] *= mY;
+
+   //GFpElement S(mX * z);
+   (*mp_worksp_gfp_el)[1].share_assign(mX);
+   (*mp_worksp_gfp_el)[1] *= (*mp_worksp_gfp_el)[0];
+
+   //GFpElement x(S + S);
+   (*mp_worksp_gfp_el)[2].share_assign((*mp_worksp_gfp_el)[1]);
+   (*mp_worksp_gfp_el)[2] += (*mp_worksp_gfp_el)[1];
+
+   //S = x + x;
+   (*mp_worksp_gfp_el)[1].share_assign((*mp_worksp_gfp_el)[2]);
+   (*mp_worksp_gfp_el)[1] += (*mp_worksp_gfp_el)[2];
+
+   if (!mAZpow4_set)
+      {
+      if (mZ == *(mC.get_mres_one()))
+         {
+         mAZpow4 = mC.get_mres_a();
+         mAZpow4_set = true;
+         }
+      else
+         {
+         if (!mZpow2_set)
+            {
+            mZpow2 = mZ;
+            mZpow2 *= mZ;
+
+            mZpow2_set = true;
+            }
+         //x = mZpow2 * mZpow2;
+         (*mp_worksp_gfp_el)[2].share_assign(mZpow2);
+         (*mp_worksp_gfp_el)[2] *= mZpow2;
+
+         //mAZpow4 = mC.get_mres_a() * x;
+         mAZpow4 = mC.get_mres_a();
+         mAZpow4 *= (*mp_worksp_gfp_el)[2];
+
+         }
+
+      }
+
+   //GFpElement y(mX * mX);
+   (*mp_worksp_gfp_el)[3].share_assign(mX);
+   (*mp_worksp_gfp_el)[3] *= mX;
+
+   //GFpElement M(y + y + y + mAZpow4);
+   (*mp_worksp_gfp_el)[4].share_assign((*mp_worksp_gfp_el)[3]);
+   (*mp_worksp_gfp_el)[4] += (*mp_worksp_gfp_el)[3];
+   (*mp_worksp_gfp_el)[4] += (*mp_worksp_gfp_el)[3];
+   (*mp_worksp_gfp_el)[4] += mAZpow4;
+
+   //x = M * M - (S+S);
+   (*mp_worksp_gfp_el)[2].share_assign((*mp_worksp_gfp_el)[4]);
+   (*mp_worksp_gfp_el)[2] *= (*mp_worksp_gfp_el)[4];
+   (*mp_worksp_gfp_el)[2] -= (*mp_worksp_gfp_el)[1];
+   (*mp_worksp_gfp_el)[2] -= (*mp_worksp_gfp_el)[1];
+
+   //y = z * z;
+   (*mp_worksp_gfp_el)[3].share_assign((*mp_worksp_gfp_el)[0]);
+   (*mp_worksp_gfp_el)[3] *= (*mp_worksp_gfp_el)[0];
+
+   //GFpElement U(y + y);
+   (*mp_worksp_gfp_el)[5].share_assign((*mp_worksp_gfp_el)[3]);
+   (*mp_worksp_gfp_el)[5] += (*mp_worksp_gfp_el)[3];
+
+   //z = U + U;
+   (*mp_worksp_gfp_el)[0].share_assign((*mp_worksp_gfp_el)[5]);
+   (*mp_worksp_gfp_el)[0] += (*mp_worksp_gfp_el)[5];
+
+   //U = z + z;
+   (*mp_worksp_gfp_el)[5].share_assign((*mp_worksp_gfp_el)[0]);
+   (*mp_worksp_gfp_el)[5] += (*mp_worksp_gfp_el)[0];
+
+   //y = M * (S - x) - U;
+   (*mp_worksp_gfp_el)[3].share_assign((*mp_worksp_gfp_el)[1]);
+   (*mp_worksp_gfp_el)[3] -= (*mp_worksp_gfp_el)[2];
+   (*mp_worksp_gfp_el)[3] *= (*mp_worksp_gfp_el)[4];
+   (*mp_worksp_gfp_el)[3] -= (*mp_worksp_gfp_el)[5];
+
+   if (mZ != *(mC.get_mres_one()))
+      {
+      //z = mY * mZ;
+      (*mp_worksp_gfp_el)[0].share_assign(mY);
+      (*mp_worksp_gfp_el)[0] *= mZ;
+
+      }
+   else
+      {
+      //z = mY;
+      (*mp_worksp_gfp_el)[0].share_assign(mY);
+
+      }
+   //z = z + z;
+   (*mp_worksp_gfp_el)[6].share_assign((*mp_worksp_gfp_el)[0]);
+   (*mp_worksp_gfp_el)[0] += (*mp_worksp_gfp_el)[6];
+
+   //mX = x;
+   //mY = y;
+   //mZ = z;
+   mX = (*mp_worksp_gfp_el)[2];
+   mY = (*mp_worksp_gfp_el)[3];
+   mZ = (*mp_worksp_gfp_el)[0];
+
+   mZpow2_set = false;
+   mZpow3_set = false;
+   mAZpow4_set = false;
+   return *this;
+
+   }
+void PointGFp::turn_on_sp_red_mul() const
+   {
+   mX.turn_on_sp_red_mul();
+   mY.turn_on_sp_red_mul();
+   mZ.turn_on_sp_red_mul();
+
+   // also pretransform, otherwise
+   // we might have bad results with respect to
+   // performance because
+   // additions/subtractions in mult2_in_place()
+   // and op+= spread untransformed GFpElements
+   mX.get_mres();
+   mY.get_mres();
+   mZ.get_mres();
+
+   mZpow2.turn_on_sp_red_mul();
+   mZpow3.turn_on_sp_red_mul();
+   mAZpow4.turn_on_sp_red_mul();
+   }
+// getters
+
+/**
+* returns a point equivalent to *this but were
+* Z has value one, i.e. x and y correspond to
+* their values in affine coordinates
+*/
+PointGFp const PointGFp::get_z_to_one() const
+   {
+   return PointGFp(*this).set_z_to_one();
+   }
+/**
+* changes the representation of *this so that
+* Z has value one, i.e. x and y correspond to
+* their values in affine coordinates.
+* returns *this.
+*/
+PointGFp const& PointGFp::set_z_to_one() const
+   {
+   if (!(mZ.get_value() == BigInt(1)) && !(mZ.get_value() == BigInt(0)))
+      {
+      GFpElement z = inverse(mZ);
+      GFpElement z2 = z * z;
+      z *= z2;
+      GFpElement x = mX * z2;
+      GFpElement y = mY * z;
+      mZ = GFpElement(mC.get_p(), BigInt(1));
+      mX = x;
+      mY = y;
+      }
+   else
+      {
+      if (mZ.get_value() == BigInt(0))
+         {
+         throw Illegal_Transformation("cannot convert Z to one");
+         }
+      }
+   return *this; // mZ = 1 already
+   }
+CurveGFp const PointGFp::get_curve() const
+   {
+   return mC;
+   }
+GFpElement const PointGFp::get_affine_x() const
+   {
+
+   if (is_zero())
+      {
+      throw Illegal_Transformation("cannot convert to affine");
+
+      }
+   /*if(!mZpow2_set)
+   {*/
+   mZpow2 = mZ * mZ;
+   mZpow2_set = true;
+   //}
+   //assert(mZpow2 == mZ*mZ);
+   GFpElement z2 = mZpow2;
+   return mX * z2.inverse_in_place();
+   }
+
+GFpElement const PointGFp::get_affine_y() const
+   {
+
+   if (is_zero())
+      {
+      throw Illegal_Transformation("cannot convert to affine");
+
+      }
+   /*if(!mZpow3_set )
+   {*/
+   mZpow3 = mZ * mZ * mZ;
+   mZpow3_set = true;
+   //}
+   //assert(mZpow3 == mZ * mZ *mZ);
+   GFpElement z3 = mZpow3;
+   return mY * z3.inverse_in_place();
+   }
+GFpElement const PointGFp::get_jac_proj_x() const
+   {
+   return GFpElement(mX);
+   }
+GFpElement const PointGFp::get_jac_proj_y() const
+   {
+   return GFpElement(mY);
+   }
+GFpElement const PointGFp::get_jac_proj_z() const
+   {
+   return GFpElement(mZ);
+   }
+
+// Is this the point at infinity?
+bool PointGFp::is_zero() const
+   {
+   return(mX.is_zero() && mZ.is_zero());
+   //NOTE: the calls to GFpElement::is_zero() instead of getting the value and
+   // and comparing it are import because they do not provoke backtransformations
+   // to the ordinary residue.
+   }
+
+// Is the point still on the curve??
+// (If everything is correct, the point is always on its curve; then the
+// function will return silently. If Oskar managed to corrupt this object's state,
+// then it will throw an exception.)
+void PointGFp::check_invariants() const
+   {
+   if (is_zero())
+      {
+      return;
+      }
+   const GFpElement y2 = mY * mY;
+   const GFpElement x3 = mX * mX * mX;
+
+   if (mZ.get_value() == BigInt(1))
+      {
+      GFpElement ax = mC.get_a() * mX;
+      if (y2 != (x3 + ax + mC.get_b()))
+         {
+         throw Illegal_Point();
+         }
+
+      }
+   /*if (!mZpow2_set)
+   {*/
+   mZpow2 = mZ * mZ;
+   mZpow2_set = true;
+   /*}
+   if (!mZpow3_set)
+   {*/
+   mZpow3 = mZpow2 * mZ;
+   mZpow3_set = true;
+   /*}
+   if(!mAZpow4_set)
+   {*/
+   mAZpow4 = mZpow3 * mZ * mC.get_a();
+   mAZpow4_set = true;
+   //}
+   const GFpElement aXZ4 = mAZpow4 * mX;
+   const GFpElement bZ6 = mC.get_b() * mZpow3 * mZpow3;
+
+   if (y2 != (x3 + aXZ4 + bZ6))
+      {
+      throw Illegal_Point();
+      }
+
+   }
+// swaps the states of *this and other, does not throw!
+void PointGFp::swap(PointGFp& other)
+   {
+   mC.swap(other.mC);
+   mX.swap(other.mX);
+   mY.swap(other.mY);
+   mZ.swap(other.mZ);
+   mZpow2.swap(other.mZpow2);
+   mZpow3.swap(other.mZpow3);
+   mAZpow4.swap(other.mAZpow4);
+   std::swap<bool>(mZpow2_set, other.mZpow2_set);
+   std::swap<bool>(mZpow3_set, other.mZpow3_set);
+   std::swap<bool>(mAZpow4_set, other.mAZpow4_set);
+   }
+
+PointGFp const mult2(PointGFp const& point)
+   {
+   return (PointGFp(point)).mult2_in_place();
+   }
+
+
+bool operator==(PointGFp const& lhs, PointGFp const& rhs)
+   {
+   if (lhs.is_zero() && rhs.is_zero())
+      {
+      return true;
+      }
+   if ((lhs.is_zero() && !rhs.is_zero()) || (!lhs.is_zero() && rhs.is_zero()))
+      {
+      return false;
+      }
+   // neither operand is zero, so we can call get_z_to_one()
+   //assert(!lhs.is_zero());
+   //assert(!rhs.is_zero());
+   PointGFp aff_lhs = lhs.get_z_to_one();
+   PointGFp aff_rhs = rhs.get_z_to_one();
+   return (aff_lhs.get_curve() == aff_rhs.get_curve() &&
+           aff_lhs.get_jac_proj_x() == aff_rhs.get_jac_proj_x() &&
+           aff_lhs.get_jac_proj_y() == aff_rhs.get_jac_proj_y());
+   }
+
+// arithmetic operators
+PointGFp operator+(PointGFp const& lhs, PointGFp const& rhs)
+   {
+   PointGFp tmp(lhs);
+   return tmp += rhs;
+   }
+
+PointGFp operator-(PointGFp const& lhs, PointGFp const& rhs)
+   {
+   PointGFp tmp(lhs);
+   return tmp -= rhs;
+   }
+
+PointGFp operator-(PointGFp const& lhs)
+   {
+   return PointGFp(lhs).negate();
+   }
+
+PointGFp operator*(const BigInt& scalar, PointGFp const& point)
+   {
+   PointGFp result(point);
+   return result *= scalar;
+   }
+
+PointGFp operator*(PointGFp const& point, const BigInt& scalar)
+   {
+   PointGFp result(point);
+   return result *= scalar;
+   }
+
+PointGFp mult_point_secure(PointGFp const& point, const BigInt& scalar, const BigInt& point_order, const BigInt& max_secret)
+   {
+   PointGFp result(point);
+   result.mult_this_secure(scalar, point_order, max_secret);
+   return result;
+   }
+
+// encoding and decoding
+SecureVector<byte> EC2OSP(PointGFp const& point, byte format)
+   {
+   SecureVector<byte> result;
+   if (format == PointGFp::UNCOMPRESSED)
+      {
+      result = encode_uncompressed(point);
+      }
+   else if (format == PointGFp::COMPRESSED)
+      {
+      result = encode_compressed(point);
+
+      }
+   else if (format == PointGFp::HYBRID)
+      {
+      result = encode_hybrid(point);
+      }
+   else
+      {
+      throw Format_Error("illegal point encoding format specification");
+      }
+   return result;
+   }
+SecureVector<byte> encode_compressed(PointGFp const& point)
+   {
+
+
+   if (point.is_zero())
+      {
+      SecureVector<byte> result (1);
+      result[0] = 0;
+      return result;
+
+      }
+   u32bit l = point.get_curve().get_p().bits();
+   int dummy = l & 7;
+   if (dummy != 0)
+      {
+      l += 8 - dummy;
+      }
+   l /= 8;
+   SecureVector<byte> result (l+1);
+   result[0] = 2;
+   BigInt x = point.get_affine_x().get_value();
+   SecureVector<byte> bX = BigInt::encode_1363(x, l);
+   result.copy(1, bX.begin(), bX.size());
+   BigInt y = point.get_affine_y().get_value();
+   if (y.get_bit(0))
+      {
+      result[0] |= 1;
+      }
+   return result;
+   }
+
+
+SecureVector<byte> encode_uncompressed(PointGFp const& point)
+   {
+   if (point.is_zero())
+      {
+      SecureVector<byte> result (1);
+      result[0] = 0;
+      return result;
+      }
+   u32bit l = point.get_curve().get_p().bits();
+   int dummy = l & 7;
+   if (dummy != 0)
+      {
+      l += 8 - dummy;
+      }
+   l /= 8;
+   SecureVector<byte> result (2*l+1);
+   result[0] = 4;
+   BigInt x = point.get_affine_x().get_value();
+   BigInt y = point.get_affine_y().get_value();
+   SecureVector<byte> bX = BigInt::encode_1363(x, l);
+   SecureVector<byte> bY = BigInt::encode_1363(y, l);
+   result.copy(1, bX.begin(), l);
+   result.copy(l+1, bY.begin(), l);
+   return result;
+
+   }
+
+SecureVector<byte> encode_hybrid(PointGFp const& point)
+   {
+   if (point.is_zero())
+      {
+      SecureVector<byte> result (1);
+      result[0] = 0;
+      return result;
+      }
+   u32bit l = point.get_curve().get_p().bits();
+   int dummy = l & 7;
+   if (dummy != 0)
+      {
+      l += 8 - dummy;
+      }
+   l /= 8;
+   SecureVector<byte> result (2*l+1);
+   result[0] = 6;
+   BigInt x = point.get_affine_x().get_value();
+   BigInt y = point.get_affine_y().get_value();
+   SecureVector<byte> bX = BigInt::encode_1363(x, l);
+   SecureVector<byte> bY = BigInt::encode_1363(y, l);
+   result.copy(1, bX.begin(), bX.size());
+   result.copy(l+1, bY.begin(), bY.size());
+   if (y.get_bit(0))
+      {
+      result[0] |= 1;
+      }
+   return result;
+   }
+
+PointGFp OS2ECP(MemoryRegion<byte> const& os, CurveGFp const& curve)
+   {
+   if (os.size() == 1 && os[0] == 0)
+      {
+      return PointGFp(curve); // return zero
+      }
+   SecureVector<byte> bX;
+   SecureVector<byte> bY;
+
+   GFpElement x(1,0);
+   GFpElement y(1,0);
+   GFpElement z(1,0);
+
+   const byte pc = os[0];
+   BigInt bi_dec_x;
+   BigInt bi_dec_y;
+   switch (pc)
+      {
+      case 2:
+      case 3:
+         //compressed form
+         bX = SecureVector<byte>(os.size() - 1);
+         bX.copy(os.begin()+1, os.size()-1);
+
+         /* Problem wäre, wenn decode() das erste bit als Vorzeichen interpretiert.
+         *---------------------
+         * AW(FS): decode() interpretiert das erste Bit nicht als Vorzeichen
+         */
+         bi_dec_x = BigInt::decode(bX, bX.size());
+         x = GFpElement(curve.get_p(), bi_dec_x);
+         bool yMod2;
+         yMod2 = (pc & 1) == 1;
+         y = PointGFp::decompress(yMod2, x, curve);
+         break;
+      case 4:
+         // uncompressed form
+         int l;
+         l = (os.size() -1)/2;
+         bX = SecureVector<byte>(l);
+         bY = SecureVector<byte>(l);
+         bX.copy(os.begin()+1, l);
+         bY.copy(os.begin()+1+l, l);
+         bi_dec_x = BigInt::decode(bX.begin(), bX.size());
+
+         bi_dec_y = BigInt::decode(bY.begin(),bY.size());
+         x = GFpElement(curve.get_p(), bi_dec_x);
+         y = GFpElement(curve.get_p(), bi_dec_y);
+         break;
+
+      case 6:
+      case 7:
+         //hybrid form
+         l = (os.size() - 1)/2;
+         bX = SecureVector<byte>(l);
+         bY = SecureVector<byte>(l);
+         bX.copy(os.begin() + 1, l);
+         bY.copy(os.begin()+1+l, l);
+         yMod2 = (pc & 0x01) == 1;
+         if (!(PointGFp::decompress(yMod2, x, curve) == y))
+            {
+            throw Illegal_Point("error during decoding hybrid format");
+            }
+         break;
+      default:
+         throw Format_Error("encountered illegal format specification while decoding point");
+      }
+   z = GFpElement(curve.get_p(), BigInt(1));
+   //assert((x.is_trf_to_mres() && x.is_use_montgm()) || !x.is_trf_to_mres());
+   //assert((y.is_trf_to_mres() && y.is_use_montgm()) || !y.is_trf_to_mres());
+   //assert((z.is_trf_to_mres() && z.is_use_montgm()) || !z.is_trf_to_mres());
+   PointGFp result(curve, x, y, z);
+   result.check_invariants();
+   //assert((result.get_jac_proj_x().is_trf_to_mres() && result.get_jac_proj_x().is_use_montgm()) || !result.get_jac_proj_x().is_trf_to_mres());
+   //assert((result.get_jac_proj_y().is_trf_to_mres() && result.get_jac_proj_y().is_use_montgm()) || !result.get_jac_proj_y().is_trf_to_mres());
+   //assert((result.get_jac_proj_z().is_trf_to_mres() && result.get_jac_proj_z().is_use_montgm()) || !result.get_jac_proj_z().is_trf_to_mres());
+   return result;
+   }
+
+GFpElement PointGFp::decompress(bool yMod2, GFpElement const& x,
+                                CurveGFp const& curve)
+   {
+   BigInt xVal = x.get_value();
+   BigInt xpow3 = xVal * xVal * xVal;
+   BigInt g = curve.get_a().get_value() * xVal;
+   g += xpow3;
+   g += curve.get_b().get_value();
+   g = g%curve.get_p();
+   BigInt z = ressol(g, curve.get_p());
+
+   if(z < 0)
+      throw Illegal_Point("error during decompression");
+
+   bool zMod2 = z.get_bit(0);
+   if ((zMod2 && ! yMod2) || (!zMod2 && yMod2))
+      {
+      z = curve.get_p() - z;
+      }
+   return GFpElement(curve.get_p(),z);
+   }
+
+PointGFp const create_random_point(RandomNumberGenerator& rng,
+                                   CurveGFp const& curve)
+   {
+
+   // create a random point
+   GFpElement mX(1,1);
+   GFpElement mY(1,1);
+   GFpElement mZ(1,1);
+   GFpElement minusOne(curve.get_p(), BigInt(BigInt::Negative,1));
+   mY = minusOne;
+   GFpElement y2(1,1);
+   GFpElement x(1,1);
+
+   while (mY == minusOne)
+      {
+      BigInt value(rng, curve.get_p().bits());
+      mX = GFpElement(curve.get_p(),value);
+      y2 = curve.get_a() * mX;
+      x = mX * mX;
+      x *= mX;
+      y2 += (x + curve.get_b());
+
+      value = ressol(y2.get_value(), curve.get_p());
+
+      if(value < 0)
+         mY = minusOne;
+      else
+         mY = GFpElement(curve.get_p(), value);
+      }
+   mZ = GFpElement(curve.get_p(), BigInt(1));
+
+   return PointGFp(curve, mX, mY, mZ);
+   }
+
+} // namespace Botan
diff --git a/src/math/gfpmath/point_gfp.h b/src/math/gfpmath/point_gfp.h
new file mode 100644
index 000000000..7e5aec379
--- /dev/null
+++ b/src/math/gfpmath/point_gfp.h
@@ -0,0 +1,307 @@
+/*************************************************
+* Arithmetic over GF(p)                          *
+*                                                *
+* (C) 2007 Martin Doering                        *
+*          Christoph Ludwig                      *
+*          Falko Strenzke                        *
+* (C) 2008 Jack Lloyd                            *
+*************************************************/
+
+#ifndef BOTAN_POINT_GFP_H__
+#define BOTAN_POINT_GFP_H__
+
+#include <botan/curve_gfp.h>
+#include <botan/gfp_element.h>
+#include <botan/bigint.h>
+#include <botan/exceptn.h>
+#include <vector>
+
+namespace Botan {
+
+struct Illegal_Point : public Exception
+   {
+   Illegal_Point(const std::string& err = "") : Exception(err) {}
+   };
+
+/**
+* This class represents one point on a curve of GF(p).
+*/
+class PointGFp
+   {
+   public:
+      /**
+      * uncompressed encoding byte value
+      */
+      static const int UNCOMPRESSED = 0;
+
+      /**
+      * compressed encoding byte value
+      */
+      static const int COMPRESSED = 1;
+
+      /**
+      * hybrid encoding byte value
+      */
+      static const int HYBRID = 2;
+
+      /**
+      * Construct the point O
+      * @param curve The base curve
+      */
+      explicit PointGFp(CurveGFp const& curve);
+
+      /**
+      * Construct a point given its affine coordinates
+      * @param curve the base curve
+      * @param x affine x coordinate
+      * @param y affine y coordinate
+      */
+      explicit PointGFp(CurveGFp const& curve, GFpElement const& x,
+                          GFpElement const& y );
+
+      /**
+      * Construct a point given its jacobian projective coordinates
+      * @param curve the base curve
+      * @param x jacobian projective x coordinate
+      * @param y jacobian projective y coordinate
+      * @param z jacobian projective y coordinate
+      */
+      explicit PointGFp(CurveGFp const& curve, GFpElement const& x,
+                          GFpElement const& y, GFpElement const& z );
+
+      /**
+      * copy constructor
+      * @param other the value to clone
+      */
+      PointGFp(PointGFp const& other );
+
+      /**
+      * assignment operator
+      * @param other The point to use as source for the assignment
+      */
+      PointGFp const& operator=(PointGFp const& other );
+
+      /**
+      * assign another point which is on the same curve as *this
+      * @param other The point to use as source for the assignment
+      */
+      PointGFp const& assign_within_same_curve(PointGFp const& other);
+
+
+
+      /**
+      * += Operator
+      * @param rhs the PointGFp to add to the local value
+      * @result resulting PointGFp
+      */
+      PointGFp& operator+=(PointGFp const& rhs );
+
+      /**
+      * -= Operator
+      * @param rhs the PointGFp to subtract from the local value
+      * @result resulting PointGFp
+      */
+      PointGFp& operator-=(PointGFp const& rhs );
+
+      /**
+      * *= Operator
+      * This function turns on the the special reduction multiplication
+      * itself for fast computation, turns it off again when finished.
+      * @param scalar the PointGFp to multiply with *this
+      * @result resulting PointGFp
+      */
+      PointGFp& operator*=(const BigInt& scalar );
+
+      /**
+      * the equivalent to operator*= with countermeasures against
+      * sidechannel attacks, using the randomized exponent
+      * and add-and-double-always
+      * countermeasures (suitable for ECDSA and ECKAEG)
+      * @param scalar the scalar to multiply the point with
+      * @param point_order a multiple of the order of the point
+      *(= n * k in the general case; k is the cofactor)
+      * @param max_secr the maximal size of the scalar
+      * (will usually be  n-1 )
+      * @result resulting PointGFp
+      */
+      PointGFp& mult_this_secure(const BigInt& scalar,
+                                 const BigInt& point_order,
+                                 const BigInt& max_secr
+         );
+
+      /**
+      * Negate internal value(*this *= -1 )
+      * @return *this
+      */
+      PointGFp& negate();
+
+      /**
+      * Multiply the point by two(*this *= 2 )
+      * @return *this
+      */
+      PointGFp& mult2_in_place();
+
+      /**
+      * Set z coordinate to one.
+      * @return *this
+      */
+      PointGFp const& set_z_to_one() const;
+
+      /**
+      * Turn on the special reduction multiplication (i.e. the
+      * Montgomery multiplication in the current implementation) for
+      * the coordinates. This enables fast execution of mult2_in_place()
+      * and operator+=().
+      */
+      void turn_on_sp_red_mul() const;
+
+      /**
+      * Return a point
+      * where the coordinates are transformed
+      * so that z equals one,
+      * thus x and y have just the affine values.
+      * @result *this
+      */
+      PointGFp const get_z_to_one() const;
+
+      /**
+      * Return base curve of this point
+      * @result the curve over GF(p) of this point
+      */
+      CurveGFp const get_curve() const;
+
+      /**
+      * get affine x coordinate
+      * @result affine x coordinate
+      */
+      GFpElement const get_affine_x() const;
+
+      /**
+      * get affine y coordinate
+      * @result affine y coordinate
+      */
+      GFpElement const get_affine_y() const;
+
+      /**
+      * get the jacobian projective x coordinate
+      * @result jacobian projective x coordinate
+      */
+      GFpElement const get_jac_proj_x() const;
+
+      /**
+      * get the jacobian projective y coordinate
+      * @result jacobian projective y coordinate
+      */
+      GFpElement const get_jac_proj_y() const;
+
+      /**
+      * get the jacobian projective z coordinate
+      * @result jacobian projective z coordinate
+      */
+      GFpElement const get_jac_proj_z() const;
+
+      /**
+      * Is this the point at infinity?
+      * @result true, if this point is at infinity, false otherwise.
+      */
+      bool is_zero() const;
+
+      /**
+      *  Checks whether the point is to be found on the underlying curve.
+      *  Throws an Invalid_Point exception in case of detecting that the point
+      *  does not satisfy the curve equation.
+      *  To be used to ensure against fault attacks.
+      */
+      void check_invariants() const;
+
+
+      /**
+      *  swaps the states of *this and other, does not throw!
+      * @param other the object to swap values with
+      */
+      void swap(PointGFp& other );
+
+      /**
+      * Sets the shared pointer to the GFpModulus that will be
+      * held in *this, specifically the various members of *this.
+      * Warning: do not use this function unless you know in detail about
+      * the implications of using
+      * the shared GFpModulus objects!
+      * Do NOT spread a shared pointer to GFpModulus over different
+      * threads!
+      * @param mod a shared pointer to a GFpModulus that will
+      * be held in the members *this
+      */
+      void set_shrd_mod(std::tr1::shared_ptr<Botan::GFpModulus> p_mod);
+
+      static GFpElement decompress(bool yMod2, GFpElement const& x, CurveGFp const& curve );
+
+   private:
+      static const u32bit GFPEL_WKSP_SIZE = 9;
+      void ensure_worksp() const;
+
+      inline std::tr1::shared_ptr<PointGFp> mult_loop(int l, const BigInt& m, std::tr1::shared_ptr<PointGFp> H, std::tr1::shared_ptr<PointGFp> tmp, PointGFp const& P);
+
+      CurveGFp mC;
+      mutable GFpElement mX;  // NOTE: these values must be mutable (affine<->proj)
+      mutable GFpElement mY;
+      mutable GFpElement mZ;
+      mutable GFpElement mZpow2;  // mZ^2
+      mutable GFpElement mZpow3;   // mZ^3
+      mutable GFpElement mAZpow4;  // mA*mZ^4
+      mutable bool mZpow2_set;
+      mutable bool mZpow3_set;
+      mutable bool mAZpow4_set;
+      mutable std::tr1::shared_ptr<std::vector<GFpElement> > mp_worksp_gfp_el;
+
+   };
+
+// relational operators
+bool operator==(PointGFp const& lhs, PointGFp const& rhs );
+inline bool operator!=(PointGFp const& lhs, PointGFp const& rhs )
+   {
+   return !operator==(lhs, rhs );
+   }
+
+// arithmetic operators
+PointGFp operator+(PointGFp const& lhs, PointGFp const& rhs );
+PointGFp operator-(PointGFp const& lhs, PointGFp const& rhs );
+PointGFp operator-(PointGFp const& lhs );
+
+PointGFp operator*(const BigInt& scalar, PointGFp const& point );
+PointGFp operator*(PointGFp const& point, const BigInt& scalar );
+PointGFp mult_point_secure(PointGFp const& point, const BigInt& scalar, const BigInt& point_order, const BigInt& max_secret);
+
+PointGFp const mult2 (PointGFp const& point);
+
+PointGFp const create_random_point(RandomNumberGenerator& rng,
+                                   CurveGFp const& curve);
+
+// encoding and decoding
+SecureVector<byte> EC2OSP(PointGFp const& point, byte format );
+PointGFp OS2ECP(MemoryRegion<byte> const& os, CurveGFp const& curve );
+
+SecureVector<byte> encode_uncompressed(PointGFp const& point ); // maybe make private
+SecureVector<byte> encode_hybrid(PointGFp const& point ); // maybe make private
+SecureVector<byte> encode_compressed(PointGFp const& point ); // maybe make private
+
+// swaps the states of point1 and point2, does not throw!
+// cf. Meyers, Item 25
+inline
+void swap(PointGFp& point1, PointGFp& point2 )
+   {
+   point1.swap(point2 );
+   }
+
+} // namespace Botan
+
+namespace std {
+
+// swaps the states of point1 and point2, does not throw!
+// cf. Meyers, Item 25
+template<> inline void
+swap<Botan::PointGFp>(Botan::PointGFp& x, Botan::PointGFp& y) { x.swap(y); }
+
+} // namespace std
+
+#endif