diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..05a9715
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,48 @@
+FLAGS_COMMON:=-Wall
+FLAGS_PROD:=-DNDEBUG -g0 -O3 -march=native
+FLAGS_DEBUG:=-DVERIFY_MAGNITUDE -ggdb3 -O1
+FLAGS_TEST:=-DVERIFY_MAGNITUDE -ggdb3 -O3 -march=native
+
+SECP256K1_FILES := num.h   field.h   group.h   ecmult.h   ecdsa.h \
+                   num.cpp field.cpp group.cpp ecmult.cpp ecdsa.cpp \
+
+ifndef CONF
+CONF := gmp
+endif
+
+ifeq ($(CONF), openssl)
+FLAGS_CONF:=-DUSE_NUM_OPENSSL -DUSE_FIELDINVERSE_BUILTIN
+LIBS:=-lcrypto
+SECP256K1_FILES := $(SECP256K1_FILES) num_openssl.h num_openssl.cpp
+else
+ifeq ($(CONF), gmp)
+FLAGS_CONF:=-DUSE_NUM_GMP
+LIBS:=-lgmp
+SECP256K1_FILES := $(SECP256K1_FILES) num_gmp.h num_gmp.cpp 
+endif
+endif
+
+all: *.cpp *.h
+	+make CONF=openssl all-openssl
+	+make CONF=gmp     all-gmp
+
+clean:
+	+make CONF=openssl clean-openssl
+	+make CONF=gmp     clean-gmp
+
+bench-any: bench-$(CONF)
+tests-any: tests-$(CONF)
+
+all-$(CONF): bench-$(CONF) tests-$(CONF)
+
+clean-$(CONF):
+	rm -f bench-$(CONF) tests-$(CONF) obj/secp256k1-$(CONF).o
+
+obj/secp256k1-$(CONF).o: $(SECP256K1_FILES)
+	$(CXX) $(FLAGS_COMMON) $(FLAGS_PROD) $(FLAGS_CONF) secp256k1.cpp -c -o obj/secp256k1-$(CONF).o
+
+bench-$(CONF): obj/secp256k1-$(CONF).o bench.cpp
+	$(CXX) $(FLAGS_COMMON) $(FLAGS_PROD) $(FLAGS_CONF) obj/secp256k1-$(CONF).o bench.cpp $(LIBS) -o bench-$(CONF)
+
+tests-$(CONF): $(SECP256K1_FILES) tests.cpp
+	$(CXX) $(FLAGS_COMMON) $(FLAGS_TEST) $(FLAGS_CONF) tests.cpp $(LIBS) -o tests-$(CONF)
diff --git a/bench.cpp b/bench.cpp
new file mode 100644
index 0000000..9450a75
--- /dev/null
+++ b/bench.cpp
@@ -0,0 +1,32 @@
+#include <stdio.h>
+
+#include "num.h"
+#include "field.h"
+#include "group.h"
+#include "ecmult.h"
+#include "ecdsa.h"
+
+using namespace secp256k1;
+
+int main() {
+    FieldElem x;
+    const Number &order = GetGroupConst().order;
+    Number r, s, m;
+    Signature sig;
+    x.SetHex("a357ae915c4a65281309edf20504740f0eb3343990216b4f81063cb65f2f7e0f");
+    int cnt = 0;
+    int good = 0;
+    for (int i=0; i<1000000; i++) {
+        r.SetPseudoRand(order);
+        s.SetPseudoRand(order);
+        m.SetPseudoRand(order);
+        sig.SetRS(r,s);
+        GroupElemJac pubkey; pubkey.SetCompressed(x, true);
+        if (pubkey.IsValid()) {
+            cnt++;
+            good += sig.Verify(pubkey, m);
+        }
+     }
+    printf("%i/%i\n", good, cnt);
+    return 0;
+}
diff --git a/ecdsa.cpp b/ecdsa.cpp
new file mode 100644
index 0000000..1b8ab69
--- /dev/null
+++ b/ecdsa.cpp
@@ -0,0 +1,82 @@
+#include "num.h"
+#include "field.h"
+#include "group.h"
+#include "ecmult.h"
+#include "ecdsa.h"
+
+namespace secp256k1 {
+
+bool ParsePubKey(GroupElemJac &elem, const unsigned char *pub, int size) {
+    if (size == 33 && (pub[0] == 0x02 || pub[0] == 0x03)) {
+        FieldElem x;
+        x.SetBytes(pub+1);
+        elem.SetCompressed(x, pub[0] == 0x03);
+    } else if (size == 65 && (pub[0] == 0x04 || pub[0] == 0x06 || pub[0] == 0x07)) {
+        FieldElem x,y;
+        x.SetBytes(pub+1);
+        y.SetBytes(pub+33);
+        elem = GroupElem(x,y);
+        if ((pub[0] == 0x06 || pub[0] == 0x07) && y.IsOdd() != (pub[0] == 0x07))
+            return false;
+    } else {
+        return false;
+    }
+    return elem.IsValid();
+}
+
+bool Signature::Parse(const unsigned char *sig, int size) {
+    if (sig[0] != 0x30) return false;
+    int lenr = sig[3];
+    if (5+lenr >= size) return false;
+    int lens = sig[lenr+5];
+    if (sig[1] != lenr+lens+4) return false;
+    if (lenr+lens+6 > size) return false;
+    if (sig[2] != 0x02) return false;
+    if (lenr == 0) return false;
+    if (sig[lenr+4] != 0x02) return false;
+    if (lens == 0) return false;
+    r.SetBytes(sig+4, lenr);
+    s.SetBytes(sig+6+lenr, lens);
+    return true;
+}
+
+bool Signature::RecomputeR(Number &r2, const GroupElemJac &pubkey, const Number &message) {
+    const GroupConstants &c = GetGroupConst();
+
+    if (r.IsNeg() || s.IsNeg())
+        return false;
+    if (r.IsZero() || s.IsZero())
+        return false;
+    if (r.Compare(c.order) >= 0 || s.Compare(c.order) >= 0)
+        return false;
+
+    Number sn, u1, u2;
+    sn.SetModInverse(s, c.order);
+    u1.SetModMul(sn, message, c.order);
+    u2.SetModMul(sn, r, c.order);
+    GroupElemJac pr; ECMult(pr, pubkey, u2, u1);
+    if (pr.IsInfinity())
+        return false;
+    FieldElem xr; pr.GetX(xr);
+    unsigned char xrb[32]; xr.GetBytes(xrb);
+    r2.SetBytes(xrb,32); r2.SetMod(r2,c.order);
+    return true;
+}
+
+bool Signature::Verify(const GroupElemJac &pubkey, const Number &message) {
+    Number r2;
+    if (!RecomputeR(r2, pubkey, message))
+        return false;
+    return r2.Compare(r) == 0;
+}
+
+void Signature::SetRS(const Number &rin, const Number &sin) {
+    r = rin;
+    s = sin;
+}
+
+std::string Signature::ToString() const {
+    return "(" + r.ToString() + "," + s.ToString() + ")";
+}
+
+}
diff --git a/ecdsa.h b/ecdsa.h
index 73fdab1..4763dc3 100644
--- a/ecdsa.h
+++ b/ecdsa.h
@@ -3,104 +3,18 @@
 
 namespace secp256k1 {
 
-bool ParsePubKey(GroupElemJac &elem, const unsigned char *pub, int size) {
-    if (size == 33 && (pub[0] == 0x02 || pub[0] == 0x03)) {
-        FieldElem x;
-        x.SetBytes(pub+1);
-        elem.SetCompressed(x, pub[0] == 0x03);
-    } else if (size == 65 && (pub[0] == 0x04 || pub[0] == 0x06 || pub[0] == 0x07)) {
-        FieldElem x,y;
-        x.SetBytes(pub+1);
-        y.SetBytes(pub+33);
-        elem = GroupElem(x,y);
-        if ((pub[0] == 0x06 || pub[0] == 0x07) && y.IsOdd() != (pub[0] == 0x07))
-            return false;
-    } else {
-        return false;
-    }
-    return elem.IsValid();
-}
-
 class Signature {
 private:
     Number r,s;
 
 public:
-    bool Parse(const unsigned char *sig, int size) {
-        if (sig[0] != 0x30) return false;
-        int lenr = sig[3];
-        if (5+lenr >= size) return false;
-        int lens = sig[lenr+5];
-        if (sig[1] != lenr+lens+4) return false;
-        if (lenr+lens+6 > size) return false;
-        if (sig[2] != 0x02) return false;
-        if (lenr == 0) return false;
-        if (sig[lenr+4] != 0x02) return false;
-        if (lens == 0) return false;
-        r.SetBytes(sig+4, lenr);
-        s.SetBytes(sig+6+lenr, lens);
-        return true;
-    }
-
-    bool RecomputeR(Number &r2, const GroupElemJac &pubkey, const Number &message) {
-        const GroupConstants &c = GetGroupConst();
-
-        if (r.IsNeg() || s.IsNeg())
-            return false;
-        if (r.IsZero() || s.IsZero())
-            return false;
-        if (r.Compare(c.order) >= 0 || s.Compare(c.order) >= 0)
-            return false;
-
-        Number sn, u1, u2;
-        sn.SetModInverse(s, c.order);
-        u1.SetModMul(sn, message, c.order);
-        u2.SetModMul(sn, r, c.order);
-        GroupElemJac pr; ECMult(pr, pubkey, u2, u1);
-        if (pr.IsInfinity())
-            return false;
-        FieldElem xr; pr.GetX(xr);
-        unsigned char xrb[32]; xr.GetBytes(xrb);
-        r2.SetBytes(xrb,32); r2.SetMod(r2,c.order);
-        return true;
-    }
-
-    bool Verify(const GroupElemJac &pubkey, const Number &message) {
-        Number r2;
-        if (!RecomputeR(r2, pubkey, message))
-            return false;
-        return r2.Compare(r) == 0;
-    }
-
-    void SetRS(const Number &rin, const Number &sin) {
-        r = rin;
-        s = sin;
-    }
-
-    std::string ToString() const {
-        return "(" + r.ToString() + "," + s.ToString() + ")";
-    }
+    bool Parse(const unsigned char *sig, int size);
+    bool RecomputeR(Number &r2, const GroupElemJac &pubkey, const Number &message);
+    bool Verify(const GroupElemJac &pubkey, const Number &message);
+    void SetRS(const Number &rin, const Number &sin);
+    std::string ToString() const;
 };
 
-int VerifyECDSA(const unsigned char *msg, int msglen, const unsigned char *sig, int siglen, const unsigned char *pubkey, int pubkeylen) {
-    Number m; 
-    Signature s;
-    GroupElemJac q;
-    m.SetBytes(msg, msglen);
-    if (!ParsePubKey(q, pubkey, pubkeylen))
-        return -1;
-    if (!s.Parse(sig, siglen)) {
-        fprintf(stderr, "Can't parse signature: ");
-        for (int i=0; i<siglen; i++) fprintf(stderr,"%02x", sig[i]);
-        fprintf(stderr, "\n");
-        return -2;
-    }
-//    fprintf(stderr, "Verifying ECDSA: msg=%s pubkey=%s sig=%s\n", m.ToString().c_str(), q.ToString().c_str(), s.ToString().c_str());
-    if (!s.Verify(q, m))
-        return 0;
-    return 1;
-}
-
 }
 
 #endif
diff --git a/ecmult.cpp b/ecmult.cpp
new file mode 100644
index 0000000..e48c09b
--- /dev/null
+++ b/ecmult.cpp
@@ -0,0 +1,185 @@
+#include <sstream>
+#include <algorithm>
+
+#include "num.h"
+#include "group.h"
+#include "ecmult.h"
+
+// optimal for 128-bit and 256-bit exponents
+#define WINDOW_A 5
+
+// larger numbers may result in slightly better performance, at the cost of
+// exponentially larger precomputed tables. WINDOW_G == 13 results in 640 KiB.
+#define WINDOW_G 14
+
+namespace secp256k1 {
+
+template<typename G, int W> class WNAFPrecomp {
+private:
+    G pre[1 << (W-2)];
+
+public:
+    WNAFPrecomp() {}
+
+    void Build(const G &base) {
+        pre[0] = base;
+        GroupElemJac x(base);
+        GroupElemJac d; d.SetDouble(x);
+        for (int i=1; i<(1 << (W-2)); i++) {
+            x.SetAdd(d,pre[i-1]);
+            pre[i].SetJac(x);
+        }
+    }
+
+    WNAFPrecomp(const G &base) {
+        Build(base);
+    }
+
+    void Get(G &out, int exp) const {
+        assert((exp & 1) == 1);
+        assert(exp >= -((1 << (W-1)) - 1));
+        assert(exp <=  ((1 << (W-1)) - 1));
+        if (exp > 0) {
+            out = pre[(exp-1)/2];
+        } else {
+            out.SetNeg(pre[(-exp-1)/2]);
+        }
+    }
+};
+
+template<int B> class WNAF {
+private:
+    int naf[B+1];
+    int used;
+
+    void PushNAF(int num, int zeroes) {
+        assert(used < B+1);
+        for (int i=0; i<zeroes; i++) {
+            naf[used++]=0;
+        }
+        naf[used++]=num;
+    }
+
+public:
+    WNAF(const Number &exp, int w) : used(0) {
+        int zeroes = 0;
+        Number x;
+        x.SetNumber(exp);
+        int sign = 1;
+        if (x.IsNeg()) {
+            sign = -1;
+            x.Negate();
+        }
+        while (!x.IsZero()) {
+            while (!x.IsOdd()) {
+                zeroes++;
+                x.Shift1();
+            }
+            int word = x.ShiftLowBits(w);
+            if (word & (1 << (w-1))) {
+                x.Inc();
+                PushNAF(sign * (word - (1 << w)), zeroes);
+            } else {
+                PushNAF(sign * word, zeroes);
+            }
+            zeroes = w-1;
+        }
+    }
+
+    int GetSize() const {
+        return used;
+    }
+
+    int Get(int pos) const {
+        assert(pos >= 0 && pos < used);
+        return naf[pos];
+    }
+
+    std::string ToString() {
+        std::stringstream ss;
+        ss << "(";
+        for (int i=0; i<GetSize(); i++) {
+            ss << Get(used-1-i);
+            if (i != used-1)
+                ss << ',';
+        }
+        ss << ")";
+        return ss.str();
+    }
+};
+
+class ECMultConsts {
+public:
+    WNAFPrecomp<GroupElem,WINDOW_G> wpg;
+    WNAFPrecomp<GroupElem,WINDOW_G> wpg128;
+
+    ECMultConsts() {
+        const GroupElem &g = GetGroupConst().g;
+        GroupElemJac g128j(g);
+        for (int i=0; i<128; i++)
+            g128j.SetDouble(g128j);
+        GroupElem g128; g128.SetJac(g128j);
+        wpg.Build(g);
+        wpg128.Build(g128);
+    }
+};
+
+const ECMultConsts &GetECMultConsts() {
+    static const ECMultConsts ecmult_consts;
+    return ecmult_consts;
+}
+
+void ECMult(GroupElemJac &out, const GroupElemJac &a, const Number &an, const Number &gn) {
+    Number an1, an2;
+    Number gn1, gn2;
+
+    SplitExp(an, an1, an2);
+//    printf("an=%s\n", an.ToString().c_str());
+//    printf("an1=%s\n", an1.ToString().c_str());
+//    printf("an2=%s\n", an2.ToString().c_str());
+//    printf("an1.len=%i\n", an1.GetBits());
+//    printf("an2.len=%i\n", an2.GetBits());
+    gn.SplitInto(128, gn1, gn2);
+
+    WNAF<128> wa1(an1, WINDOW_A);
+    WNAF<128> wa2(an2, WINDOW_A);
+    WNAF<128> wg1(gn1, WINDOW_G);
+    WNAF<128> wg2(gn2, WINDOW_G);
+    GroupElemJac a2; a2.SetMulLambda(a);
+    WNAFPrecomp<GroupElemJac,WINDOW_A> wpa1(a);
+    WNAFPrecomp<GroupElemJac,WINDOW_A> wpa2(a2);
+    const ECMultConsts &c = GetECMultConsts();
+
+    int size_a1 = wa1.GetSize();
+    int size_a2 = wa2.GetSize();
+    int size_g1 = wg1.GetSize();
+    int size_g2 = wg2.GetSize();
+    int size = std::max(std::max(size_a1, size_a2), std::max(size_g1, size_g2));
+
+    out = GroupElemJac();
+    GroupElemJac tmpj;
+    GroupElem tmpa;
+
+    for (int i=size-1; i>=0; i--) {
+        out.SetDouble(out);
+        int nw;
+        if (i < size_a1 && (nw = wa1.Get(i))) {
+            wpa1.Get(tmpj, nw);
+            out.SetAdd(out, tmpj);
+        }
+        if (i < size_a2 && (nw = wa2.Get(i))) {
+            wpa2.Get(tmpj, nw);
+            out.SetAdd(out, tmpj);
+        }
+        if (i < size_g1 && (nw = wg1.Get(i))) {
+            c.wpg.Get(tmpa, nw);
+            out.SetAdd(out, tmpa);
+        }
+        if (i < size_g2 && (nw = wg2.Get(i))) {
+            c.wpg128.Get(tmpa, nw);
+            out.SetAdd(out, tmpa);
+        }
+    }
+}
+
+}
diff --git a/ecmult.h b/ecmult.h
index b145bcf..0ad57e7 100644
--- a/ecmult.h
+++ b/ecmult.h
@@ -1,188 +1,12 @@
 #ifndef _SECP256K1_ECMULT_
 #define _SECP256K1_ECMULT_
 
-#include <sstream>
-#include <algorithm>
-
 #include "group.h"
 #include "num.h"
 
-// optimal for 128-bit and 256-bit exponents
-#define WINDOW_A 5
-
-// larger numbers may result in slightly better performance, at the cost of
-// exponentially larger precomputed tables. WINDOW_G == 13 results in 640 KiB.
-#define WINDOW_G 14
-
 namespace secp256k1 {
 
-template<typename G, int W> class WNAFPrecomp {
-private:
-    G pre[1 << (W-2)];
-
-public:
-    WNAFPrecomp() {}
-
-    void Build(const G &base) {
-        pre[0] = base;
-        GroupElemJac x(base);
-        GroupElemJac d; d.SetDouble(x);
-        for (int i=1; i<(1 << (W-2)); i++) {
-            x.SetAdd(d,pre[i-1]);
-            pre[i].SetJac(x);
-        }
-    }
-
-    WNAFPrecomp(const G &base) {
-        Build(base);
-    }
-
-    void Get(G &out, int exp) const {
-        assert((exp & 1) == 1);
-        assert(exp >= -((1 << (W-1)) - 1));
-        assert(exp <=  ((1 << (W-1)) - 1));
-        if (exp > 0) {
-            out = pre[(exp-1)/2];
-        } else {
-            out.SetNeg(pre[(-exp-1)/2]);
-        }
-    }
-};
-
-template<int B> class WNAF {
-private:
-    int naf[B+1];
-    int used;
-
-    void PushNAF(int num, int zeroes) {
-        assert(used < B+1);
-        for (int i=0; i<zeroes; i++) {
-            naf[used++]=0;
-        }
-        naf[used++]=num;
-    }
-
-public:
-    WNAF(const Number &exp, int w) : used(0) {
-        int zeroes = 0;
-        Number x;
-        x.SetNumber(exp);
-        int sign = 1;
-        if (x.IsNeg()) {
-            sign = -1;
-            x.Negate();
-        }
-        while (!x.IsZero()) {
-            while (!x.IsOdd()) {
-                zeroes++;
-                x.Shift1();
-            }
-            int word = x.ShiftLowBits(w);
-            if (word & (1 << (w-1))) {
-                x.Inc();
-                PushNAF(sign * (word - (1 << w)), zeroes);
-            } else {
-                PushNAF(sign * word, zeroes);
-            }
-            zeroes = w-1;
-        }
-    }
-
-    int GetSize() const {
-        return used;
-    }
-
-    int Get(int pos) const {
-        assert(pos >= 0 && pos < used);
-        return naf[pos];
-    }
-
-    std::string ToString() {
-        std::stringstream ss;
-        ss << "(";
-        for (int i=0; i<GetSize(); i++) {
-            ss << Get(used-1-i);
-            if (i != used-1)
-                ss << ',';
-        }
-        ss << ")";
-        return ss.str();
-    }
-};
-
-class ECMultConsts {
-public:
-    WNAFPrecomp<GroupElem,WINDOW_G> wpg;
-    WNAFPrecomp<GroupElem,WINDOW_G> wpg128;
-
-    ECMultConsts() {
-        const GroupElem &g = GetGroupConst().g;
-        GroupElemJac g128j(g);
-        for (int i=0; i<128; i++)
-            g128j.SetDouble(g128j);
-        GroupElem g128; g128.SetJac(g128j);
-        wpg.Build(g);
-        wpg128.Build(g128);
-    }
-};
-
-const ECMultConsts &GetECMultConsts() {
-    static const ECMultConsts ecmult_consts;
-    return ecmult_consts;
-}
-
-void ECMult(GroupElemJac &out, const GroupElemJac &a, const Number &an, const Number &gn) {
-    Number an1, an2;
-    Number gn1, gn2;
-
-    SplitExp(an, an1, an2);
-//    printf("an=%s\n", an.ToString().c_str());
-//    printf("an1=%s\n", an1.ToString().c_str());
-//    printf("an2=%s\n", an2.ToString().c_str());
-//    printf("an1.len=%i\n", an1.GetBits());
-//    printf("an2.len=%i\n", an2.GetBits());
-    gn.SplitInto(128, gn1, gn2);
-
-    WNAF<128> wa1(an1, WINDOW_A);
-    WNAF<128> wa2(an2, WINDOW_A);
-    WNAF<128> wg1(gn1, WINDOW_G);
-    WNAF<128> wg2(gn2, WINDOW_G);
-    GroupElemJac a2; a2.SetMulLambda(a);
-    WNAFPrecomp<GroupElemJac,WINDOW_A> wpa1(a);
-    WNAFPrecomp<GroupElemJac,WINDOW_A> wpa2(a2);
-    const ECMultConsts &c = GetECMultConsts();
-
-    int size_a1 = wa1.GetSize();
-    int size_a2 = wa2.GetSize();
-    int size_g1 = wg1.GetSize();
-    int size_g2 = wg2.GetSize();
-    int size = std::max(std::max(size_a1, size_a2), std::max(size_g1, size_g2));
-
-    out = GroupElemJac();
-    GroupElemJac tmpj;
-    GroupElem tmpa;
-
-    for (int i=size-1; i>=0; i--) {
-        out.SetDouble(out);
-        int nw;
-        if (i < size_a1 && (nw = wa1.Get(i))) {
-            wpa1.Get(tmpj, nw);
-            out.SetAdd(out, tmpj);
-        }
-        if (i < size_a2 && (nw = wa2.Get(i))) {
-            wpa2.Get(tmpj, nw);
-            out.SetAdd(out, tmpj);
-        }
-        if (i < size_g1 && (nw = wg1.Get(i))) {
-            c.wpg.Get(tmpa, nw);
-            out.SetAdd(out, tmpa);
-        }
-        if (i < size_g2 && (nw = wg2.Get(i))) {
-            c.wpg128.Get(tmpa, nw);
-            out.SetAdd(out, tmpa);
-        }
-    }
-}
+void ECMult(GroupElemJac &out, const GroupElemJac &a, const Number &an, const Number &gn);
 
 }
 
diff --git a/field.cpp b/field.cpp
new file mode 100644
index 0000000..5f7714c
--- /dev/null
+++ b/field.cpp
@@ -0,0 +1,385 @@
+using namespace std;
+
+#include <assert.h>
+#include <stdint.h>
+#include <string>
+
+#include "num.h"
+#include "field.h"
+
+// #define VERIFY_MAGNITUDE 1
+
+namespace secp256k1 {
+
+/** Implements arithmetic modulo FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFE FFFFFC2F,
+ *  represented as 5 uint64_t's in base 2^52. The values are allowed to contain >52 each. In particular,
+ *  each FieldElem has a 'magnitude' associated with it. Internally, a magnitude M means each element
+ *  is at most M*(2^53-1), except the most significant one, which is limited to M*(2^49-1). All operations
+ *  accept any input with magnitude at most M, and have different rules for propagating magnitude to their
+ *  output.
+ */
+
+FieldElem::FieldElem(int x) {
+    n[0] = x;
+    n[1] = n[2] = n[3] = n[4] = 0;
+#ifdef VERIFY_MAGNITUDE
+    magnitude = 1;
+#endif
+}
+
+FieldElem::FieldElem(const unsigned char *b32) {
+    SetBytes(b32);
+}
+
+void FieldElem::Normalize() {
+    uint64_t c;
+    if (n[0] > 0xFFFFFFFFFFFFFULL || n[1] > 0xFFFFFFFFFFFFFULL || n[2] > 0xFFFFFFFFFFFFFULL || n[3] > 0xFFFFFFFFFFFFFULL || n[4] > 0xFFFFFFFFFFFFULL) {
+        c = n[0];
+        uint64_t t0 = c & 0xFFFFFFFFFFFFFULL;
+        c = (c >> 52) + n[1];
+        uint64_t t1 = c & 0xFFFFFFFFFFFFFULL;
+        c = (c >> 52) + n[2];
+        uint64_t t2 = c & 0xFFFFFFFFFFFFFULL;
+        c = (c >> 52) + n[3];
+        uint64_t t3 = c & 0xFFFFFFFFFFFFFULL;
+        c = (c >> 52) + n[4];
+        uint64_t t4 = c & 0x0FFFFFFFFFFFFULL;
+        c >>= 48;
+        if (c) {
+            c = c * 0x1000003D1ULL + t0;
+            t0 = c & 0xFFFFFFFFFFFFFULL;
+            c = (c >> 52) + t1;
+            t1 = c & 0xFFFFFFFFFFFFFULL;
+            c = (c >> 52) + t2;
+            t2 = c & 0xFFFFFFFFFFFFFULL;
+            c = (c >> 52) + t3;
+            t3 = c & 0xFFFFFFFFFFFFFULL;
+            c = (c >> 52) + t4;
+            t4 = c & 0x0FFFFFFFFFFFFULL;
+            c >>= 48;
+        }
+        n[0] = t0; n[1] = t1; n[2] = t2; n[3] = t3; n[4] = t4;
+    }
+    if (n[4] == 0xFFFFFFFFFFFFULL && n[3] == 0xFFFFFFFFFFFFFULL && n[2] == 0xFFFFFFFFFFFFFULL && n[1] == 0xFFFFFFFFFFFFF && n[0] >= 0xFFFFEFFFFFC2FULL) {
+        n[4] = 0;
+        n[3] = 0;
+        n[2] = 0;
+        n[1] = 0;
+        n[0] -= 0xFFFFEFFFFFC2FULL;
+    }
+#ifdef VERIFY_MAGNITUDE
+    magnitude = 1;
+#endif
+}
+
+bool inline FieldElem::IsZero() {
+    Normalize();
+    return (n[0] == 0 && n[1] == 0 && n[2] == 0 && n[3] == 0 && n[4] == 0);
+}
+
+bool inline operator==(FieldElem &a, FieldElem &b) {
+    a.Normalize();
+    b.Normalize();
+    return (a.n[0] == b.n[0] && a.n[1] == b.n[1] && a.n[2] == b.n[2] && a.n[3] == b.n[3] && a.n[4] == b.n[4]);
+}
+
+void FieldElem::GetBytes(unsigned char *o) {
+    Normalize();
+    for (int i=0; i<32; i++) {
+        int c = 0;
+        for (int j=0; j<2; j++) {
+            int limb = (8*i+4*j)/52;
+            int shift = (8*i+4*j)%52;
+            c |= ((n[limb] >> shift) & 0xF) << (4 * j);
+        }
+        o[31-i] = c;
+    }
+}
+
+void FieldElem::SetBytes(const unsigned char *in) {
+    n[0] = n[1] = n[2] = n[3] = n[4] = 0;
+    for (int i=0; i<32; i++) {
+        for (int j=0; j<2; j++) {
+            int limb = (8*i+4*j)/52;
+            int shift = (8*i+4*j)%52;
+            n[limb] |= (uint64_t)((in[31-i] >> (4*j)) & 0xF) << shift;
+        }
+    }
+#ifdef VERIFY_MAGNITUDE
+    magnitude = 1;
+#endif
+}
+
+void inline FieldElem::SetNeg(const FieldElem &a, int magnitudeIn) {
+#ifdef VERIFY_MAGNITUDE
+    assert(a.magnitude <= magnitudeIn);
+    magnitude = magnitudeIn + 1;
+#endif
+    n[0] = 0xFFFFEFFFFFC2FULL * (magnitudeIn + 1) - a.n[0];
+    n[1] = 0xFFFFFFFFFFFFFULL * (magnitudeIn + 1) - a.n[1];
+    n[2] = 0xFFFFFFFFFFFFFULL * (magnitudeIn + 1) - a.n[2];
+    n[3] = 0xFFFFFFFFFFFFFULL * (magnitudeIn + 1) - a.n[3];
+    n[4] = 0x0FFFFFFFFFFFFULL * (magnitudeIn + 1) - a.n[4];
+}
+
+void inline FieldElem::operator*=(int v) {
+#ifdef VERIFY_MAGNITUDE
+    magnitude *= v;
+#endif
+    n[0] *= v;
+    n[1] *= v;
+    n[2] *= v;
+    n[3] *= v;
+    n[4] *= v;
+}
+
+void inline FieldElem::operator+=(const FieldElem &a) {
+#ifdef VERIFY_MAGNITUDE
+    magnitude += a.magnitude;
+#endif
+    n[0] += a.n[0];
+    n[1] += a.n[1];
+    n[2] += a.n[2];
+    n[3] += a.n[3];
+    n[4] += a.n[4];
+}
+
+void FieldElem::SetMult(const FieldElem &a, const FieldElem &b) {
+#ifdef VERIFY_MAGNITUDE
+    assert(a.magnitude <= 8);
+    assert(b.magnitude <= 8);
+#endif
+    __int128 c = (__int128)a.n[0] * b.n[0];
+    uint64_t t0 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0FFFFFFFFFFFFFE0
+    c = c + (__int128)a.n[0] * b.n[1] +
+            (__int128)a.n[1] * b.n[0];
+    uint64_t t1 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 20000000000000BF
+    c = c + (__int128)a.n[0] * b.n[2] +
+            (__int128)a.n[1] * b.n[1] +
+            (__int128)a.n[2] * b.n[0];
+    uint64_t t2 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 30000000000001A0
+    c = c + (__int128)a.n[0] * b.n[3] +
+            (__int128)a.n[1] * b.n[2] +
+            (__int128)a.n[2] * b.n[1] +
+            (__int128)a.n[3] * b.n[0];
+    uint64_t t3 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 4000000000000280
+    c = c + (__int128)a.n[0] * b.n[4] +
+            (__int128)a.n[1] * b.n[3] +
+            (__int128)a.n[2] * b.n[2] +
+            (__int128)a.n[3] * b.n[1] +
+            (__int128)a.n[4] * b.n[0];
+    uint64_t t4 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 320000000000037E
+    c = c + (__int128)a.n[1] * b.n[4] +
+            (__int128)a.n[2] * b.n[3] +
+            (__int128)a.n[3] * b.n[2] +
+            (__int128)a.n[4] * b.n[1];
+    uint64_t t5 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 22000000000002BE
+    c = c + (__int128)a.n[2] * b.n[4] +
+            (__int128)a.n[3] * b.n[3] +
+            (__int128)a.n[4] * b.n[2];
+    uint64_t t6 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 12000000000001DE
+    c = c + (__int128)a.n[3] * b.n[4] +
+            (__int128)a.n[4] * b.n[3];
+    uint64_t t7 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 02000000000000FE
+    c = c + (__int128)a.n[4] * b.n[4];
+    uint64_t t8 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 001000000000001E
+    uint64_t t9 = c;
+    c = t0 + (__int128)t5 * 0x1000003D10ULL;
+
+    t0 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
+    c = c + t1 + (__int128)t6 * 0x1000003D10ULL;
+    t1 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
+    c = c + t2 + (__int128)t7 * 0x1000003D10ULL;
+    n[2] = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
+    c = c + t3 + (__int128)t8 * 0x1000003D10ULL;
+    n[3] = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
+    c = c + t4 + (__int128)t9 * 0x1000003D10ULL;
+    n[4] = c & 0x0FFFFFFFFFFFFULL; c = c >> 48; // c max 000001000003D110
+    c = t0 + (__int128)c * 0x1000003D1ULL;
+    n[0] = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 1000008
+    n[1] = t1 + c;
+#ifdef VERIFY_MAGNITUDE
+    magnitude = 1;
+#endif
+}
+
+void FieldElem::SetSquare(const FieldElem &a) {
+#ifdef VERIFY_MAGNITUDE
+    assert(a.magnitude <= 8);
+#endif
+    __int128 c = (__int128)a.n[0] * a.n[0];
+    uint64_t t0 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0FFFFFFFFFFFFFE0
+    c = c + (__int128)(a.n[0]*2) * a.n[1];
+    uint64_t t1 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 20000000000000BF
+    c = c + (__int128)(a.n[0]*2) * a.n[2] +
+            (__int128)a.n[1] * a.n[1];
+    uint64_t t2 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 30000000000001A0
+    c = c + (__int128)(a.n[0]*2) * a.n[3] +
+            (__int128)(a.n[1]*2) * a.n[2];
+    uint64_t t3 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 4000000000000280
+    c = c + (__int128)(a.n[0]*2) * a.n[4] +
+            (__int128)(a.n[1]*2) * a.n[3] +
+            (__int128)a.n[2] * a.n[2];
+    uint64_t t4 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 320000000000037E
+    c = c + (__int128)(a.n[1]*2) * a.n[4] +
+            (__int128)(a.n[2]*2) * a.n[3];
+    uint64_t t5 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 22000000000002BE
+    c = c + (__int128)(a.n[2]*2) * a.n[4] +
+            (__int128)a.n[3] * a.n[3];
+    uint64_t t6 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 12000000000001DE
+    c = c + (__int128)(a.n[3]*2) * a.n[4];
+    uint64_t t7 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 02000000000000FE
+    c = c + (__int128)a.n[4] * a.n[4];
+    uint64_t t8 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 001000000000001E
+    uint64_t t9 = c;
+    c = t0 + (__int128)t5 * 0x1000003D10ULL;
+    t0 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
+    c = c + t1 + (__int128)t6 * 0x1000003D10ULL;
+    t1 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
+    c = c + t2 + (__int128)t7 * 0x1000003D10ULL;
+    n[2] = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
+    c = c + t3 + (__int128)t8 * 0x1000003D10ULL;
+    n[3] = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
+    c = c + t4 + (__int128)t9 * 0x1000003D10ULL;
+    n[4] = c & 0x0FFFFFFFFFFFFULL; c = c >> 48; // c max 000001000003D110
+    c = t0 + (__int128)c * 0x1000003D1ULL;
+    n[0] = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 1000008
+    n[1] = t1 + c;
+#ifdef VERIFY_MAGNITUDE
+    assert(a.magnitude <= 8);
+#endif
+}
+
+void FieldElem::SetSquareRoot(const FieldElem &a) {
+    // calculate a^p, with p={15,780,1022,1023}
+    FieldElem a2; a2.SetSquare(a);
+    FieldElem a3; a3.SetMult(a2,a);
+    FieldElem a6; a6.SetSquare(a3);
+    FieldElem a12; a12.SetSquare(a6);
+    FieldElem a15; a15.SetMult(a12,a3);
+    FieldElem a30; a30.SetSquare(a15);
+    FieldElem a60; a60.SetSquare(a30);
+    FieldElem a120; a120.SetSquare(a60);
+    FieldElem a240; a240.SetSquare(a120);
+    FieldElem a255; a255.SetMult(a240,a15);
+    FieldElem a510; a510.SetSquare(a255);
+    FieldElem a750; a750.SetMult(a510,a240);
+    FieldElem a780; a780.SetMult(a750,a30);
+    FieldElem a1020; a1020.SetSquare(a510);
+    FieldElem a1022; a1022.SetMult(a1020,a2);
+    FieldElem a1023; a1023.SetMult(a1022,a);
+    FieldElem x = a15;
+    for (int i=0; i<21; i++) {
+        for (int j=0; j<10; j++) x.SetSquare(x);
+        x.SetMult(x,a1023);
+    }
+    for (int j=0; j<10; j++) x.SetSquare(x);
+    x.SetMult(x,a1022);
+    for (int i=0; i<2; i++) {
+        for (int j=0; j<10; j++) x.SetSquare(x);
+        x.SetMult(x,a1023);
+    }
+    for (int j=0; j<10; j++) x.SetSquare(x);
+    SetMult(x,a780);
+}
+
+bool FieldElem::IsOdd() {
+    Normalize();
+    return n[0] & 1;
+}
+
+std::string FieldElem::ToString() {
+    unsigned char tmp[32];
+    GetBytes(tmp);
+    std::string ret;
+    for (int i=0; i<32; i++) {
+        static const char *c = "0123456789ABCDEF";
+        ret += c[(tmp[i] >> 4) & 0xF];
+        ret += c[(tmp[i]) & 0xF];
+    }
+    return ret;
+}
+
+void FieldElem::SetHex(const std::string &str) {
+    unsigned char tmp[32] = {};
+    static const int cvt[256] = {0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
+                                 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
+                                 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
+                                 0, 1, 2, 3, 4, 5, 6,7,8,9,0,0,0,0,0,0,
+                                 0,10,11,12,13,14,15,0,0,0,0,0,0,0,0,0,
+                                 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
+                                 0,10,11,12,13,14,15,0,0,0,0,0,0,0,0,0,
+                                 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
+                                 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
+                                 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
+                                 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
+                                 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
+                                 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
+                                 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
+                                 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
+                                 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0};
+    for (unsigned int i=0; i<32; i++) {
+        if (str.length() > i*2)
+            tmp[32 - str.length()/2 + i] = (cvt[(unsigned char)str[2*i]] << 4) + cvt[(unsigned char)str[2*i+1]];
+    }
+    SetBytes(tmp);
+}
+
+static const unsigned char field_p_[] = {0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
+                                         0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
+                                         0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
+                                         0xFF,0xFF,0xFF,0xFE,0xFF,0xFF,0xFC,0x2F};
+
+FieldConstants::FieldConstants() : field_p(field_p_, sizeof(field_p_)) {}
+
+const FieldConstants &GetFieldConst() {
+    static const FieldConstants field_const;
+    return field_const;
+}
+
+void FieldElem::SetInverse(FieldElem &a) {
+#if defined(USE_FIELDINVERSE_BUILTIN)
+    // calculate a^p, with p={45,63,1019,1023}
+    FieldElem a2; a2.SetSquare(a);
+    FieldElem a3; a3.SetMult(a2,a);
+    FieldElem a4; a4.SetSquare(a2);
+    FieldElem a5; a5.SetMult(a4,a);
+    FieldElem a10; a10.SetSquare(a5);
+    FieldElem a11; a11.SetMult(a10,a);
+    FieldElem a21; a21.SetMult(a11,a10);
+    FieldElem a42; a42.SetSquare(a21);
+    FieldElem a45; a45.SetMult(a42,a3);
+    FieldElem a63; a63.SetMult(a42,a21);
+    FieldElem a126; a126.SetSquare(a63);
+    FieldElem a252; a252.SetSquare(a126);
+    FieldElem a504; a504.SetSquare(a252);
+    FieldElem a1008; a1008.SetSquare(a504);
+    FieldElem a1019; a1019.SetMult(a1008,a11);
+    FieldElem a1023; a1023.SetMult(a1019,a4);
+    FieldElem x = a63;
+    for (int i=0; i<21; i++) {
+        for (int j=0; j<10; j++) x.SetSquare(x);
+        x.SetMult(x,a1023);
+    }
+    for (int j=0; j<10; j++) x.SetSquare(x);
+    x.SetMult(x,a1019);
+    for (int i=0; i<2; i++) {
+        for (int j=0; j<10; j++) x.SetSquare(x);
+        x.SetMult(x,a1023);
+    }
+    for (int j=0; j<10; j++) x.SetSquare(x);
+    SetMult(x,a45);
+#else
+    unsigned char b[32];
+    a.GetBytes(b);
+    {
+        const Number &p = GetFieldConst().field_p;
+        Number n; n.SetBytes(b, 32);
+        n.SetModInverse(n, p);
+        n.GetBytes(b, 32);
+    }
+    SetBytes(b);
+#endif
+}
+
+}
diff --git a/field.h b/field.h
index 6b21ce7..0584dba 100644
--- a/field.h
+++ b/field.h
@@ -3,7 +3,6 @@
 
 using namespace std;
 
-#include <assert.h>
 #include <stdint.h>
 #include <string>
 
@@ -22,7 +21,7 @@ namespace secp256k1 {
  */
 class FieldElem {
 private:
-    // X = sum(i=0..4, elem[i]*2^52)
+    // X = sum(i=0..4, elem[i]*2^52) mod n
     uint64_t n[5];
 #ifdef VERIFY_MAGNITUDE
     int magnitude;
@@ -31,384 +30,58 @@ private:
 public:
 
     /** Creates a constant field element. Magnitude=1 */
-    FieldElem(int x = 0) {
-        n[0] = x;
-        n[1] = n[2] = n[3] = n[4] = 0;
-#ifdef VERIFY_MAGNITUDE
-        magnitude = 1;
-#endif
-    }
+    FieldElem(int x = 0);
 
-    FieldElem(const unsigned char *b32) {
-        SetBytes(b32);
-    }
+    FieldElem(const unsigned char *b32);
 
     /** Normalizes the internal representation entries. Magnitude=1 */
-    void Normalize() {
-        uint64_t c;
-        if (n[0] > 0xFFFFFFFFFFFFFULL || n[1] > 0xFFFFFFFFFFFFFULL || n[2] > 0xFFFFFFFFFFFFFULL || n[3] > 0xFFFFFFFFFFFFFULL || n[4] > 0xFFFFFFFFFFFFULL) {
-            c = n[0];
-            uint64_t t0 = c & 0xFFFFFFFFFFFFFULL;
-            c = (c >> 52) + n[1];
-            uint64_t t1 = c & 0xFFFFFFFFFFFFFULL;
-            c = (c >> 52) + n[2];
-            uint64_t t2 = c & 0xFFFFFFFFFFFFFULL;
-            c = (c >> 52) + n[3];
-            uint64_t t3 = c & 0xFFFFFFFFFFFFFULL;
-            c = (c >> 52) + n[4];
-            uint64_t t4 = c & 0x0FFFFFFFFFFFFULL;
-            c >>= 48;
-            if (c) {
-                c = c * 0x1000003D1ULL + t0;
-                t0 = c & 0xFFFFFFFFFFFFFULL;
-                c = (c >> 52) + t1;
-                t1 = c & 0xFFFFFFFFFFFFFULL;
-                c = (c >> 52) + t2;
-                t2 = c & 0xFFFFFFFFFFFFFULL;
-                c = (c >> 52) + t3;
-                t3 = c & 0xFFFFFFFFFFFFFULL;
-                c = (c >> 52) + t4;
-                t4 = c & 0x0FFFFFFFFFFFFULL;
-                c >>= 48;
-            }
-            n[0] = t0; n[1] = t1; n[2] = t2; n[3] = t3; n[4] = t4;
-        }
-        if (n[4] == 0xFFFFFFFFFFFFULL && n[3] == 0xFFFFFFFFFFFFFULL && n[2] == 0xFFFFFFFFFFFFFULL && n[1] == 0xFFFFFFFFFFFFF && n[0] >= 0xFFFFEFFFFFC2FULL) {
-            n[4] = 0;
-            n[3] = 0;
-            n[2] = 0;
-            n[1] = 0;
-            n[0] -= 0xFFFFEFFFFFC2FULL;
-        }
-#ifdef VERIFY_MAGNITUDE
-        magnitude = 1;
-#endif
-    }
+    void Normalize();
 
-    bool IsZero() {
-        Normalize();
-        return (n[0] == 0 && n[1] == 0 && n[2] == 0 && n[3] == 0 && n[4] == 0);
-    }
+    bool IsZero();
 
-    bool friend operator==(FieldElem &a, FieldElem &b) {
-        a.Normalize();
-        b.Normalize();
-        return (a.n[0] == b.n[0] && a.n[1] == b.n[1] && a.n[2] == b.n[2] && a.n[3] == b.n[3] && a.n[4] == b.n[4]);
-    }
+    bool friend operator==(FieldElem &a, FieldElem &b);
 
     /** extract as 32-byte big endian array */
-    void GetBytes(unsigned char *o) {
-        Normalize();
-        for (int i=0; i<32; i++) {
-            int c = 0;
-            for (int j=0; j<2; j++) {
-                int limb = (8*i+4*j)/52;
-                int shift = (8*i+4*j)%52;
-                c |= ((n[limb] >> shift) & 0xF) << (4 * j);
-            }
-            o[31-i] = c;
-        }
-    }
+    void GetBytes(unsigned char *o);
 
     /** set value of 32-byte big endian array */
-    void SetBytes(const unsigned char *in) {
-        n[0] = n[1] = n[2] = n[3] = n[4] = 0;
-        for (int i=0; i<32; i++) {
-            for (int j=0; j<2; j++) {
-                int limb = (8*i+4*j)/52;
-                int shift = (8*i+4*j)%52;
-                n[limb] |= (uint64_t)((in[31-i] >> (4*j)) & 0xF) << shift;
-            }
-        }
-#ifdef VERIFY_MAGNITUDE
-        magnitude = 1;
-#endif
-    }
+    void SetBytes(const unsigned char *in);
 
     /** Set a FieldElem to be the negative of another. Increases magnitude by one. */
-    void SetNeg(const FieldElem &a, int magnitudeIn) {
-#ifdef VERIFY_MAGNITUDE
-        assert(a.magnitude <= magnitudeIn);
-        magnitude = magnitudeIn + 1;
-#endif
-        n[0] = 0xFFFFEFFFFFC2FULL * (magnitudeIn + 1) - a.n[0];
-        n[1] = 0xFFFFFFFFFFFFFULL * (magnitudeIn + 1) - a.n[1];
-        n[2] = 0xFFFFFFFFFFFFFULL * (magnitudeIn + 1) - a.n[2];
-        n[3] = 0xFFFFFFFFFFFFFULL * (magnitudeIn + 1) - a.n[3];
-        n[4] = 0x0FFFFFFFFFFFFULL * (magnitudeIn + 1) - a.n[4];
-    }
+    void SetNeg(const FieldElem &a, int magnitudeIn);
 
     /** Multiplies this FieldElem with an integer constant. Magnitude is multiplied by v */
-    void operator*=(int v) {
-#ifdef VERIFY_MAGNITUDE
-        magnitude *= v;
-#endif
-        n[0] *= v;
-        n[1] *= v;
-        n[2] *= v;
-        n[3] *= v;
-        n[4] *= v;
-    }
+    void operator*=(int v);
 
-    void operator+=(const FieldElem &a) {
-#ifdef VERIFY_MAGNITUDE
-        magnitude += a.magnitude;
-#endif
-        n[0] += a.n[0];
-        n[1] += a.n[1];
-        n[2] += a.n[2];
-        n[3] += a.n[3];
-        n[4] += a.n[4];
-    }
+    void operator+=(const FieldElem &a);
 
     /** Set this FieldElem to be the multiplication of two others. Magnitude=1 */
-    void SetMult(const FieldElem &a, const FieldElem &b) {
-#ifdef VERIFY_MAGNITUDE
-        assert(a.magnitude <= 8);
-        assert(b.magnitude <= 8);
-#endif
-        __int128 c = (__int128)a.n[0] * b.n[0];
-        uint64_t t0 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0FFFFFFFFFFFFFE0
-        c = c + (__int128)a.n[0] * b.n[1] +
-                (__int128)a.n[1] * b.n[0];
-        uint64_t t1 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 20000000000000BF
-        c = c + (__int128)a.n[0] * b.n[2] +
-                (__int128)a.n[1] * b.n[1] +
-                (__int128)a.n[2] * b.n[0];
-        uint64_t t2 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 30000000000001A0
-        c = c + (__int128)a.n[0] * b.n[3] +
-                (__int128)a.n[1] * b.n[2] +
-                (__int128)a.n[2] * b.n[1] +
-                (__int128)a.n[3] * b.n[0];
-        uint64_t t3 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 4000000000000280
-        c = c + (__int128)a.n[0] * b.n[4] +
-                (__int128)a.n[1] * b.n[3] +
-                (__int128)a.n[2] * b.n[2] +
-                (__int128)a.n[3] * b.n[1] +
-                (__int128)a.n[4] * b.n[0];
-        uint64_t t4 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 320000000000037E
-        c = c + (__int128)a.n[1] * b.n[4] +
-                (__int128)a.n[2] * b.n[3] +
-                (__int128)a.n[3] * b.n[2] +
-                (__int128)a.n[4] * b.n[1];
-        uint64_t t5 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 22000000000002BE
-        c = c + (__int128)a.n[2] * b.n[4] +
-                (__int128)a.n[3] * b.n[3] +
-                (__int128)a.n[4] * b.n[2];
-        uint64_t t6 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 12000000000001DE
-        c = c + (__int128)a.n[3] * b.n[4] +
-                (__int128)a.n[4] * b.n[3];
-        uint64_t t7 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 02000000000000FE
-        c = c + (__int128)a.n[4] * b.n[4];
-        uint64_t t8 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 001000000000001E
-        uint64_t t9 = c;
-        c = t0 + (__int128)t5 * 0x1000003D10ULL;
-        t0 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
-        c = c + t1 + (__int128)t6 * 0x1000003D10ULL;
-        t1 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
-        c = c + t2 + (__int128)t7 * 0x1000003D10ULL;
-        n[2] = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
-        c = c + t3 + (__int128)t8 * 0x1000003D10ULL;
-        n[3] = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
-        c = c + t4 + (__int128)t9 * 0x1000003D10ULL;
-        n[4] = c & 0x0FFFFFFFFFFFFULL; c = c >> 48; // c max 000001000003D110
-        c = t0 + (__int128)c * 0x1000003D1ULL;
-        n[0] = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 1000008
-        n[1] = t1 + c;
-#ifdef VERIFY_MAGNITUDE
-        magnitude = 1;
-#endif
-    }
+    void SetMult(const FieldElem &a, const FieldElem &b);
 
     /** Set this FieldElem to be the square of another. Magnitude=1 */
-    void SetSquare(const FieldElem &a) {
-#ifdef VERIFY_MAGNITUDE
-        assert(a.magnitude <= 8);
-#endif
-        __int128 c = (__int128)a.n[0] * a.n[0];
-        uint64_t t0 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0FFFFFFFFFFFFFE0
-        c = c + (__int128)(a.n[0]*2) * a.n[1];
-        uint64_t t1 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 20000000000000BF
-        c = c + (__int128)(a.n[0]*2) * a.n[2] +
-                (__int128)a.n[1] * a.n[1];
-        uint64_t t2 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 30000000000001A0
-        c = c + (__int128)(a.n[0]*2) * a.n[3] +
-                (__int128)(a.n[1]*2) * a.n[2];
-        uint64_t t3 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 4000000000000280
-        c = c + (__int128)(a.n[0]*2) * a.n[4] +
-                (__int128)(a.n[1]*2) * a.n[3] +
-                (__int128)a.n[2] * a.n[2];
-        uint64_t t4 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 320000000000037E
-        c = c + (__int128)(a.n[1]*2) * a.n[4] +
-                (__int128)(a.n[2]*2) * a.n[3];
-        uint64_t t5 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 22000000000002BE
-        c = c + (__int128)(a.n[2]*2) * a.n[4] +
-                (__int128)a.n[3] * a.n[3];
-        uint64_t t6 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 12000000000001DE
-        c = c + (__int128)(a.n[3]*2) * a.n[4];
-        uint64_t t7 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 02000000000000FE
-        c = c + (__int128)a.n[4] * a.n[4];
-        uint64_t t8 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 001000000000001E
-        uint64_t t9 = c;
-        c = t0 + (__int128)t5 * 0x1000003D10ULL;
-        t0 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
-        c = c + t1 + (__int128)t6 * 0x1000003D10ULL;
-        t1 = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
-        c = c + t2 + (__int128)t7 * 0x1000003D10ULL;
-        n[2] = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
-        c = c + t3 + (__int128)t8 * 0x1000003D10ULL;
-        n[3] = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 0000001000003D10
-        c = c + t4 + (__int128)t9 * 0x1000003D10ULL;
-        n[4] = c & 0x0FFFFFFFFFFFFULL; c = c >> 48; // c max 000001000003D110
-        c = t0 + (__int128)c * 0x1000003D1ULL;
-        n[0] = c & 0xFFFFFFFFFFFFFULL; c = c >> 52; // c max 1000008
-        n[1] = t1 + c;
-#ifdef VERIFY_MAGNITUDE
-        assert(a.magnitude <= 8);
-#endif
-    }
+    void SetSquare(const FieldElem &a);
 
     /** Set this to be the (modular) square root of another FieldElem. Magnitude=1 */
-    void SetSquareRoot(const FieldElem &a) {
-        // calculate a^p, with p={15,780,1022,1023}
-        FieldElem a2; a2.SetSquare(a);
-        FieldElem a3; a3.SetMult(a2,a);
-        FieldElem a6; a6.SetSquare(a3);
-        FieldElem a12; a12.SetSquare(a6);
-        FieldElem a15; a15.SetMult(a12,a3);
-        FieldElem a30; a30.SetSquare(a15);
-        FieldElem a60; a60.SetSquare(a30);
-        FieldElem a120; a120.SetSquare(a60);
-        FieldElem a240; a240.SetSquare(a120);
-        FieldElem a255; a255.SetMult(a240,a15);
-        FieldElem a510; a510.SetSquare(a255);
-        FieldElem a750; a750.SetMult(a510,a240);
-        FieldElem a780; a780.SetMult(a750,a30);
-        FieldElem a1020; a1020.SetSquare(a510);
-        FieldElem a1022; a1022.SetMult(a1020,a2);
-        FieldElem a1023; a1023.SetMult(a1022,a);
-        FieldElem x = a15;
-        for (int i=0; i<21; i++) {
-            for (int j=0; j<10; j++) x.SetSquare(x);
-            x.SetMult(x,a1023);
-        }
-        for (int j=0; j<10; j++) x.SetSquare(x);
-        x.SetMult(x,a1022);
-        for (int i=0; i<2; i++) {
-            for (int j=0; j<10; j++) x.SetSquare(x);
-            x.SetMult(x,a1023);
-        }
-        for (int j=0; j<10; j++) x.SetSquare(x);
-        SetMult(x,a780);
-    }
+    void SetSquareRoot(const FieldElem &a);
 
-    bool IsOdd() {
-        Normalize();
-        return n[0] & 1;
-    }
+    bool IsOdd();
 
     /** Set this to be the (modular) inverse of another FieldElem. Magnitude=1 */
     void SetInverse(FieldElem &a);
 
-    std::string ToString() {
-        unsigned char tmp[32];
-        GetBytes(tmp);
-        std::string ret;
-        for (int i=0; i<32; i++) {
-            static const char *c = "0123456789ABCDEF";
-            ret += c[(tmp[i] >> 4) & 0xF];
-            ret += c[(tmp[i]) & 0xF];
-        }
-        return ret;
-    }
+    std::string ToString();
 
-    void SetHex(const std::string &str) {
-        unsigned char tmp[32] = {};
-        static const int cvt[256] = {0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
-                                     0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
-                                     0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
-                                     0, 1, 2, 3, 4, 5, 6,7,8,9,0,0,0,0,0,0,
-                                     0,10,11,12,13,14,15,0,0,0,0,0,0,0,0,0,
-                                     0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
-                                     0,10,11,12,13,14,15,0,0,0,0,0,0,0,0,0,
-                                     0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
-                                     0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
-                                     0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
-                                     0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
-                                     0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
-                                     0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
-                                     0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
-                                     0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0,
-                                     0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0,0,0,0};
-        for (unsigned int i=0; i<32; i++) {
-            if (str.length() > i*2)
-                tmp[32 - str.length()/2 + i] = (cvt[(unsigned char)str[2*i]] << 4) + cvt[(unsigned char)str[2*i+1]];
-        }
-        SetBytes(tmp);
-    }
+    void SetHex(const std::string &str);
 };
 
-static const unsigned char field_p_[] = {0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
-                                         0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
-                                         0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
-                                         0xFF,0xFF,0xFF,0xFE,0xFF,0xFF,0xFC,0x2F};
-
 class FieldConstants {
 public:
     const Number field_p;
 
-    FieldConstants() : field_p(field_p_, sizeof(field_p_)) {}
+    FieldConstants();
 };
 
-const FieldConstants &GetFieldConst() {
-    static const FieldConstants field_const;
-    return field_const;
-}
-
-void FieldElem::SetInverse(FieldElem &a) {
-#if defined(USE_FIELDINVERSE_BUILTIN)
-    // calculate a^p, with p={45,63,1019,1023}
-    FieldElem a2; a2.SetSquare(a);
-    FieldElem a3; a3.SetMult(a2,a);
-    FieldElem a4; a4.SetSquare(a2);
-    FieldElem a5; a5.SetMult(a4,a);
-    FieldElem a10; a10.SetSquare(a5);
-    FieldElem a11; a11.SetMult(a10,a);
-    FieldElem a21; a21.SetMult(a11,a10);
-    FieldElem a42; a42.SetSquare(a21);
-    FieldElem a45; a45.SetMult(a42,a3);
-    FieldElem a63; a63.SetMult(a42,a21);
-    FieldElem a126; a126.SetSquare(a63);
-    FieldElem a252; a252.SetSquare(a126);
-    FieldElem a504; a504.SetSquare(a252);
-    FieldElem a1008; a1008.SetSquare(a504);
-    FieldElem a1019; a1019.SetMult(a1008,a11);
-    FieldElem a1023; a1023.SetMult(a1019,a4);
-    FieldElem x = a63;
-    for (int i=0; i<21; i++) {
-        for (int j=0; j<10; j++) x.SetSquare(x);
-        x.SetMult(x,a1023);
-    }
-    for (int j=0; j<10; j++) x.SetSquare(x);
-    x.SetMult(x,a1019);
-    for (int i=0; i<2; i++) {
-        for (int j=0; j<10; j++) x.SetSquare(x);
-        x.SetMult(x,a1023);
-    }
-    for (int j=0; j<10; j++) x.SetSquare(x);
-    SetMult(x,a45);
-#else
-    unsigned char b[32];
-    a.GetBytes(b);
-    {
-        const Number &p = GetFieldConst().field_p;
-        Number n; n.SetBytes(b, 32);
-        n.SetModInverse(n, p);
-        n.GetBytes(b, 32);
-    }
-    SetBytes(b);
-#endif
-}
+const FieldConstants &GetFieldConst();
 
 }
 
diff --git a/group.cpp b/group.cpp
new file mode 100644
index 0000000..af38b5d
--- /dev/null
+++ b/group.cpp
@@ -0,0 +1,317 @@
+#include <string>
+
+#include "num.h"
+#include "field.h"
+#include "group.h"
+
+namespace secp256k1 {
+
+GroupElem::GroupElem() {
+    fInfinity = true;
+}
+
+GroupElem::GroupElem(const FieldElem &xin, const FieldElem &yin) {
+    fInfinity = false;
+    x = xin;
+    y = yin;
+}
+
+bool GroupElem::IsInfinity() const {
+    return fInfinity;
+}
+
+void GroupElem::SetNeg(const GroupElem &p) {
+    *this = p;
+    y.Normalize();
+    y.SetNeg(y, 1);
+}
+
+void GroupElem::GetX(FieldElem &xout) {
+    xout = x;
+}
+
+void GroupElem::GetY(FieldElem &yout) {
+    yout = y;
+}
+
+std::string GroupElem::ToString() const {
+    if (fInfinity)
+        return "(inf)";
+    FieldElem xc = x, yc = y;
+    return "(" + xc.ToString() + "," + yc.ToString() + ")";
+}
+
+GroupElemJac::GroupElemJac() : GroupElem(), z(1) {}
+
+GroupElemJac::GroupElemJac(const FieldElem &xin, const FieldElem &yin) : GroupElem(xin,yin), z(1) {}
+
+GroupElemJac::GroupElemJac(const GroupElem &in) : GroupElem(in), z(1) {}
+
+void GroupElemJac::SetJac(GroupElemJac &jac) {
+    *this = jac;
+}
+
+bool GroupElemJac::IsValid() const {
+    if (IsInfinity())
+        return false;
+    // y^2 = x^3 + 7
+    // (Y/Z^3)^2 = (X/Z^2)^3 + 7
+    // Y^2 / Z^6 = X^3 / Z^6 + 7
+    // Y^2 = X^3 + 7*Z^6
+    FieldElem y2; y2.SetSquare(y);
+    FieldElem x3; x3.SetSquare(x); x3.SetMult(x3,x);
+    FieldElem z2; z2.SetSquare(z);
+    FieldElem z6; z6.SetSquare(z2); z6.SetMult(z6,z2);
+    z6 *= 7;
+    x3 += z6;
+    return y2 == x3;
+}
+
+void GroupElemJac::GetAffine(GroupElem &aff) {
+    z.SetInverse(z);
+    FieldElem z2;
+    z2.SetSquare(z);
+    FieldElem z3;
+    z3.SetMult(z,z2);
+    x.SetMult(x,z2);
+    y.SetMult(y,z3);
+    z = FieldElem(1);
+    aff.fInfinity = fInfinity;
+    aff.x = x;
+    aff.y = y;
+}
+
+void GroupElemJac::GetX(FieldElem &xout) {
+    FieldElem zi;
+    zi.SetInverse(z);
+    zi.SetSquare(zi);
+    xout.SetMult(x, zi);
+}
+
+void GroupElemJac::GetY(FieldElem &yout) {
+    FieldElem zi;
+    zi.SetInverse(z);
+    FieldElem zi3; zi3.SetSquare(zi); zi3.SetMult(zi, zi3);
+    yout.SetMult(y, zi3);
+}
+
+bool GroupElemJac::IsInfinity() const {
+    return fInfinity;
+}
+
+
+void GroupElemJac::SetNeg(const GroupElemJac &p) {
+    *this = p;
+    y.Normalize();
+    y.SetNeg(y, 1);
+}
+
+void GroupElemJac::SetCompressed(const FieldElem &xin, bool fOdd) {
+    x = xin;
+    FieldElem x2; x2.SetSquare(x);
+    FieldElem x3; x3.SetMult(x,x2);
+    fInfinity = false;
+    FieldElem c(7);
+    c += x3;
+    y.SetSquareRoot(c);
+    z = FieldElem(1);
+    if (y.IsOdd() != fOdd)
+        y.SetNeg(y,1);
+}
+
+void GroupElemJac::SetDouble(const GroupElemJac &p) {
+    FieldElem t5 = p.y;
+    if (p.fInfinity || t5.IsZero()) {
+        fInfinity = true;
+        return;
+    }
+
+    FieldElem t1,t2,t3,t4;
+    z.SetMult(t5,p.z);
+    z *= 2;                // Z' = 2*Y*Z (2)
+    t1.SetSquare(p.x);
+    t1 *= 3;               // T1 = 3*X^2 (3)
+    t2.SetSquare(t1);      // T2 = 9*X^4 (1)
+    t3.SetSquare(t5);
+    t3 *= 2;               // T3 = 2*Y^2 (2)
+    t4.SetSquare(t3);
+    t4 *= 2;               // T4 = 8*Y^4 (2)
+    t3.SetMult(p.x,t3);      // T3 = 2*X*Y^2 (1)
+    x = t3;
+    x *= 4;                // X' = 8*X*Y^2 (4)
+    x.SetNeg(x,4);         // X' = -8*X*Y^2 (5)
+    x += t2;               // X' = 9*X^4 - 8*X*Y^2 (6)
+    t2.SetNeg(t2,1);       // T2 = -9*X^4 (2)
+    t3 *= 6;               // T3 = 12*X*Y^2 (6)
+    t3 += t2;              // T3 = 12*X*Y^2 - 9*X^4 (8)
+    y.SetMult(t1,t3);      // Y' = 36*X^3*Y^2 - 27*X^6 (1)
+    t2.SetNeg(t4,2);       // T2 = -8*Y^4 (3)
+    y += t2;               // Y' = 36*X^3*Y^2 - 27*X^6 - 8*Y^4 (4)
+    fInfinity = false;
+}
+
+void GroupElemJac::SetAdd(const GroupElemJac &p, const GroupElemJac &q) {
+    if (p.fInfinity) {
+        *this = q;
+        return;
+    }
+    if (q.fInfinity) {
+        *this = p;
+        return;
+    }
+    fInfinity = false;
+    const FieldElem &x1 = p.x, &y1 = p.y, &z1 = p.z, &x2 = q.x, &y2 = q.y, &z2 = q.z;
+    FieldElem z22; z22.SetSquare(z2);
+    FieldElem z12; z12.SetSquare(z1);
+    FieldElem u1; u1.SetMult(x1, z22);
+    FieldElem u2; u2.SetMult(x2, z12);
+    FieldElem s1; s1.SetMult(y1, z22); s1.SetMult(s1, z2);
+    FieldElem s2; s2.SetMult(y2, z12); s2.SetMult(s2, z1);
+    if (u1 == u2) {
+        if (s1 == s2) {
+            SetDouble(p);
+        } else {
+            fInfinity = true;
+        }
+        return;
+    }
+    FieldElem h; h.SetNeg(u1,1); h += u2;
+    FieldElem r; r.SetNeg(s1,1); r += s2;
+    FieldElem r2; r2.SetSquare(r);
+    FieldElem h2; h2.SetSquare(h);
+    FieldElem h3; h3.SetMult(h,h2);
+    z.SetMult(z1,z2); z.SetMult(z, h);
+    FieldElem t; t.SetMult(u1,h2);
+    x = t; x *= 2; x += h3; x.SetNeg(x,3); x += r2;
+    y.SetNeg(x,5); y += t; y.SetMult(y,r);
+    h3.SetMult(h3,s1); h3.SetNeg(h3,1);
+    y += h3;
+}
+
+void GroupElemJac::SetAdd(const GroupElemJac &p, const GroupElem &q) {
+    if (p.fInfinity) {
+        x = q.x;
+        y = q.y;
+        fInfinity = q.fInfinity;
+        z = FieldElem(1);
+        return;
+    }
+    if (q.fInfinity) {
+        *this = p;
+        return;
+    }
+    fInfinity = false;
+    const FieldElem &x1 = p.x, &y1 = p.y, &z1 = p.z, &x2 = q.x, &y2 = q.y;
+    FieldElem z12; z12.SetSquare(z1);
+    FieldElem u1 = x1; u1.Normalize();
+    FieldElem u2; u2.SetMult(x2, z12);
+    FieldElem s1 = y1; s1.Normalize(); 
+    FieldElem s2; s2.SetMult(y2, z12); s2.SetMult(s2, z1);
+    if (u1 == u2) {
+        if (s1 == s2) {
+            SetDouble(p);
+        } else {
+            fInfinity = true;
+        }
+        return;
+    }
+    FieldElem h; h.SetNeg(u1,1); h += u2;
+    FieldElem r; r.SetNeg(s1,1); r += s2;
+    FieldElem r2; r2.SetSquare(r);
+    FieldElem h2; h2.SetSquare(h);
+    FieldElem h3; h3.SetMult(h,h2);
+    z = p.z; z.SetMult(z, h);
+    FieldElem t; t.SetMult(u1,h2);
+    x = t; x *= 2; x += h3; x.SetNeg(x,3); x += r2;
+    y.SetNeg(x,5); y += t; y.SetMult(y,r);
+    h3.SetMult(h3,s1); h3.SetNeg(h3,1);
+    y += h3;
+}
+
+std::string GroupElemJac::ToString() const {
+    GroupElemJac cop = *this;
+    GroupElem aff;
+    cop.GetAffine(aff);
+    return aff.ToString();
+}
+
+void GroupElem::SetJac(GroupElemJac &jac) {
+    jac.GetAffine(*this);
+}
+
+static const unsigned char order_[] = {0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
+                                       0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFE,
+                                       0xBA,0xAE,0xDC,0xE6,0xAF,0x48,0xA0,0x3B,
+                                       0xBF,0xD2,0x5E,0x8C,0xD0,0x36,0x41,0x41};
+
+static const unsigned char g_x_[] = {0x79,0xBE,0x66,0x7E,0xF9,0xDC,0xBB,0xAC,
+                                     0x55,0xA0,0x62,0x95,0xCE,0x87,0x0B,0x07,
+                                     0x02,0x9B,0xFC,0xDB,0x2D,0xCE,0x28,0xD9,
+                                     0x59,0xF2,0x81,0x5B,0x16,0xF8,0x17,0x98};
+
+static const unsigned char g_y_[] = {0x48,0x3A,0xDA,0x77,0x26,0xA3,0xC4,0x65,
+                                     0x5D,0xA4,0xFB,0xFC,0x0E,0x11,0x08,0xA8,
+                                     0xFD,0x17,0xB4,0x48,0xA6,0x85,0x54,0x19,
+                                     0x9C,0x47,0xD0,0x8F,0xFB,0x10,0xD4,0xB8};
+
+// properties of secp256k1's efficiently computable endomorphism
+static const unsigned char lambda_[] = {0x53,0x63,0xad,0x4c,0xc0,0x5c,0x30,0xe0,
+                                        0xa5,0x26,0x1c,0x02,0x88,0x12,0x64,0x5a,
+                                        0x12,0x2e,0x22,0xea,0x20,0x81,0x66,0x78,
+                                        0xdf,0x02,0x96,0x7c,0x1b,0x23,0xbd,0x72};
+static const unsigned char beta_[] =   {0x7a,0xe9,0x6a,0x2b,0x65,0x7c,0x07,0x10,
+                                        0x6e,0x64,0x47,0x9e,0xac,0x34,0x34,0xe9,
+                                        0x9c,0xf0,0x49,0x75,0x12,0xf5,0x89,0x95,
+                                        0xc1,0x39,0x6c,0x28,0x71,0x95,0x01,0xee};
+static const unsigned char a1b2_[] =   {0x30,0x86,0xd2,0x21,0xa7,0xd4,0x6b,0xcd,
+                                        0xe8,0x6c,0x90,0xe4,0x92,0x84,0xeb,0x15};
+static const unsigned char b1_[]   =   {0xe4,0x43,0x7e,0xd6,0x01,0x0e,0x88,0x28,
+                                        0x6f,0x54,0x7f,0xa9,0x0a,0xbf,0xe4,0xc3};
+static const unsigned char a2_[]   =   {0x01,
+                                        0x14,0xca,0x50,0xf7,0xa8,0xe2,0xf3,0xf6,
+                                        0x57,0xc1,0x10,0x8d,0x9d,0x44,0xcf,0xd8};
+
+GroupConstants::GroupConstants() : g_x(g_x_), g_y(g_y_),
+                                   order(order_, sizeof(order_)),
+                                   g(g_x,g_y),
+                                   beta(beta_),
+                                   lambda(lambda_, sizeof(lambda_)),
+                                   a1b2(a1b2_, sizeof(a1b2_)),
+                                   b1(b1_, sizeof(b1_)),
+                                   a2(a2_, sizeof(a2_)) {}
+
+const GroupConstants &GetGroupConst() {
+    static const GroupConstants group_const;
+    return group_const;
+}
+
+void GroupElemJac::SetMulLambda(const GroupElemJac &p) {
+    FieldElem beta = GetGroupConst().beta;
+    *this = p;
+    x.SetMult(x, beta);
+}
+
+void SplitExp(const Number &exp, Number &exp1, Number &exp2) {
+    const GroupConstants &c = GetGroupConst();
+    Number bnc1, bnc2, bnt1, bnt2, bnn2;
+    bnn2.SetNumber(c.order);
+    bnn2.Shift1();
+
+    bnc1.SetMult(exp, c.a1b2);
+    bnc1.SetAdd(bnc1, bnn2);
+    bnc1.SetDiv(bnc1, c.order);
+
+    bnc2.SetMult(exp, c.b1);
+    bnc2.SetAdd(bnc2, bnn2);
+    bnc2.SetDiv(bnc2, c.order);
+
+    bnt1.SetMult(bnc1, c.a1b2);
+    bnt2.SetMult(bnc2, c.a2);
+    bnt1.SetAdd(bnt1, bnt2);
+    exp1.SetSub(exp, bnt1);
+    bnt1.SetMult(bnc1, c.b1);
+    bnt2.SetMult(bnc2, c.a1b2);
+    exp2.SetSub(bnt1, bnt2);
+}
+
+}
diff --git a/group.h b/group.h
index 3dc233e..d95e08c 100644
--- a/group.h
+++ b/group.h
@@ -1,6 +1,9 @@
 #ifndef _SECP256K1_GROUP_
 #define _SECP256K1_GROUP_
 
+#include <string>
+
+#include "num.h"
 #include "field.h"
 
 namespace secp256k1 {
@@ -17,42 +20,21 @@ protected:
 public:
 
     /** Creates the point at infinity */
-    GroupElem() {
-        fInfinity = true;
-    }
+    GroupElem();
 
     /** Creates the point with given affine coordinates */
-    GroupElem(const FieldElem &xin, const FieldElem &yin) {
-        fInfinity = false;
-        x = xin;
-        y = yin;
-    }
+    GroupElem(const FieldElem &xin, const FieldElem &yin);
 
     /** Checks whether this is the point at infinity */
-    bool IsInfinity() const {
-        return fInfinity;
-    }
+    bool IsInfinity() const;
 
-    void SetNeg(const GroupElem &p) {
-        *this = p;
-        y.Normalize();
-        y.SetNeg(y, 1);
-    }
+    void SetNeg(const GroupElem &p);
 
-    void GetX(FieldElem &xout) {
-        xout = x;
-    }
+    void GetX(FieldElem &xout);
 
-    void GetY(FieldElem &yout) {
-        yout = y;
-    }
+    void GetY(FieldElem &yout);
 
-    std::string ToString() const {
-        if (fInfinity)
-            return "(inf)";
-        FieldElem xc = x, yc = y;
-        return "(" + xc.ToString() + "," + yc.ToString() + ")";
-    }
+    std::string ToString() const;
 
     void SetJac(GroupElemJac &jac);
 
@@ -66,244 +48,45 @@ protected:
 
 public:
     /** Creates the point at infinity */
-    GroupElemJac() : GroupElem(), z(1) {}
+    GroupElemJac();
 
     /** Creates the point with given affine coordinates */
-    GroupElemJac(const FieldElem &xin, const FieldElem &yin) : GroupElem(xin,yin), z(1) {}
+    GroupElemJac(const FieldElem &xin, const FieldElem &yin);
 
-    GroupElemJac(const GroupElem &in) : GroupElem(in), z(1) {}
+    GroupElemJac(const GroupElem &in);
 
-    void SetJac(GroupElemJac &jac) {
-        *this = jac;
-    }
+    void SetJac(GroupElemJac &jac);
 
     /** Checks whether this is a non-infinite point on the curve */
-    bool IsValid() const {
-        if (IsInfinity())
-            return false;
-        // y^2 = x^3 + 7
-        // (Y/Z^3)^2 = (X/Z^2)^3 + 7
-        // Y^2 / Z^6 = X^3 / Z^6 + 7
-        // Y^2 = X^3 + 7*Z^6
-        FieldElem y2; y2.SetSquare(y);
-        FieldElem x3; x3.SetSquare(x); x3.SetMult(x3,x);
-        FieldElem z2; z2.SetSquare(z);
-        FieldElem z6; z6.SetSquare(z2); z6.SetMult(z6,z2);
-        z6 *= 7;
-        x3 += z6;
-        return y2 == x3;
-    }
+    bool IsValid() const;
 
     /** Returns the affine coordinates of this point */
-    void GetAffine(GroupElem &aff) {
-        z.SetInverse(z);
-        FieldElem z2;
-        z2.SetSquare(z);
-        FieldElem z3;
-        z3.SetMult(z,z2);
-        x.SetMult(x,z2);
-        y.SetMult(y,z3);
-        z = FieldElem(1);
-        aff.fInfinity = fInfinity;
-        aff.x = x;
-        aff.y = y;
-    }
+    void GetAffine(GroupElem &aff);
 
-    void GetX(FieldElem &xout) {
-        FieldElem zi;
-        zi.SetInverse(z);
-        zi.SetSquare(zi);
-        xout.SetMult(x, zi);
-    }
+    void GetX(FieldElem &xout);
+    void GetY(FieldElem &yout);
 
-    bool IsInfinity() const {
-        return fInfinity;
-    }
+    bool IsInfinity() const;
 
-    void GetY(FieldElem &yout) {
-        FieldElem zi;
-        zi.SetInverse(z);
-        FieldElem zi3; zi3.SetSquare(zi); zi3.SetMult(zi, zi3);
-        yout.SetMult(y, zi3);
-    }
-
-    void SetNeg(const GroupElemJac &p) {
-        *this = p;
-        y.Normalize();
-        y.SetNeg(y, 1);
-    }
+    void SetNeg(const GroupElemJac &p);
 
     /** Sets this point to have a given X coordinate & given Y oddness */
-    void SetCompressed(const FieldElem &xin, bool fOdd) {
-        x = xin;
-        FieldElem x2; x2.SetSquare(x);
-        FieldElem x3; x3.SetMult(x,x2);
-        fInfinity = false;
-        FieldElem c(7);
-        c += x3;
-        y.SetSquareRoot(c);
-        z = FieldElem(1);
-        if (y.IsOdd() != fOdd)
-            y.SetNeg(y,1);
-    }
+    void SetCompressed(const FieldElem &xin, bool fOdd);
 
     /** Sets this point to be the EC double of another */
-    void SetDouble(const GroupElemJac &p) {
-        FieldElem t5 = p.y;
-        if (p.fInfinity || t5.IsZero()) {
-            fInfinity = true;
-            return;
-        }
-
-        FieldElem t1,t2,t3,t4;
-        z.SetMult(t5,p.z);
-        z *= 2;                // Z' = 2*Y*Z (2)
-        t1.SetSquare(p.x);
-        t1 *= 3;               // T1 = 3*X^2 (3)
-        t2.SetSquare(t1);      // T2 = 9*X^4 (1)
-        t3.SetSquare(t5);
-        t3 *= 2;               // T3 = 2*Y^2 (2)
-        t4.SetSquare(t3);
-        t4 *= 2;               // T4 = 8*Y^4 (2)
-        t3.SetMult(p.x,t3);      // T3 = 2*X*Y^2 (1)
-        x = t3;
-        x *= 4;                // X' = 8*X*Y^2 (4)
-        x.SetNeg(x,4);         // X' = -8*X*Y^2 (5)
-        x += t2;               // X' = 9*X^4 - 8*X*Y^2 (6)
-        t2.SetNeg(t2,1);       // T2 = -9*X^4 (2)
-        t3 *= 6;               // T3 = 12*X*Y^2 (6)
-        t3 += t2;              // T3 = 12*X*Y^2 - 9*X^4 (8)
-        y.SetMult(t1,t3);      // Y' = 36*X^3*Y^2 - 27*X^6 (1)
-        t2.SetNeg(t4,2);       // T2 = -8*Y^4 (3)
-        y += t2;               // Y' = 36*X^3*Y^2 - 27*X^6 - 8*Y^4 (4)
-        fInfinity = false;
-    }
+    void SetDouble(const GroupElemJac &p);
 
     /** Sets this point to be the EC addition of two others */
-    void SetAdd(const GroupElemJac &p, const GroupElemJac &q) {
-        if (p.fInfinity) {
-            *this = q;
-            return;
-        }
-        if (q.fInfinity) {
-            *this = p;
-            return;
-        }
-        fInfinity = false;
-        const FieldElem &x1 = p.x, &y1 = p.y, &z1 = p.z, &x2 = q.x, &y2 = q.y, &z2 = q.z;
-        FieldElem z22; z22.SetSquare(z2);
-        FieldElem z12; z12.SetSquare(z1);
-        FieldElem u1; u1.SetMult(x1, z22);
-        FieldElem u2; u2.SetMult(x2, z12);
-        FieldElem s1; s1.SetMult(y1, z22); s1.SetMult(s1, z2);
-        FieldElem s2; s2.SetMult(y2, z12); s2.SetMult(s2, z1);
-        if (u1 == u2) {
-            if (s1 == s2) {
-                SetDouble(p);
-            } else {
-                fInfinity = true;
-            }
-            return;
-        }
-        FieldElem h; h.SetNeg(u1,1); h += u2;
-        FieldElem r; r.SetNeg(s1,1); r += s2;
-        FieldElem r2; r2.SetSquare(r);
-        FieldElem h2; h2.SetSquare(h);
-        FieldElem h3; h3.SetMult(h,h2);
-        z.SetMult(z1,z2); z.SetMult(z, h);
-        FieldElem t; t.SetMult(u1,h2);
-        x = t; x *= 2; x += h3; x.SetNeg(x,3); x += r2;
-        y.SetNeg(x,5); y += t; y.SetMult(y,r);
-        h3.SetMult(h3,s1); h3.SetNeg(h3,1);
-        y += h3;
-    }
+    void SetAdd(const GroupElemJac &p, const GroupElemJac &q);
 
     /** Sets this point to be the EC addition of two others (one of which is in affine coordinates) */
-    void SetAdd(const GroupElemJac &p, const GroupElem &q) {
-        if (p.fInfinity) {
-            x = q.x;
-            y = q.y;
-            fInfinity = q.fInfinity;
-            z = FieldElem(1);
-            return;
-        }
-        if (q.fInfinity) {
-            *this = p;
-            return;
-        }
-        fInfinity = false;
-        const FieldElem &x1 = p.x, &y1 = p.y, &z1 = p.z, &x2 = q.x, &y2 = q.y;
-        FieldElem z12; z12.SetSquare(z1);
-        FieldElem u1 = x1; u1.Normalize();
-        FieldElem u2; u2.SetMult(x2, z12);
-        FieldElem s1 = y1; s1.Normalize();
-        FieldElem s2; s2.SetMult(y2, z12); s2.SetMult(s2, z1);
-        if (u1 == u2) {
-            if (s1 == s2) {
-                SetDouble(p);
-            } else {
-                fInfinity = true;
-            }
-            return;
-        }
-        FieldElem h; h.SetNeg(u1,1); h += u2;
-        FieldElem r; r.SetNeg(s1,1); r += s2;
-        FieldElem r2; r2.SetSquare(r);
-        FieldElem h2; h2.SetSquare(h);
-        FieldElem h3; h3.SetMult(h,h2);
-        z = p.z; z.SetMult(z, h);
-        FieldElem t; t.SetMult(u1,h2);
-        x = t; x *= 2; x += h3; x.SetNeg(x,3); x += r2;
-        y.SetNeg(x,5); y += t; y.SetMult(y,r);
-        h3.SetMult(h3,s1); h3.SetNeg(h3,1);
-        y += h3;
-    }
+    void SetAdd(const GroupElemJac &p, const GroupElem &q);
 
-    std::string ToString() const {
-        GroupElemJac cop = *this;
-        GroupElem aff;
-        cop.GetAffine(aff);
-        return aff.ToString();
-    }
+    std::string ToString() const;
 
     void SetMulLambda(const GroupElemJac &p);
 };
 
-void GroupElem::SetJac(GroupElemJac &jac) {
-    jac.GetAffine(*this);
-}
-
-static const unsigned char order_[] = {0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,
-                                       0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFE,
-                                       0xBA,0xAE,0xDC,0xE6,0xAF,0x48,0xA0,0x3B,
-                                       0xBF,0xD2,0x5E,0x8C,0xD0,0x36,0x41,0x41};
-
-static const unsigned char g_x_[] = {0x79,0xBE,0x66,0x7E,0xF9,0xDC,0xBB,0xAC,
-                                     0x55,0xA0,0x62,0x95,0xCE,0x87,0x0B,0x07,
-                                     0x02,0x9B,0xFC,0xDB,0x2D,0xCE,0x28,0xD9,
-                                     0x59,0xF2,0x81,0x5B,0x16,0xF8,0x17,0x98};
-
-static const unsigned char g_y_[] = {0x48,0x3A,0xDA,0x77,0x26,0xA3,0xC4,0x65,
-                                     0x5D,0xA4,0xFB,0xFC,0x0E,0x11,0x08,0xA8,
-                                     0xFD,0x17,0xB4,0x48,0xA6,0x85,0x54,0x19,
-                                     0x9C,0x47,0xD0,0x8F,0xFB,0x10,0xD4,0xB8};
-
-// properties of secp256k1's efficiently computable endomorphism
-static const unsigned char lambda_[] = {0x53,0x63,0xad,0x4c,0xc0,0x5c,0x30,0xe0,
-                                        0xa5,0x26,0x1c,0x02,0x88,0x12,0x64,0x5a,
-                                        0x12,0x2e,0x22,0xea,0x20,0x81,0x66,0x78,
-                                        0xdf,0x02,0x96,0x7c,0x1b,0x23,0xbd,0x72};
-static const unsigned char beta_[] =   {0x7a,0xe9,0x6a,0x2b,0x65,0x7c,0x07,0x10,
-                                        0x6e,0x64,0x47,0x9e,0xac,0x34,0x34,0xe9,
-                                        0x9c,0xf0,0x49,0x75,0x12,0xf5,0x89,0x95,
-                                        0xc1,0x39,0x6c,0x28,0x71,0x95,0x01,0xee};
-static const unsigned char a1b2_[] =   {0x30,0x86,0xd2,0x21,0xa7,0xd4,0x6b,0xcd,
-                                        0xe8,0x6c,0x90,0xe4,0x92,0x84,0xeb,0x15};
-static const unsigned char b1_[]   =   {0xe4,0x43,0x7e,0xd6,0x01,0x0e,0x88,0x28,
-                                        0x6f,0x54,0x7f,0xa9,0x0a,0xbf,0xe4,0xc3};
-static const unsigned char a2_[]   =   {0x01,
-                                        0x14,0xca,0x50,0xf7,0xa8,0xe2,0xf3,0xf6,
-                                        0x57,0xc1,0x10,0x8d,0x9d,0x44,0xcf,0xd8};
 class GroupConstants {
 private:
     const FieldElem g_x;
@@ -315,49 +98,12 @@ public:
     const FieldElem beta;
     const Number lambda, a1b2, b1, a2;
 
-    GroupConstants() : g_x(g_x_), g_y(g_y_),
-                       order(order_, sizeof(order_)),
-                       g(g_x,g_y),
-                       beta(beta_),
-                       lambda(lambda_, sizeof(lambda_)),
-                       a1b2(a1b2_, sizeof(a1b2_)),
-                       b1(b1_, sizeof(b1_)),
-                       a2(a2_, sizeof(a2_)) {}
+    GroupConstants();
 };
 
-const GroupConstants &GetGroupConst() {
-    static const GroupConstants group_const;
-    return group_const;
-}
+const GroupConstants &GetGroupConst();
 
-void GroupElemJac::SetMulLambda(const GroupElemJac &p) {
-    FieldElem beta = GetGroupConst().beta;
-    *this = p;
-    x.SetMult(x, beta);
-}
-
-void SplitExp(const Number &exp, Number &exp1, Number &exp2) {
-    const GroupConstants &c = GetGroupConst();
-    Number bnc1, bnc2, bnt1, bnt2, bnn2;
-    bnn2.SetNumber(c.order);
-    bnn2.Shift1();
-
-    bnc1.SetMult(exp, c.a1b2);
-    bnc1.SetAdd(bnc1, bnn2);
-    bnc1.SetDiv(bnc1, c.order);
-
-    bnc2.SetMult(exp, c.b1);
-    bnc2.SetAdd(bnc2, bnn2);
-    bnc2.SetDiv(bnc2, c.order);
-
-    bnt1.SetMult(bnc1, c.a1b2);
-    bnt2.SetMult(bnc2, c.a2);
-    bnt1.SetAdd(bnt1, bnt2);
-    exp1.SetSub(exp, bnt1);
-    bnt1.SetMult(bnc1, c.b1);
-    bnt2.SetMult(bnc2, c.a1b2);
-    exp2.SetSub(bnt1, bnt2);
-}
+void SplitExp(const Number &exp, Number &exp1, Number &exp2);
 
 }
 
diff --git a/num.cpp b/num.cpp
new file mode 100644
index 0000000..02afeb8
--- /dev/null
+++ b/num.cpp
@@ -0,0 +1,7 @@
+#if defined(USE_NUM_GMP)
+#include "num_gmp.cpp"
+#elif defined(USE_NUM_OPENSSL)
+#include "num_openssl.cpp"
+#else
+#error "Please select num implementation"
+#endif
diff --git a/num_builtin.h b/num_builtin.h
new file mode 100644
index 0000000..a77e5ff
--- /dev/null
+++ b/num_builtin.h
@@ -0,0 +1,155 @@
+#ifndef _SECP256K1_NUM_GMP_
+#define _SECP256K1_NUM_GMP_
+
+#include <assert.h>
+#include <string>
+#include <string.h>
+#include <stdlib.h>
+
+namespace secp256k1 {
+
+class Number {
+private:
+    uint64_t n[8]; // 512 bits ought to be enough for everyone
+    int top; // number of used entries in n
+
+    void FixTop() {
+        while (top > 0 && n[top-1] == 0)
+            top--;
+    }
+
+    int GetNumBytes() {
+        if (top==0)
+            return 0;
+        int ret = 8*(top-1);
+        uint64_t h = n[top-1];
+        while (h>0) {
+            ret++;
+            h >>= 8;
+        }
+        return ret;
+    }
+
+    Number(const Number &c) {}
+public:
+    Number() {
+        top = 0;
+    }
+    Number(const unsigned char *bin, int len) {
+        top = 0;
+        SetBytes(bin, len);
+    }
+    Number &operator=(const Number &x) {
+        for (int pos = 0; pos < x.top; pos++)
+            n[pos] = x.n[pos];
+        top = x.top;
+    }
+    void SetNumber(const Number &x) {
+        *this = x;
+    }
+    void SetBytes(const unsigned char *bin, unsigned int len) {
+        n[0] = n[1] = n[2] = n[3] = n[4] = n[5] = n[6] = n[7] = 0;
+        assert(len<=64);
+        for (int pos = 0; pos < len; pos++)
+            n[(len-1-pos)/8] |= ((uint64_t)(bin[pos])) << (((len-1-pos)%8)*8);
+        top = 8;
+        FixTop();
+    }
+
+    void GetBytes(unsigned char *bin, unsigned int len) {
+        unsigned int size = GetNumBytes();
+        assert(size <= len);
+        memset(bin,0,len);
+        for (int i=0; i<size; i++)
+            bin[i] = (n[(size-i-1)/8] >> (((size-1-i)%8)*8)) & 0xFF;
+    }
+
+    void SetInt(int x) {
+        n[0] = x;
+        top = 1;
+        FixTop();
+    }
+
+    void SetModInverse(const Number &x, const Number &m) {
+        mpz_invert(bn, x.bn, m.bn);
+    }
+
+    void SetModMul(const Number &a, const Number &b, const Number &m) {
+        mpz_mul(bn, a.bn, b.bn);
+        mpz_mod(bn, bn, m.bn);
+    }
+    void SetAdd(const Number &a1, const Number &a2) {
+        
+        mpz_add(bn, a1.bn, a2.bn);
+    }
+    void SetSub(const Number &a1, const Number &a2) {
+        mpz_sub(bn, a1.bn, a2.bn);
+    }
+    void SetMult(const Number &a1, const Number &a2) {
+        mpz_mul(bn, a1.bn, a2.bn);
+    }
+    void SetDiv(const Number &a1, const Number &a2) {
+        mpz_tdiv_q(bn, a1.bn, a2.bn);
+    }
+    void SetMod(const Number &a, const Number &m) {
+        mpz_mod(bn, a.bn, m.bn);
+    }
+    int Compare(const Number &a) const {
+        return mpz_cmp(bn, a.bn);
+    }
+    int GetBits() const {
+        return mpz_sizeinbase(bn,2);
+    }
+    // return the lowest (rightmost) bits bits, and rshift them away
+    int ShiftLowBits(int bits) {
+        int ret = mpz_get_ui(bn) & ((1 << bits) - 1);
+        mpz_fdiv_q_2exp(bn, bn, bits);
+        return ret;
+    }
+    // check whether number is 0,
+    bool IsZero() const {
+        return mpz_size(bn) == 0;
+    }
+    bool IsOdd() const {
+        return mpz_get_ui(bn) & 1;
+    }
+    bool IsNeg() const {
+        return mpz_sgn(bn) < 0;
+    }
+    void Negate() {
+        mpz_neg(bn, bn);
+    }
+    void Shift1() {
+        mpz_fdiv_q_2exp(bn, bn, 1);
+    }
+    void Inc() {
+        mpz_add_ui(bn, bn, 1);
+    }
+    void SetHex(const std::string &str) {
+        mpz_set_str(bn, str.c_str(), 16);
+    }
+    void SetPseudoRand(const Number &max) {
+        number_state.gen(bn, max.bn);
+    }
+    void SplitInto(int bits, Number &low, Number &high) const {
+        mpz_t tmp;
+        mpz_init_set_ui(tmp,1);
+        mpz_mul_2exp(tmp,tmp,bits);
+        mpz_sub_ui(tmp,tmp,1);
+        mpz_and(low.bn, bn, tmp);
+        mpz_clear(tmp);
+        mpz_fdiv_q_2exp(high.bn, bn, bits);
+    }
+
+    std::string ToString() const {
+        char *str = (char*)malloc(mpz_sizeinbase(bn,16) + 2);
+        mpz_get_str(str, 16, bn);
+        std::string ret(str);
+        free(str);
+        return ret;
+    }
+};
+
+}
+
+#endif
diff --git a/num_gmp.cpp b/num_gmp.cpp
new file mode 100644
index 0000000..7945d60
--- /dev/null
+++ b/num_gmp.cpp
@@ -0,0 +1,167 @@
+#include <assert.h>
+#include <string>
+#include <string.h>
+#include <stdlib.h>
+#include <gmp.h>
+
+#include "num_gmp.h"
+
+namespace secp256k1 {
+
+class NumberState {
+private:
+    gmp_randstate_t rng;
+
+public:
+    NumberState() {
+        gmp_randinit_default(rng);
+    }
+
+    ~NumberState() {
+        gmp_randclear(rng);
+    }
+
+    void gen(mpz_t out, mpz_t size) {
+        mpz_urandomm(out, rng, size);
+    }
+};
+
+static NumberState number_state;
+
+Number::Number(const Number &x) {
+    mpz_init_set(bn, x.bn);
+}
+
+Number::Number() {
+    mpz_init(bn);
+}
+
+Number::~Number() {
+    mpz_clear(bn);
+}
+
+Number &Number::operator=(const Number &x) {
+    mpz_set(bn, x.bn);
+    return *this;
+}
+
+void Number::SetNumber(const Number &x) {
+    mpz_set(bn, x.bn);
+}
+
+Number::Number(const unsigned char *bin, int len) {
+    mpz_init(bn);
+    SetBytes(bin,len);
+}
+
+void Number::SetBytes(const unsigned char *bin, unsigned int len) {
+    mpz_import(bn, len, 1, 1, 1, 0, bin);
+}
+
+void Number::GetBytes(unsigned char *bin, unsigned int len) {
+    unsigned int size = (mpz_sizeinbase(bn,2)+7)/8;
+    assert(size <= len);
+    memset(bin,0,len);
+    size_t count = 0;
+    mpz_export(bin + len - size, &count, 1, 1, 1, 0, bn);
+    assert(count == 0 || size == count);
+}
+
+void Number::SetInt(int x) {
+    mpz_set_si(bn, x);
+}
+
+void Number::SetModInverse(const Number &x, const Number &m) {
+    mpz_invert(bn, x.bn, m.bn);
+}
+
+void Number::SetModMul(const Number &a, const Number &b, const Number &m) {
+    mpz_mul(bn, a.bn, b.bn);
+    mpz_mod(bn, bn, m.bn);
+}
+
+void Number::SetAdd(const Number &a1, const Number &a2) {
+    mpz_add(bn, a1.bn, a2.bn);
+}
+
+void Number::SetSub(const Number &a1, const Number &a2) {
+    mpz_sub(bn, a1.bn, a2.bn);
+}
+
+void Number::SetMult(const Number &a1, const Number &a2) {
+    mpz_mul(bn, a1.bn, a2.bn);
+}
+
+void Number::SetDiv(const Number &a1, const Number &a2) {
+    mpz_tdiv_q(bn, a1.bn, a2.bn);
+}
+
+void Number::SetMod(const Number &a, const Number &m) {
+    mpz_mod(bn, a.bn, m.bn);
+}
+
+int Number::Compare(const Number &a) const {
+    return mpz_cmp(bn, a.bn);
+}
+
+int Number::GetBits() const {
+    return mpz_sizeinbase(bn,2);
+}
+
+int Number::ShiftLowBits(int bits) {
+    int ret = mpz_get_ui(bn) & ((1 << bits) - 1);
+    mpz_fdiv_q_2exp(bn, bn, bits);
+    return ret;
+}
+
+bool Number::IsZero() const {
+    return mpz_size(bn) == 0;
+}
+
+bool Number::IsOdd() const {
+    return mpz_get_ui(bn) & 1;
+}
+
+bool Number::IsNeg() const {
+    return mpz_sgn(bn) < 0;
+}
+
+void Number::Negate() {
+    mpz_neg(bn, bn);
+}
+
+void Number::Shift1() {
+    mpz_fdiv_q_2exp(bn, bn, 1);
+}
+
+void Number::Inc() {
+    mpz_add_ui(bn, bn, 1);
+}
+
+void Number::SetHex(const std::string &str) {
+    mpz_set_str(bn, str.c_str(), 16);
+}
+
+void Number::SetPseudoRand(const Number &max) {
+    number_state.gen(bn, max.bn);
+}
+
+void Number::SplitInto(int bits, Number &low, Number &high) const {
+    mpz_t tmp;
+    mpz_init_set_ui(tmp,1);
+    mpz_mul_2exp(tmp,tmp,bits);
+    mpz_sub_ui(tmp,tmp,1);
+    mpz_and(low.bn, bn, tmp);
+    mpz_clear(tmp);
+    mpz_fdiv_q_2exp(high.bn, bn, bits);
+}
+
+std::string Number::ToString() const {
+    char *str = (char*)malloc(mpz_sizeinbase(bn,16) + 2);
+    mpz_get_str(str, 16, bn);
+    std::string ret(str);
+    free(str);
+    return ret;
+}
+
+}
diff --git a/num_gmp.h b/num_gmp.h
index 954e3ff..ea330b4 100644
--- a/num_gmp.h
+++ b/num_gmp.h
@@ -1,147 +1,45 @@
 #ifndef _SECP256K1_NUM_GMP_
 #define _SECP256K1_NUM_GMP_
 
-#include <assert.h>
 #include <string>
-#include <string.h>
-#include <stdlib.h>
 #include <gmp.h>
 
 namespace secp256k1 {
 
-class NumberState {
-private:
-    gmp_randstate_t rng;
-
-public:
-    NumberState() {
-        gmp_randinit_default(rng);
-    }
-
-    ~NumberState() {
-        gmp_randclear(rng);
-    }
-
-    void gen(mpz_t out, mpz_t size) {
-        mpz_urandomm(out, rng, size);
-    }
-};
-
-static NumberState number_state;
-
 class Number {
 private:
     mutable mpz_t bn;
-    Number(const Number &x) {
-    }
-public:
-    Number() {
-        mpz_init(bn);
-    }
-    ~Number() {
-        mpz_clear(bn);
-    }
-    Number(const unsigned char *bin, int len) {
-        mpz_init(bn);
-        SetBytes(bin,len);
-    }
-    Number &operator=(const Number &x) {
-        mpz_set(bn, x.bn);
-        return *this;
-    }
-    void SetNumber(const Number &x) {
-        mpz_set(bn, x.bn);
-    }
-    void SetBytes(const unsigned char *bin, unsigned int len) {
-        mpz_import(bn, len, 1, 1, 1, 0, bin);
-    }
-    void GetBytes(unsigned char *bin, unsigned int len) {
-        unsigned int size = (mpz_sizeinbase(bn,2)+7)/8;
-        assert(size <= len);
-        memset(bin,0,len);
-        size_t count = 0;
-        mpz_export(bin + len - size, &count, 1, 1, 1, 0, bn);
-        assert(count == 0 || size == count);
-    }
-    void SetInt(int x) {
-        mpz_set_si(bn, x);
-    }
-    void SetModInverse(const Number &x, const Number &m) {
-        mpz_invert(bn, x.bn, m.bn);
-    }
-    void SetModMul(const Number &a, const Number &b, const Number &m) {
-        mpz_mul(bn, a.bn, b.bn);
-        mpz_mod(bn, bn, m.bn);
-    }
-    void SetAdd(const Number &a1, const Number &a2) {
-        mpz_add(bn, a1.bn, a2.bn);
-    }
-    void SetSub(const Number &a1, const Number &a2) {
-        mpz_sub(bn, a1.bn, a2.bn);
-    }
-    void SetMult(const Number &a1, const Number &a2) {
-        mpz_mul(bn, a1.bn, a2.bn);
-    }
-    void SetDiv(const Number &a1, const Number &a2) {
-        mpz_tdiv_q(bn, a1.bn, a2.bn);
-    }
-    void SetMod(const Number &a, const Number &m) {
-        mpz_mod(bn, a.bn, m.bn);
-    }
-    int Compare(const Number &a) const {
-        return mpz_cmp(bn, a.bn);
-    }
-    int GetBits() const {
-        return mpz_sizeinbase(bn,2);
-    }
-    // return the lowest (rightmost) bits bits, and rshift them away
-    int ShiftLowBits(int bits) {
-        int ret = mpz_get_ui(bn) & ((1 << bits) - 1);
-        mpz_fdiv_q_2exp(bn, bn, bits);
-        return ret;
-    }
-    // check whether number is 0,
-    bool IsZero() const {
-        return mpz_size(bn) == 0;
-    }
-    bool IsOdd() const {
-        return mpz_get_ui(bn) & 1;
-    }
-    bool IsNeg() const {
-        return mpz_sgn(bn) < 0;
-    }
-    void Negate() {
-        mpz_neg(bn, bn);
-    }
-    void Shift1() {
-        mpz_fdiv_q_2exp(bn, bn, 1);
-    }
-    void Inc() {
-        mpz_add_ui(bn, bn, 1);
-    }
-    void SetHex(const std::string &str) {
-        mpz_set_str(bn, str.c_str(), 16);
-    }
-    void SetPseudoRand(const Number &max) {
-        number_state.gen(bn, max.bn);
-    }
-    void SplitInto(int bits, Number &low, Number &high) const {
-        mpz_t tmp;
-        mpz_init_set_ui(tmp,1);
-        mpz_mul_2exp(tmp,tmp,bits);
-        mpz_sub_ui(tmp,tmp,1);
-        mpz_and(low.bn, bn, tmp);
-        mpz_clear(tmp);
-        mpz_fdiv_q_2exp(high.bn, bn, bits);
-    }
+    Number(const Number &x);
 
-    std::string ToString() const {
-        char *str = (char*)malloc(mpz_sizeinbase(bn,16) + 2);
-        mpz_get_str(str, 16, bn);
-        std::string ret(str);
-        free(str);
-        return ret;
-    }
+public:
+    Number();
+    ~Number();
+    Number(const unsigned char *bin, int len);
+    Number &operator=(const Number &x);
+    void SetNumber(const Number &x);
+    void SetBytes(const unsigned char *bin, unsigned int len);
+    void GetBytes(unsigned char *bin, unsigned int len);
+    void SetInt(int x);
+    void SetModInverse(const Number &x, const Number &m);
+    void SetModMul(const Number &a, const Number &b, const Number &m);
+    void SetAdd(const Number &a1, const Number &a2);
+    void SetSub(const Number &a1, const Number &a2);
+    void SetMult(const Number &a1, const Number &a2);
+    void SetDiv(const Number &a1, const Number &a2);
+    void SetMod(const Number &a, const Number &m);
+    int Compare(const Number &a) const;
+    int GetBits() const;
+    int ShiftLowBits(int bits);
+    bool IsZero() const;
+    bool IsOdd() const;
+    bool IsNeg() const;
+    void Negate();
+    void Shift1();
+    void Inc();
+    void SetHex(const std::string &str);
+    void SetPseudoRand(const Number &max);
+    void SplitInto(int bits, Number &low, Number &high) const;
+    std::string ToString() const;
 };
 
 }
diff --git a/num_openssl.cpp b/num_openssl.cpp
new file mode 100644
index 0000000..c95a733
--- /dev/null
+++ b/num_openssl.cpp
@@ -0,0 +1,179 @@
+#include <assert.h>
+#include <string>
+#include <string.h>
+#include <openssl/bn.h>
+#include <openssl/crypto.h>
+
+#include "num_openssl.h"
+
+namespace secp256k1 {
+
+class Context {
+private:
+    BN_CTX *ctx;
+
+    operator BN_CTX*() {
+        return ctx;
+    }
+
+    friend class Number;
+public:
+    Context() {
+        ctx = BN_CTX_new();
+    }
+
+    ~Context() {
+        BN_CTX_free(ctx);
+    }
+};
+
+Number::operator const BIGNUM*() const {
+    return &b;
+}
+
+Number::operator BIGNUM*() {
+    return &b;
+}
+
+Number::Number() {
+    BN_init(*this);
+}
+
+Number::~Number() {
+    BN_free(*this);
+}
+
+Number::Number(const unsigned char *bin, int len) {
+    BN_init(*this);
+    SetBytes(bin,len);
+}
+
+void Number::SetNumber(const Number &x) {
+    BN_copy(*this, x);
+}
+
+Number::Number(const Number &x) {
+    BN_init(*this);
+    BN_copy(*this, x);
+}
+
+Number &Number::operator=(const Number &x) {
+    BN_copy(*this, x);
+    return *this;
+}
+
+void Number::SetBytes(const unsigned char *bin, int len) {
+    BN_bin2bn(bin, len, *this);
+}
+
+void Number::GetBytes(unsigned char *bin, int len) {
+    int size = BN_num_bytes(*this);
+    assert(size <= len);
+    memset(bin,0,len);
+    BN_bn2bin(*this, bin + len - size);
+}
+
+void Number::SetInt(int x) {
+    if (x >= 0) {
+        BN_set_word(*this, x);
+    } else {
+        BN_set_word(*this, -x);
+        BN_set_negative(*this, 1);
+    }
+}
+
+void Number::SetModInverse(const Number &x, const Number &m) {
+    Context ctx;
+    BN_mod_inverse(*this, x, m, ctx);
+}
+
+void Number::SetModMul(const Number &a, const Number &b, const Number &m) {
+    Context ctx;
+    BN_mod_mul(*this, a, b, m, ctx);
+}
+
+void Number::SetAdd(const Number &a1, const Number &a2) {
+    BN_add(*this, a1, a2);
+}
+
+void Number::SetSub(const Number &a1, const Number &a2) {
+    BN_sub(*this, a1, a2);
+}
+
+void Number::SetMult(const Number &a1, const Number &a2) {
+    Context ctx;
+    BN_mul(*this, a1, a2, ctx);
+}
+
+void Number::SetDiv(const Number &a1, const Number &a2) {
+    Context ctx;
+    BN_div(*this, NULL, a1, a2, ctx);
+}
+
+void Number::SetMod(const Number &a, const Number &m) {
+    Context ctx;
+    BN_nnmod(*this, a, m, ctx);
+}
+
+int Number::Compare(const Number &a) const {
+    return BN_cmp(*this, a);
+}
+
+int Number::GetBits() const {
+    return BN_num_bits(*this);
+}
+
+int Number::ShiftLowBits(int bits) {
+    BIGNUM *bn = *this;
+    int ret = BN_is_zero(bn) ? 0 : bn->d[0] & ((1 << bits) - 1);
+    BN_rshift(*this, *this, bits);
+    return ret;
+}
+
+bool Number::IsZero() const {
+    return BN_is_zero((const BIGNUM*)*this);
+}
+
+bool Number::IsOdd() const {
+    return BN_is_odd((const BIGNUM*)*this);
+}
+
+bool Number::IsNeg() const {
+    return BN_is_negative((const BIGNUM*)*this);
+}
+
+void Number::Negate() {
+    BN_set_negative(*this, !IsNeg());
+}
+
+void Number::Shift1() {
+    BN_rshift1(*this,*this);
+}
+
+void Number::Inc() {
+    BN_add_word(*this,1);
+}
+
+void Number::SetHex(const std::string &str) {
+    BIGNUM *bn = *this;
+    BN_hex2bn(&bn, str.c_str());
+}
+
+void Number::SetPseudoRand(const Number &max) {
+    BN_pseudo_rand_range(*this, max);
+}
+
+void Number::SplitInto(int bits, Number &low, Number &high) const {
+    BN_copy(low, *this);
+    BN_mask_bits(low, bits);
+    BN_rshift(high, *this, bits);
+}
+
+std::string Number::ToString() const {
+    char *str = BN_bn2hex(*this);
+    std::string ret(str);
+    OPENSSL_free(str);
+    return ret;
+}
+
+}
diff --git a/num_openssl.h b/num_openssl.h
index acb92ed..ddd73b9 100644
--- a/num_openssl.h
+++ b/num_openssl.h
@@ -1,161 +1,47 @@
 #ifndef _SECP256K1_NUM_OPENSSL_
 #define _SECP256K1_NUM_OPENSSL_
 
-#include <assert.h>
 #include <string>
-#include <string.h>
 #include <openssl/bn.h>
-#include <openssl/crypto.h>
-
-#include <stdint.h>
 
 namespace secp256k1 {
 
-class Context {
-private:
-    BN_CTX *ctx;
-
-    operator BN_CTX*() {
-        return ctx;
-    }
-
-    friend class Number;
-public:
-    Context() {
-        ctx = BN_CTX_new();
-    }
-
-    ~Context() {
-        BN_CTX_free(ctx);
-    }
-};
-
 class Number {
 private:
     BIGNUM b;
-    Number(const Number &x) {}
+    Number(const Number &x);
 
-    operator const BIGNUM*() const {
-        return &b;
-    }
-
-    operator BIGNUM*() {
-        return &b;
-    }
+    operator const BIGNUM*() const;
+    operator BIGNUM*();
 public:
-    Number() {
-        BN_init(*this);
-    }
-
-    ~Number() {
-        BN_free(*this);
-    }
-
-    Number(const unsigned char *bin, int len) {
-        BN_init(*this);
-        SetBytes(bin,len);
-    }
-    void SetNumber(const Number &x) {
-        BN_copy(*this, x);
-    }
-    Number &operator=(const Number &x) {
-        BN_copy(*this, x);
-        return *this;
-    }
-    void SetBytes(const unsigned char *bin, int len) {
-        BN_bin2bn(bin, len, *this);
-    }
-    void GetBytes(unsigned char *bin, int len) {
-        int size = BN_num_bytes(*this);
-        assert(size <= len);
-        memset(bin,0,len);
-        BN_bn2bin(*this, bin + len - size);
-    }
-    void SetInt(int x) {
-        if (x >= 0) {
-            BN_set_word(*this, x);
-        } else {
-            BN_set_word(*this, -x);
-            BN_set_negative(*this, 1);
-        }
-    }
-    void SetModInverse(const Number &x, const Number &m) {
-        Context ctx;
-        BN_mod_inverse(*this, x, m, ctx);
-    }
-    void SetModMul(const Number &a, const Number &b, const Number &m) {
-        Context ctx;
-        BN_mod_mul(*this, a, b, m, ctx);
-    }
-    void SetAdd(const Number &a1, const Number &a2) {
-        BN_add(*this, a1, a2);
-    }
-    void SetSub(const Number &a1, const Number &a2) {
-        BN_sub(*this, a1, a2);
-    }
-    void SetMult(const Number &a1, const Number &a2) {
-        Context ctx;
-        BN_mul(*this, a1, a2, ctx);
-    }
-    void SetDiv(const Number &a1, const Number &a2) {
-        Context ctx;
-        BN_div(*this, NULL, a1, a2, ctx);
-    }
-    void SetMod(const Number &a, const Number &m) {
-        Context ctx;
-        BN_nnmod(*this, a, m, ctx);
-    }
-    int Compare(const Number &a) const {
-        return BN_cmp(*this, a);
-    }
-    int GetBits() const {
-        return BN_num_bits(*this);
-    }
-    // return the lowest (rightmost) bits bits, and rshift them away
-    int ShiftLowBits(int bits) {
-        BIGNUM *bn = *this;
-        int ret = BN_is_zero(bn) ? 0 : bn->d[0] & ((1 << bits) - 1);
-        BN_rshift(*this, *this, bits);
-        return ret;
-    }
-    // check whether number is 0,
-    bool IsZero() const {
-        return BN_is_zero((const BIGNUM*)*this);
-    }
-    bool IsOdd() const {
-        return BN_is_odd((const BIGNUM*)*this);
-    }
-    bool IsNeg() const {
-        return BN_is_negative((const BIGNUM*)*this);
-    }
-    void Negate() {
-        BN_set_negative(*this, !IsNeg());
-    }
-    void Shift1() {
-        BN_rshift1(*this,*this);
-    }
-    void Inc() {
-        BN_add_word(*this,1);
-    }
-    void SetHex(const std::string &str) {
-        BIGNUM *bn = *this;
-        BN_hex2bn(&bn, str.c_str());
-    }
-    void SetPseudoRand(const Number &max) {
-        BN_pseudo_rand_range(*this, max);
-    }
-    void SplitInto(int bits, Number &low, Number &high) const {
-        BN_copy(low, *this);
-        BN_mask_bits(low, bits);
-        BN_rshift(high, *this, bits);
-    }
-
-    std::string ToString() const {
-        char *str = BN_bn2hex(*this);
-        std::string ret(str);
-        OPENSSL_free(str);
-        return ret;
-    }
+    Number();
+    ~Number();
+    Number(const unsigned char *bin, int len);
+    void SetNumber(const Number &x);
+    Number &operator=(const Number &x);
+    void SetBytes(const unsigned char *bin, int len);
+    void GetBytes(unsigned char *bin, int len);
+    void SetInt(int x);
+    void SetModInverse(const Number &x, const Number &m);
+    void SetModMul(const Number &a, const Number &b, const Number &m);
+    void SetAdd(const Number &a1, const Number &a2);
+    void SetSub(const Number &a1, const Number &a2);
+    void SetMult(const Number &a1, const Number &a2);
+    void SetDiv(const Number &a1, const Number &a2);
+    void SetMod(const Number &a, const Number &m);
+    int Compare(const Number &a) const;
+    int GetBits() const;
+    int ShiftLowBits(int bits);
+    bool IsZero() const;
+    bool IsOdd() const;
+    bool IsNeg() const;
+    void Negate();
+    void Shift1();
+    void Inc();
+    void SetHex(const std::string &str);
+    void SetPseudoRand(const Number &max);
+    void SplitInto(int bits, Number &low, Number &high) const;
+    std::string ToString() const;
 };
 
 }
diff --git a/secp256k1.cpp b/secp256k1.cpp
index dbf83cc..0ea273a 100644
--- a/secp256k1.cpp
+++ b/secp256k1.cpp
@@ -1,5 +1,28 @@
-#include "num.h"
-#include "field.h"
-#include "group.h"
-#include "ecmult.h"
-#include "ecdsa.h"
+#include "num.cpp"
+#include "field.cpp"
+#include "group.cpp"
+#include "ecmult.cpp"
+#include "ecdsa.cpp"
+
+namespace secp256k1 {
+
+int VerifyECDSA(const unsigned char *msg, int msglen, const unsigned char *sig, int siglen, const unsigned char *pubkey, int pubkeylen) {
+    Number m; 
+    Signature s;
+    GroupElemJac q;
+    m.SetBytes(msg, msglen);
+    if (!ParsePubKey(q, pubkey, pubkeylen))
+        return -1;
+    if (!s.Parse(sig, siglen)) {
+        fprintf(stderr, "Can't parse signature: ");
+        for (int i=0; i<siglen; i++) fprintf(stderr,"%02x", sig[i]);
+        fprintf(stderr, "\n");
+        return -2;
+    }
+//    fprintf(stderr, "Verifying ECDSA: msg=%s pubkey=%s sig=%s\n", m.ToString().c_str(), q.ToString().c_str(), s.ToString().c_str());
+    if (!s.Verify(q, m))
+        return 0;
+    return 1;
+}
+
+}
\ No newline at end of file
diff --git a/tests.cpp b/tests.cpp
index 20273ed..2f342b6 100644
--- a/tests.cpp
+++ b/tests.cpp
@@ -1,10 +1,10 @@
 #include <assert.h>
 
-#include "num.h"
-#include "field.h"
-#include "group.h"
-#include "ecmult.h"
-#include "ecdsa.h"
+#include "num.cpp"
+#include "field.cpp"
+#include "group.cpp"
+#include "ecmult.cpp"
+#include "ecdsa.cpp"
 
 using namespace secp256k1;