> with(numtheory); [GIgcd, bigomega, cfrac, cfracpol, cyclotomic, divisors, factorEQ, factorset, fermat, imagunit, index, integral_basis, invcfrac, invphi, iscyclotomic, issqrfree, ithrational, jacobi, kronecker, lambda, legendre, mcombine, mersenne, migcdex, minkowski, mipolys, mlog, mobius, mroot, msqrt, nearestp, nthconver, nthdenom, nthnumer, nthpow, order, pdexpand, phi, pi, pprimroot, primroot, quadres, rootsunity, safeprime, sigma, sq2factor, sum2sqr, tau, thue] > p:=nextprime(10^16); p := 10000000000000061 > for a from 2 to 20 do print(a,jacobi(a,p), a &^ ((p-1)/4) mod p); od; 2, -1, 3041443415194581 3, -1, 6958556584805480 4, 1, 10000000000000060 5, 1, 10000000000000060 6, 1, 1 7, 1, 10000000000000060 8, -1, 6958556584805480 9, 1, 10000000000000060 10, -1, 6958556584805480 11, -1, 6958556584805480 12, -1, 3041443415194581 13, 1, 10000000000000060 14, -1, 6958556584805480 15, -1, 3041443415194581 16, 1, 1 17, -1, 6958556584805480 18, -1, 6958556584805480 19, -1, 3041443415194581 20, 1, 1 > a:=6; a &^ ((p-1)/4) mod p; a := 6 1 > oddd:=(p-1)/4; oddd := 2500000000000015 > b := a &^ ((oddd+1)/2) mod p; b := 2434710841166537 > b ^ 2 mod p; 6 > aa := 5; bb := aa &^ ((p+3)/8) * 2 &^ ((p-1)/4) mod p; aa := 5 bb := 6073147792550418 > bb^2 mod p; 5 > for a from 1 to 14 do print(a, a^2 mod 15); od; 1, 1 2, 4 3, 9 4, 1 5, 10 6, 6 7, 4 8, 4 9, 6 10, 10 11, 1 12, 9 13, 4 14, 1 >