CUDA椭圆曲线点加法部分点对运算失效求助

我在CUDA内核中基于CGBN库实现椭圆曲线点运算，使用的椭圆曲线参数由SageMath定义如下：

椭圆曲线参数（SageMath定义） #

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)

目前实现的加法运算对多数点对有效，但部分特定点对（如A=Q*1001与B=Q*20相加、A=P*15与B=Q相加）返回错误结果。多次检查代码未发现问题，现附上实现函数与测试内核代码：

点加法实现函数 #

__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};
}

测试内核代码 #

__global__ void ker_add_points(cgbn_error_report_t *report, EC_point *points, EC_parameters *parameters, int32_t instances)
{
    int32_t instance;
    int32_t points_index;

    instance = (blockIdx.x * blockDim.x + threadIdx.x) / TPI;
    points_index = instance * 2;

    if (instance >= instances)
    {
        return;
    }

    context_t bn_context(cgbn_report_monitor, report, instance); // construct a context
    env192_t bn192_env(bn_context.env<env192_t>());

    dev_EC_point P0, P1;
    dev_Parameters params;

    cgbn_load(bn192_env, P0.x, &(points[points_index].x));
    cgbn_load(bn192_env, P0.y, &(points[points_index].y));

    cgbn_load(bn192_env, P1.x, &(points[points_index + 1].x));
    cgbn_load(bn192_env, P1.y, &(points[points_index + 1].y));

    cgbn_load(bn192_env, params.Pmod, &(parameters->Pmod));
    cgbn_load(bn192_env, params.a, &(parameters->a));

    dev_EC_point result = add_points(bn192_env, P0, P1, params);

    cgbn_store(bn192_env, &(points[points_index].x), result.x);
    cgbn_store(bn192_env, &(points[points_index].y), result.y);
}

所有测试通过SageMath和pytest通过ctypes调用测试函数执行，请求协助排查部分点对加法结果错误的原因。

内容的提问来源于stack exchange，提问作者drzejus