OpenSSL implementation of Shamir Secret Sharing

I'm attempting to implement Shamir Secret Sharing using OpenSSL. I'm having a lot of trouble getting the message to decrypt!

I have tried several implementations for decryption, both this one: http://www.cs.cornell.edu/courses/cs754/2001fa/307.pdf

And this one (which references the paper above, but uses a different method):

(Note i am aware of the typo for calculating u2)

I believe my problem might be a loss of precision when doing the Lagrange Interpolation on ik/(ik-ij) using BIGNUM values to represent the key numbers. I wrote a small function that converts an int like a key number to a BIGNUM value.

I'll spare you my key generating code because i'm pretty certain that's working correctly. The following is my implementation of the JPEG that i linked (not the PDF):

int i;
for (i = 0; i < 3; i++) {
    aij[i] = BN_new();
    BN_mod_exp(aij[i], c.c1, key[i].key, p, ctx);
}


BIGNUM *ik = BN_new();
BIGNUM *ij = BN_new();
BIGNUM *denomTmp = BN_new();
BIGNUM *numTmp = BN_new();
BIGNUM *divTmp = BN_new();
BIGNUM *accum = BN_new();

int j, k;

/* Lagrange Interpolation */
for (j = 0; j < 3; j++) {
    BN_one(accum);
    for (k = 0; k < 3; k++) {
        if(j == k) continue;
        else {
            ik = int2BN(key[k].keynum); //int2BN is my function for converting ints to BNs
            ij = int2BN(key[j].keynum);

            BN_sub(denomTmp, ik, ij);

            BN_div(divTmp, NULL, ik, denomTmp, ctx);
            BN_mul(accum, accum, divTmp, ctx);
        }
    }
    cij[j] = BN_new();
    BN_mod(cij[j], accum, q, ctx); // accum % q = cij[j]
}

// Now for the second half...
int a;
u1 = BN_new();
BIGNUM *u1tmp = BN_new();
BN_one(u1);
for (a = 0; a < 3; a++) {
    BN_mod_exp(u1tmp, aij[a], cij[a], p, ctx);
    BN_mod_mul(u1, u1, u1tmp, p, ctx);
}

When i spit out the calculated values in cij[], i'm getting: 2, -5, 0 when using keys [2, 4, 5]. According to some math i did by hand it should actually be 10/3, 5, and 8/3. Is there any way to get around this?

Are there other problems with my code that you see? Thanks in advance.

You need to be doing modular arithmetic when you do sub, div, mul. i.e. BN_mod_sub, BN_mod_inverse on the denominator, and BN_mod_mul.