libgmpxx+libpari code to determine sum of four squares for semiprime p*q with qfsolve()

foursquare.cc

// p,q -> a,b,c,d with a^2+b^2+c^2+d^2==p*q
//
// f=foursquare
//
// g++ $f.cc -lgmp -lgmpxx -O3 -o $f -lpari -DPARI -Wall -Wextra -pedantic
//
// $ ./$f 999999999999999999999999999989 999999999999999999999999999983
// R = qfsolve(M)
// 2.85703s
// 223605138439657340220247124809 147944036070869714718268156597 314020169936458430849348767399 910771451642793118995923538064
// sum of 4 squares verified for 61-digit number
// $
//
#include <time.h>
#include <gmpxx.h>
#include <assert.h>

#include <iostream>

#include <pari/pari.h>


// from 2007 posting: https://pari.math.u-bordeaux.fr/archives/pari-users-0712/msg00001.html
#define LIMBS(x)  (reinterpret_cast<mp_limb_t *>((x)+2))
#define NLIMBS(x) (lgefint(x)-2)
void
GEN2mpz(mpz_t X, GEN x) {
    int64_t l = NLIMBS(x);
    X->_mp_alloc = l;
    X->_mp_size = signe(x) > 0? l: -l;
    X->_mp_d = LIMBS(x); /* shallow! */
}
GEN
mpz2GEN(mpz_t X) {
    int64_t l = X->_mp_size, lx = labs(l)+2;
    GEN x = cgeti(lx);
    x[1] = evalsigne(l > 0? 1: -1) | evallgefint(lx);
    for (int i = 2; i < lx; i++) x[i] = X->_mp_d[i-2];
    return x;
}


int main(int argc, char *argv[]) {
    assert(argc==3);
    mpz_class s, e, n, p(argv[1]), q(argv[2]);

    n = -p * q;

    pari_init(2000000000, 0);

    GEN N = mpz2GEN(n.get_mpz_t());

    GEN X = const_vec(5, gen_1);
    gel(X, 5) = N;
    assert(typ(X) == t_VEC);

    GEN M = diagonal(X);
    assert(typ(M) == t_MAT);

    std::cout << "R = qfsolve(M)" << "\n";
    clock_t start = clock();
    GEN R = qfsolve(M);
    std::cout << static_cast<float>(clock() - start) / CLOCKS_PER_SEC << "s\n";

    if (typ(R) == t_COL) {
        for(int i=1; i<=4; ++i) {
            GEN2mpz(e.get_mpz_t(), gel(R, i));
	    std::cout << abs(e) << " ";
        }
	std::cout << "\n";
        GEN S = vecsum(vecslice(vecmul(R, R), 1, 4));	
        GEN2mpz(s.get_mpz_t(), S);
	assert(s == -n);
	std::cout << "sum of 4 squares verified for "
		  << mpz_sizeinbase(s.get_mpz_t(), 10) << "-digit number\n";
    } else if (typ(R) == t_MAT) {
	std::cout << "t_MAT\n";
    } else {
	assert(typ(R) == t_INT);
	GEN2mpz(e.get_mpz_t(), R);
	std::cout << "t_INT: " << e << "\n";
    }

    exit(0);  // avoid "munmap_chunk(): invalid pointer" for "return 0;"
}

In response to this posting:
https://pari.math.u-bordeaux.fr/archives/pari-users-2310/msg00006.html

$ f=foursquare 
$ g++ $f.cc -lgmp -lgmpxx -O3 -o $f -lpari -DPARI -Wall -Wextra -pedantic
$ ./$f 999999999999999999999999999989 999999999999999999999999999983
R = qfsolve(M)
2.85703s
223605138439657340220247124809 147944036070869714718268156597 314020169936458430849348767399 910771451642793118995923538064 
sum of 4 squares verified for 61-digit number
$