> 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

>