Manejo de números enormes

Pow en C++ da error porque trabaja con floats, es solo una aproximación.

#include <iostream>
#include <math.h>
using namespace std;
constexpr long int_pow(long b, long e)
{
    return (e == 0) ? 1 : b * int_pow(b, e - 1);
}


long numero =   int_pow(1564557354, 2); //2447839713955481316;
int modulo = 23;

struct xgcd_return
{
   long x;
   long y;
 long GCD;
};


int invmod( long a_Dvd, int value_m) {
    /*
    Return 1 / a (mod m).
    @a and @m must be co-primes --> a * x + m * y = 1
    */
    if (value_m < 2) {
        printf("modulus must be greater than 1");
        return 0;
    }

    /*
    Extented Euclid GCD algorithm.
    Return (x, y, GCD) : a * x + b * y = gcd(a, b)
    */
    xgcd_return solution = { 0, 0, 0 };

    if (a_Dvd == 0) {
        //solution = { 0, 1, value_m };
        printf("\nsolution = { 0, 1, %i };", value_m);
    }
    else if (value_m == 0) {
        //solution = { 1, 0, a_Dvd };
        printf("\nsolution = { 1, 0, %lu };", a_Dvd);
    }
    else {
        long s_pp = 1, t_pp = 0; //s_previoprevio, t_previoprevio
        long m_dvr = value_m, s_pv = 0, t_pv = 1; //Divisor, s_previo, t_previo
        printf("q, r, s, t");
        printf("\n   %ld, %ld, %ld", a_Dvd, s_pp, t_pp);
        printf("\n   %ld, %ld, %ld", m_dvr, s_pv, t_pv);

        while(m_dvr) {
            long qtt = a_Dvd / m_dvr; //Nos quedamos con la parte entera
            //int rmr = a_Dvd % m_dvr; //Remainder
             long s_c = ( long)s_pp -  ( long)qtt * ( long) s_pv;
          long  t_c =( long) t_pp - ( long)qtt * ( long)t_pv;
            printf("\n%ld, %ld, %ld,%ld ", qtt, a_Dvd%m_dvr, s_c, t_c);

            //en la siguiente iteracion

            long a = a_Dvd % m_dvr; //Remainder
            a_Dvd=m_dvr;
            m_dvr=a;
            s_pp = s_pv;
            t_pp = t_pv;
            s_pv = s_c;
            t_pv = t_c;

            //Solucion
            solution.x = s_pp;
            solution.y = t_pp;
            solution.GCD = a_Dvd;
        }
    }

    //printf("\n(x, y, g): %i, %i, %i", solution.x,solution.y,solution.GCD);
    if (solution.GCD != 1) {
        printf("\nno invmod for given @a and @n");
        return 0;
    }
    else {
        if (solution.x<0) { //si el resultado es negativo C++ no hace bien el modulo
            //printf("\nx: %i",solution.x);
            //printf("\nm: %i",value_m);
            //printf("\nResto de %i / %i = %i",solution.x,value_m,(solution.x+value_m));
            return ((solution.x+value_m));
        }
        else {
            //printf("\nx: %i",solution.x);
            //printf("\nm: %i",value_m);
            //printf("\nResto de %i / %i = %i",solution.x,value_m,solution.x % value_m);
            return (solution.x % value_m);
        }
    }
}

int main()
{
    int resultado = invmod (numero, modulo);
    printf("\ninvmod (%lu, %i)= ", numero, modulo);
    printf("%i",resultado);

    return 0;
}

screen04