Created
October 4, 2023 00:13
-
-
Save Hermann-SW/fc96e5b1148c72ce07cfa4e55b6c0cd3 to your computer and use it in GitHub Desktop.
libgmpxx+libpari code to determine sum of four squares for semiprime p*q with qfsolve()
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// 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;" | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
In response to this posting:
https://pari.math.u-bordeaux.fr/archives/pari-users-2310/msg00006.html