Elliptic curve addtion not working in some specific cases

Question

I posted similar question on the StackOverflow (https://stackoverflow.com/questions/78828295/cuda-elliptic-curve-addtion-not-working-in-some-specific-cases) but it might be better suited here.

I am implementing elliptic curve addition in CUDA. My current implementation works for most of point pairs but there are some specific pairs that get wrong results. I assume there might be some issue with math in my algorithm but I can't see it.

Here is function that is doing the addition operation:

__device__ dev_EC_point add_points(env192_t bn_env, const dev_EC_point &P1, const dev_EC_point &P2, const dev_Parameters &params)
{
    if (cgbn_equals(bn_env, P1.x, P2.x) && cgbn_equals(bn_env, P1.y, P2.y))
    {
        return double_point(bn_env, P1, params);
    }
env192_t::cgbn_t t2;
if (cgbn_sub(bn_env, t2, P1.x, P2.x)) // x1 - x2 mod Pmod
{
    printf("BREAKPOINT 1\n");
    cgbn_sub(bn_env, t2, P2.x, P1.x);
    cgbn_sub(bn_env, t2, params.Pmod, t2);
}


cgbn_modular_inverse(bn_env, t2, t2, params.Pmod); // 1/(x1-x2) mod Pmod

// Montgomery space

env192_t::cgbn_t x1, y1, x2, y2;

uint32_t np0;
np0 = cgbn_bn2mont(bn_env, x1, P1.x, params.Pmod);
cgbn_bn2mont(bn_env, y1, P1.y, params.Pmod);
cgbn_bn2mont(bn_env, x2, P2.x, params.Pmod);
cgbn_bn2mont(bn_env, y2, P2.y, params.Pmod);
cgbn_bn2mont(bn_env, t2, t2, params.Pmod);

env192_t::cgbn_t t1;

if (cgbn_sub(bn_env, t1, y1, y2)) // y0 - y1 mod Pmod
{
    printf("BREAKPOINT 2\n");
    cgbn_sub(bn_env, t1, y2, y1);
    cgbn_sub(bn_env, t1, params.Pmod, t1);
}

env192_t::cgbn_t s, s_sq, x3, y3, t3;

cgbn_mont_mul(bn_env, s, t1, t2, params.Pmod, np0); // s = (y1-y2)/(x1-x2) mod Pmod // tested

cgbn_mont_sqr(bn_env, s_sq, s, params.Pmod, np0); // s^2 mod Pmod // tested

cgbn_add(bn_env, t3, x1, x2); // x1 + x2

if (cgbn_sub(bn_env, x3, s_sq, t3)) // x3 = s^2 - x1 - x2 // mod Pmod
{
    printf("BREAKPOINT 4\n");
    cgbn_sub(bn_env, x3, t3, s_sq);
    cgbn_sub(bn_env, x3, params.Pmod, x3);
}

if (cgbn_sub(bn_env, t3, x1, x3)) // t3 = x1 - x3 // mod Pmod
{
    printf("BREAKPOINT 5\n");
    cgbn_sub(bn_env, t3, x3, x1);
    cgbn_sub(bn_env, t3, params.Pmod, t3);
}

cgbn_mont_mul(bn_env, t3, t3, s, params.Pmod, np0);

if (cgbn_sub(bn_env, y3, t3, y1))
{
    printf("BREAKPOINT 6\n");
    cgbn_sub(bn_env, y3, y1, t3);
    cgbn_sub(bn_env, y3, params.Pmod, y3);
}

cgbn_mont2bn(bn_env, x3, x3, params.Pmod, np0);
cgbn_mont2bn(bn_env, y3, y3, params.Pmod, np0);

// cgbn_sub(bn_env, x3, s_sq, t1);

return dev_EC_point{x3, y3};

}

CGBN library API for reference: https://github.com/NVlabs/CGBN/blob/master/docs/CGBN.md

The use-case is for very specific curve parameters:

p = 0x62CE5177412ACA899CF5
r = 0x1CE4AF36EED8DE22B99D
a = 0x39C95E6DDDB1BC45733C
b = 0x1F16D880E89D5A1C0ED1
n = 0x62CE5177407B7258DC31
P_x = 0x315D4B201C208475057D
P_y = 0x035F3DF5AB370252450A
Q_x = 0x0679834CEFB7215DC365
Q_y = 0x4084BC50388C4E6FDFAB
F = GF(p)
E = EllipticCurve(F, [a, b])
P = E(P_x, P_y)
Q = E(Q_x, Q_y)

The pairs for which I get wrong results with my implementation:

A = Q * 1001
B = Q * 20
R = A + B # got wrong answer during tests
A = P * 15
B = Q
R = A + B # got wrong answer during tests 

Some working pairs:

A = P * 10
B = Q
R = A + B
A = P * 10
B = Q * 10
R = A + B
A = P * 20
B = Q * 10
R = A + B
```

Answer 1

One guess: at the line:

 cgbn_add(bn_env, t3, x1, x2); // x1 + x2

the resulting t3 may be greater than $p$, which may lead to other intermediate values being out of the range $[0,p)$. This might cause the latter cgbn_mont_mul to misbehave.

I haven't worked with this bignum library - however I find it plausible that their Montgomery Multiplication implementation makes that assumption (and it should be easy to test...)