Point addition on secp256k1. Exploring other alternatives for practical purposes

Question

#include "gmpxx.h"
#include <thread>
#include <string>
using namespace std;

// g++ -o a2.out testing.cpp -lgmp -lgmpxx -lpthread -O3
// time ./a2.out

int doer(int num)
{
    mpz_class data_x;
    mpz_class data_y;
    
        if (num == 61)
    {
        data_x = "72916064921647370085016635355356231938499058193543903073090487639385317536445";
        data_y = "57938825378424166518393844761713365209225024129681734023905584584640010510916";
    }

    else if (num == 64)
    {
        data_x = "8389021145572919717682435077323038275979373933951718417694102392981087113933";
        data_y = "8033071068995636772243200172423887768231133338759301226154161344325502516588";
    }

    mpz_class p("115792089237316195423570985008687907853269984665640564039457584007908834671663");

    mpz_class s_x("55066263022277343669578718895168534326250603453777594175500187360389116729240");

    mpz_class s_y("32670510020758816978083085130507043184471273380659243275938904335757337482424");

    mpz_class m, m1, k1, point_x, point_y;
    point_x = data_x;
    point_y = data_y;
    int j;

    for (j = 1; j < 2000000001; j += 1)
        {
                k1 = point_x - s_x;
                mpz_invert(m1.get_mpz_t(), k1.get_mpz_t(), p.get_mpz_t());
                m = (point_y - s_y) * m1;
                point_x = (m * m - s_x - point_x) % p;
                point_y = (-(s_y + m * (point_x - s_x))) % p;
                if (point_y < 0)
                {
                    point_y = p + point_y;
                }
                //Other Basic if-else with point_x and point_y ~zero exec time
        }
    return 0;
}

int main()
{
    thread th[4];
    th[0] = thread(doer,61);
    th[1] = thread(doer,64);
    th[0].join();
    th[1].join();

    return 0;
}

The above code increments a Point(point_x, point_y) on the secp256k1 curve by 1(s_x, s_y) for each loop iteration. GMP Library is used to handle very large numbers.

Execution time is around 150 minutes. I removed the unnecessary checks done while point addition, shifted operations from python to c++ (which worked a treat), implemented multithreading with O3 optimization (which increased speed by 20%). Performance analysis with perf revealed that the issue seems to be with mpz_invert() which takes around 65% of the total time.

I tried to write a function for mpz_invert() (see below) but it is 10x times slower than mpz_invert().

        k1 = point_x - s_x;
        //mpz_invert(m1.get_mpz_t(), k1.get_mpz_t(), p.get_mpz_t());
        p_temp = p;
        if (k1 < 0)
        {
            k1 = k1 + p;
        }
        nik_2 = -(p_temp / k1);

        nik_temp_2 = k1;
        //k1 = p_temp - (p_temp / k1) * k1
        k1 = p_temp + (nik_2) * k1;
        p_temp = nik_temp_2;

        nik_4 = 1-(nik_2 * (p_temp / k1));
        nik_temp_2 = k1;
        k1 = p_temp - (p_temp / k1) * k1;
        p_temp = nik_temp_2;

        while(k1 != 1)
        {
            q = p_temp / k1;
            nik_temp_2 = nik_2 - nik_4*(q);
            nik_2 = nik_4;
            nik_4 = nik_temp_2;
            nik_temp_2 = k1;
            k1 = p_temp - (q) * k1;
            p_temp = nik_temp_2;
        }

        if (nik_4 < 0)
        {
            nik_4 = p + nik_4;
        }

Is there anything I can do here or have I reached a dead end where I have to put my resources on the hardware side? Putting up multiple cloud instances does not work as the operation is needed to be scaled-up and the required cloud instances are just not practical with the current speed of execution. Any other alternatives to GMP which works faster? Anything with Parallel processing like CUDA? Can OpenCL have a huge impact on execution time here?