Logging to /Users/s/tmp/pari-16.09
GPRC Done.

GP/PARI CALCULATOR Version 2.12.1 (development 25476-78c5f8d2c)
 i386 running darwin (x86-64/GMP-6.2.0 kernel) 64-bit version
compiled: Jun 21 2020, Apple clang version 11.0.0 (clang-1100.0.33.17)
                   threading engine: single
        (readline v8.0 enabled, extended help enabled)

            Copyright (C) 2000-2020 The PARI Group

PARI/GP is free software, covered by the GNU General Public 
License, and comes WITHOUT ANY WARRANTY WHATSOEVER.

Type ? for help, \q to quit.
Type ?17 for how to get moral (and possibly technical) support.

parisizemax = 900001792, primelimit = 1000000
(15:17) gp > ?moebius
moebius(x): Moebius function of x.

(15:17) gp > ??moebius
moebius(x):

   Moebius mu-function of |x|;  x must be a non-zero integer.

   The library syntax is long moebius(GEN x).

(15:17) gp > ?intnum
intnum(X=a,b,expr,{tab}): numerical integration of expr from a to b with 
respect to X. Plus/minus infinity is coded as +oo/-oo. Finally tab is 
either omitted (let the program choose the integration step), a 
non-negative integer m (divide integration step by 2^m), or data 
precomputed with intnuminit.

(15:27) gp > intnum(x=-oo,oo,exp(x^2))
  ***   at top-level: intnum(x=-oo,oo,exp(x^2))
  ***                                 ^---------
  *** exp: overflow in expo().
  ***   Break loop: type 'break' to go back to GP prompt
break> break

(15:28) gp > intnum(x=-oo,+oo,exp(x^2))
  ***   at top-level: intnum(x=-oo,+oo,exp(x^2))
  ***                                  ^---------
  *** exp: overflow in expo().
  ***   Break loop: type 'break' to go back to GP prompt
break> break

(15:28) gp > \p2
   realprecision = 19 significant digits (2 digits displayed)
(15:28) gp > intnum(x=-oo,+oo,exp(x^2))
  ***   at top-level: intnum(x=-oo,+oo,exp(x^2))
  ***                                  ^---------
  *** exp: overflow in expo().
  ***   Break loop: type 'break' to go back to GP prompt
break> break

(15:28) gp > ??intnum
intnum(X = a,b,expr,{tab}):

   Numerical  integration  of expr on ]a,b[ with respect to X,  using the
double-exponential method, and thus O(Dlog D) evaluation of the integrand
in precision D. The integrand may have values belonging to a vector space
over  the  real  numbers;    in  particular,  it can be complex-valued or
vector-valued.   But it is assumed that the function is regular on ]a,b[.
If  the  endpoints  a and b are finite and the function is regular there,
the situation is simple:

   ? intnum(x = 0,1, x^2)
   %1 = 0.3333333333333333333333333333
   ? intnum(x = 0,Pi/2, [cos(x), sin(x)])
   %2 = [1.000000000000000000000000000, 1.000000000000000000000000000]

An  endpoint  equal  to  +/-  oo  is  coded  as +oo or -oo,  as expected:

   ? intnum(x = 1,+oo, 1/x^2)
   %3 = 1.000000000000000000000000000
/*-- (type RETURN to continue) --*/

In  basic  usage,   it  is  assumed  that  the function does not decrease
exponentially fast at infinity:

   ? intnum(x=0,+oo, exp(-x))
     ***   at top-level: intnum(x=0,+oo,exp(-
     ***                 ^--------------------
     *** exp: overflow in expo().

We shall see in a moment how to avoid that last problem, after describing
the last optional argument tab.

The  tab.  argument  The  routine  uses  weights  w_i,   which are mostly
independent of the function being integrated,  evaluated at many sampling
points  x_i  and approximates the integral by sum w_i f(x_i).   If tab is

* a non-negative integer m,  we multiply the number of sampling points by
2^m, hopefully increasing accuracy.  Note that the running time increases
roughly by a factor 2^m.   One may try consecutive values of m until they
/*-- (type RETURN to continue) --*/
give the same value up to an accepted error.

*  a  set  of  integration  tables  containing precomputed x_i and w_i as
output by intnuminit.  This is useful if several integrations of the same
type  are  performed  (on the same kind of interval and functions,  for a
given  accuracy):  we  skip  a  precomputation  of  O(Dlog  D) elementary
functions  in  accuracy  D,   whose  running  time  has the same order of
magnitude  as  the  evaluation  of the integrand.   This is in particular
useful for multivariate integrals.

Specifying  the  behavior  at  endpoints.  This  is done as follows.   An
endpoint a is either given as such (a scalar,  real or complex, oo or -oo
for  +/-  oo ),  or as a two component vector [a,alpha],  to indicate the
behavior of the integrand in a neighborhood of a.

   If a is finite, the code [a,alpha] means the function has a
singularity of the form (x-a)^{alpha}, up to logarithms. (If alpha \ge 0,
we only assume the function is regular, which is the default assumption.)
If  a wrong singularity exponent is used,  the result will lose decimals:
/*-- (type RETURN to continue) --*/

   ? c = -9/10;
   ? intnum(x=0, 1, x^c)         \\ assume x^{-9/10} is regular at 0
   %1 = 9.9999839078827082322596783301939063944
   ? intnum(x=[0,c], 1, x^c)  \\ no, it's not
   %2 = 10.000000000000000000000000000000000000
   ? intnum(x=[0,c/2], 1, x^c) \\ using a wrong exponent is bad
   %3 = 9.9999999997122749095442279375719919769

   If a is +/- oo ,  which is coded as +oo or -oo,  the situation is more
complicated, and [+/-oo,alpha] means:

* alpha = 0  (or no alpha at all,  i.e.   simply +/-oo)  assumes that the
integrand tends to zero moderately quickly, at least as O(x^{-2}) but not
exponentially fast.

*  alpha  >  0 assumes that the function tends to zero exponentially fast
approximately  as  exp(-alpha x).   This includes oscillating but quickly
decreasing functions such as exp(-x)sin(x).
/*-- (type RETURN to continue) --*/

   ? intnum(x=0, +oo, exp(-2*x))
     ***   at top-level: intnum(x=0,+oo,exp(-
     ***                 ^--------------------
     *** exp: exponent (expo) overflow
   ? intnum(x=0, [+oo, 2], exp(-2*x))  \\ OK!
   %1 = 0.50000000000000000000000000000000000000
   ? intnum(x=0, [+oo, 3], exp(-2*x))  \\ imprecise exponent, still OK !
   %2 = 0.50000000000000000000000000000000000000
   ? intnum(x=0, [+oo, 10], exp(-2*x)) \\ wrong exponent  ==>  disaster
   %3 = 0.49999999999952372962457451698256707393

As  the  last  exemple  shows,   the  exponential  decrease  rate must be
indicated to avoid overflow,  but the method is robust enough for a rough
guess to be acceptable.

* alpha < -1 assumes that the function tends to 0 slowly, like x^{alpha}.
Here  the  algorithm  is  less robust and it is essential to give a sharp
alpha,   unless alpha <= -2 in which case we use the default algorithm as
/*-- (type RETURN to continue) --*/
if alpha were missing (or equal to 0).

   ? intnum(x=1, +oo, x^(-3/2))         \\ default
   %1 = 1.9999999999999999999999999999646391207
   ? intnum(x=1, [+oo,-3/2], x^(-3/2))  \\ precise decrease rate
   %2 = 2.0000000000000000000000000000000000000
   ? intnum(x=1, [+oo,-11/10], x^(-3/2)) \\ worse than default
   %3 = 2.0000000000000000000000000089298011973

   The last two codes are reserved for oscillating functions.   Let k > 0
real, and g(x) a non-oscillating function tending slowly to 0 (e.g.  like
a negative power of x), then

*  alpha  =  k  *  I  assumes that the function behaves like cos(kx)g(x).

*  alpha  =  -k*  I  assumes  that the function behaves like sin(kx)g(x).

Here it is critical to give the exact value of k. If the oscillating part
is  not a pure sine or cosine,  one must expand it into a Fourier series,
/*-- (type RETURN to continue) --*/
use the above codings, and sum the resulting contributions. Otherwise you
will get nonsense.   Note that cos(kx), and similarly sin(kx), means that
very function, and not a translated version such as cos(kx+a).

Note.  If f(x)  = cos(kx)g(x) where g(x) tends to zero exponentially fast
as  exp(-alpha  x),  it is up to the user to choose between [+/-oo,alpha]
and  [+/-oo,k*  I],  but a good rule of thumb is that if the oscillations
are weaker than the exponential decrease, choose [+/-oo,alpha], otherwise
choose  [+/-oo,k*I],   although  the latter can reasonably be used in all
cases,   while  the former cannot.   To take a specific example,  in most
inverse Mellin transforms, the integrand is a product of an exponentially
decreasing and an oscillating factor.   If we choose the oscillating type
of integral we perhaps obtain the best results,  at the expense of having
to recompute our functions for a different value of the variable z giving
the transform,  preventing us to use a function such as intfuncinit.   On
the  other  hand  using the exponential type of integral,  we obtain less
accurate results, but we skip expensive recomputations.   See intfuncinit
for more explanations.

/*-- (type RETURN to continue) --*/
Power series limits. The limits a and b can be power series of
non-negative valuation,  giving a power series expansion for the integral
-- provided it exists.

   ? intnum(t=0,X + O(X^3), exp(t))
   %4 = 1.000...*X - 0.5000...*X^2 + O(X^3)
   ? bestappr( intnum(t=0,X + O(X^17), exp(t)) )- exp(X) + 1
   %5 = O(X^17)

The valuation of the limit cannot be negative since
int_0^{1/X}(1+t^2)^{-1} dt = pi/2 - sign(X)+O(X^2).

   Polynomials  and  rational functions are also allowed and converted to
power series using current seriesprecision:

   ? bestappr( intnum(t=1,1+X, 1/t) )
   %6 = X - 1/2*X^2 + 1/3*X^3 - 1/4*X^4 + [...] + 1/15*X^15 + O(X^16)

The  function does not work if the integral is singular with the constant
/*-- (type RETURN to continue) --*/
coefficient of the series as limit:

   ? intnum(t=X^2+O(X^4),1, 1/sqrt(t))
   %8 = 2.000... - 6.236608109630992528 E28*X^2 + O(X^4)

however you can use

   ? intnum(t=[X^2+O(X^4),-1/2],1, 1/sqrt(t))
   %10 = 2.000000000000000000000000000-2.000000000000000000000000000*X^2+O(X^4)

whis is translated internally to

   ? intnum(t=[0,-1/2],1, 1/sqrt(t))-intnum(t=[0,-1/2],X^2+O(X^4), 1/sqrt(t))

For  this  form  the  argument tab can be used only as an integer,  not a
table precomputed by intnuminit.

   We  shall  now see many examples to get a feeling for what the various
parameters  achieve.    All examples below assume precision is set to 115
/*-- (type RETURN to continue) --*/
decimal digits. We first type

   ? \p 115

Apparent  singularities.  In  many  cases,  apparent singularities can be
ignored.   For instance, if f(x) = 1 /(exp(x)-1) - exp(-x)/x, then int_0^
oo f(x)dx = gamma, Euler's constant Euler. But

   ? f(x) = 1/(exp(x)-1) - exp(-x)/x
   ? intnum(x = 0, [oo,1],  f(x)) - Euler
   %1 = 0.E-115

But  close  to  0  the  function  f  is computed with an enormous loss of
accuracy,   and  we  are  in fact lucky that it get multiplied by weights
which are sufficiently close to 0 to hide this:

   ? f(1e-200)
   %2 = -3.885337784451458142 E84

/*-- (type RETURN to continue) --*/
   A  more  robust  solution  is  to define the function differently near
special points, e.g. by a Taylor expansion

   ? F = truncate( f(t + O(t^10)) ); \\ expansion around t = 0
   ? poldegree(F)
   %4 = 7
   ? g(x) = if (x > 1e-18, f(x), subst(F,t,x)); \\ note that 7.18 > 105
   ? intnum(x = 0, [oo,1],  g(x)) - Euler
   %2 = 0.E-115

It  is  up to the user to determine constants such as the 10^{-18} and 10
used above.

True  singularities.  With  true singularities the result is worse.   For
instance

   ? intnum(x = 0, 1,  x^(-1/2)) - 2
   %1 = -3.5... E-68 \\ only 68 correct decimals
   
/*-- (type RETURN to continue) --*/
   ? intnum(x = [0,-1/2], 1,  x^(-1/2)) - 2
   %2 = 0.E-114 \\ better

Oscillating functions.

   ? intnum(x = 0, oo, sin(x) / x) - Pi/2
   %1 = 16.19.. \\ nonsense
   ? intnum(x = 0, [oo,1], sin(x)/x) - Pi/2
   %2 = -0.006.. \\ bad
   ? intnum(x = 0, [oo,-I], sin(x)/x) - Pi/2
   %3 = 0.E-115 \\ perfect
   ? intnum(x = 0, [oo,-I], sin(2*x)/x) - Pi/2  \\ oops, wrong k
   %4 = 0.06...
   ? intnum(x = 0, [oo,-2*I], sin(2*x)/x) - Pi/2
   %5 = 0.E-115 \\ perfect
   
   ? intnum(x = 0, [oo,-I], sin(x)^3/x) - Pi/4
   %6 = -0.0008... \\ bad
   ? sin(x)^3 - (3*sin(x)-sin(3*x))/4
/*-- (type RETURN to continue) --*/
   %7 = O(x^17)

We  may use the above linearization and compute two oscillating integrals
with  endpoints  [oo,   -I]  and [oo,  -3*I] respectively,  or notice the
obvious change of variable, and reduce to the single integral (1/2)int_0^
oo sin(x)/xdx. We finish with some more complicated examples:

   ? intnum(x = 0, [oo,-I], (1-cos(x))/x^2) - Pi/2
   %1 = -0.0003... \\ bad
   ? intnum(x = 0, 1, (1-cos(x))/x^2) \
   + intnum(x = 1, oo, 1/x^2) - intnum(x = 1, [oo,I], cos(x)/x^2) - Pi/2
   %2 = 0.E-115 \\ perfect
   
   ? intnum(x = 0, [oo, 1], sin(x)^3*exp(-x)) - 0.3
   %3 = -7.34... E-55 \\ bad
   ? intnum(x = 0, [oo,-I], sin(x)^3*exp(-x)) - 0.3
   %4 = 8.9... E-103 \\ better. Try higher m
   ? tab = intnuminit(0,[oo,-I], 1); \\ double number of sampling points
   ? intnum(x = 0, oo, sin(x)^3*exp(-x), tab) - 0.3
/*-- (type RETURN to continue) --*/
   %6 = 0.E-115 \\ perfect

Warning.  Like  sumalt,   intnum  often  assigns  a  reasonable  value to
diverging  integrals.   Use these values at your own risk!   For example:

   ? intnum(x = 0, [oo, -I], x^2*sin(x))
   %1 = -2.0000000000...

Note  the  formula  int_0^ oo sin(x)/x^sdx = cos(pi s/2)  Gamma(1-s)  , a
priori valid only for 0 < Re(s)  < 2, but the right hand side provides an
analytic continuation which may be evaluated at s = -2...

Multivariate  integration.  Using  successive univariate integration with
respect  to  different  formal  parameters,   it is immediate to do naive
multivariate integration. But it is important to use a suitable
intnuminit  to  precompute  data  for the internal integrations at least!

   For  example,  to compute the double integral on the unit disc x^2+y^2
<= 1 of the function x^2+y^2, we can write
/*-- (type RETURN to continue) --*/

   ? tab = intnuminit(-1,1);
   ? intnum(x=-1,1, intnum(y=-sqrt(1-x^2),sqrt(1-x^2), x^2+y^2, tab),tab) - Pi/2
   %2 = -7.1... E-115 \\ OK
   

The first tab is essential, the second optional. Compare:

   ? tab = intnuminit(-1,1);
   time = 4 ms.
   ? intnum(x=-1,1, intnum(y=-sqrt(1-x^2),sqrt(1-x^2), x^2+y^2));
   time = 3,092 ms. \\ slow
   ? intnum(x=-1,1, intnum(y=-sqrt(1-x^2),sqrt(1-x^2), x^2+y^2, tab), tab);
   time = 252 ms.  \\ faster
   ? intnum(x=-1,1, intnum(y=-sqrt(1-x^2),sqrt(1-x^2), x^2+y^2, tab));
   time = 261 ms.  \\ the internal integral matters most

   The  library  syntax is intnum(void *E,  GEN  (*eval)(void*,GEN),  GEN
a,GEN  b,GEN  tab,   long  prec),  where an omitted tab is coded as NULL.
/*-- (type RETURN to continue) --*/

(15:28) gp > intnum(x=[-oo,2],[+oo,2],exp(x^2))
%1 = 4.9 E925
(15:29) gp > \p5
(15:29) gp > intnum(x=[-oo,2],[+oo,2],exp(x^2))
%2 = 4.9482 E925
(15:29) gp > intnum(x=[-oo,2],[+oo,2],exp(-x^2))
%3 = 1.7725
(15:29) gp > intnum(x=-oo,+oo,exp(-x^2))
  ***   at top-level: intnum(x=-oo,+oo,exp(-x^2))
  ***                                  ^----------
  *** exp: overflow in expo().
  ***   Break loop: type 'break' to go back to GP prompt
break> break

(15:29) gp > intnum(x=[-oo,1],[+oo,1],exp(-x^2))
%4 = 1.7725
(15:29) gp > \p30
   realprecision = 38 significant digits (30 digits displayed)
(15:30) gp > intnum(x=[-oo,1],[+oo,1],exp(-x^2))
%5 = 1.77245385090551602729816748334
(15:30) gp > Pi
%6 = 3.14159265358979323846264338328
(15:30) gp > Pi/2
%7 = 1.57079632679489661923132169164
(15:30) gp > intnum(x=[-oo,1],[+oo,1],exp(-x^2/2))
%8 = 2.50662827463100050241576528481
(15:30) gp > X=intnum(x=[-oo,1],[+oo,1],exp(-x^2))
%9 = 1.77245385090551602729816748334
(15:30) gp > lindep([X,Pi])
%10 = [-74232514175457, 41881211258354]~
(15:30) gp > lindep([X,Pi,sqrt(Pi)])
%11 = [-1, 0, 1]~
(15:30) gp > X^2
%12 = 3.14159265358979323846264338328
(15:30) gp > Pi
%13 = 3.14159265358979323846264338328
(15:31) gp > gs(p)=sum(k=1,p-1,exp(2*Pi*k^2/p))
%14 = (p)->sum(k=1,p-1,exp(2*Pi*k^2/p))
(15:32) gp > gs(prime(2))
%15 = 4356.59518677357620374776678996
(15:32) gp > prime(2)
%16 = 3
(15:32) gp > gs(p)=sum(k=1,p-1,exp(2*Pi*I*k^2/p))
%17 = (p)->sum(k=1,p-1,exp(2*Pi*I*k^2/p))
(15:32) gp > prime(2)
%18 = 3
(15:32) gp > prime(3)
%19 = 5
(15:32) gp > gs(p)=[p,sum(k=1,p-1,exp(2*Pi*I*k^2/p))]
%20 = (p)->[p,sum(k=1,p-1,exp(2*Pi*I*k^2/p))]
(15:33) gp > gs(prime(1))
%21 = [2, -1.00000000000000000000000000000 + 1.88304107766078511674590954846 E-39*I]
(15:33) gp > \p5
   realprecision = 19 significant digits (5 digits displayed)
(15:33) gp > gs(prime(1))
%22 = [2, -1.0000 - 5.0166 E-20*I]
(15:33) gp > gs(prime(2))
%23 = [3, -1.0000 + 1.7321*I]
(15:33) gp > %23[2]
%24 = -1.0000 + 1.7321*I
(15:33) gp > algdep(%)
  ***   too few arguments: algdep(%)
  ***                              ^-
(15:33) gp > algdep(%,2)
%25 = x^2 + 2*x + 4
(15:33) gp > gs(prime(3))
%26 = [5, 1.2361 + 2.1684 E-19*I]
(15:33) gp > %[2]
%27 = 1.2361 + 2.1684 E-19*I
(15:34) gp > algdep(%,1)
%28 = 31622993*x - 39088169
(15:34) gp > algdep(%27,2)
%29 = x^2 + 2*x - 4
(15:34) gp > gs(prime(4))
%30 = [7, -1.0000 + 2.6458*I]
(15:34) gp > ?gs
gs =
  (p)->[p,sum(k=1,p-1,exp(2*Pi*I*k^2/p))]

(15:35) gp > algdep(%30[2])
  ***   too few arguments: algdep(%30[2])
  ***                                   ^-
(15:35) gp > algdep(%30[2],2)
%31 = x^2 + 2*x + 8
(15:35) gp > gss(n)=p=prime(n);algdep(sum(k=1,p-1,exp(2*Pi*I*k^2/p)),2)
%32 = (n)->p=prime(n);algdep(sum(k=1,p-1,exp(2*Pi*I*k^2/p)),2)
(15:36) gp > gss(1)
%33 = x^2 + 2*x + 1
(15:36) gp > gss(2)
%34 = x^2 + 2*x + 4
(15:36) gp > gss(3)
%35 = x^2 + 2*x - 4
(15:36) gp > gss(4)
%36 = x^2 + 2*x + 8
(15:36) gp > gss(5)
%37 = x^2 + 2*x + 12
(15:36) gp > gss(6)
%38 = x^2 + 2*x - 12
(15:36) gp > gss(7)
%39 = x^2 + 2*x - 16
(15:36) gp > ?gs
gs =
  (p)->[p,sum(k=1,p-1,exp(2*Pi*I*k^2/p))]

(15:37) gp > ?gss
gss =
  (n)->p=prime(n);algdep(sum(k=1,p-1,exp(2*Pi*I*k^2/p)),2)

(15:37) gp > \p
   realprecision = 19 significant digits (5 digits displayed)
(15:38) gp > gss(n)=p=prime(n);algdep(sum(k=0,p-1,exp(2*Pi*I*k^2/p)),2)
%40 = (n)->p=prime(n);algdep(sum(k=0,p-1,exp(2*Pi*I*k^2/p)),2)
(15:40) gp > gss(2)
%41 = x^2 + 3
(15:40) gp > gss(3)
%42 = x^2 - 5
(15:40) gp > gss(4)
%43 = x^2 + 7
(15:40) gp > gss(5)
%44 = x^2 + 11
(15:40) gp > gss(6)
%45 = x^2 - 13
(15:40) gp > gss(7)
%46 = x^2 - 17
(15:40) gp > gss(n)=p=n;algdep(sum(k=0,p-1,exp(2*Pi*I*k^2/p)),2)
%47 = (n)->p=n;algdep(sum(k=0,p-1,exp(2*Pi*I*k^2/p)),2)
(15:40) gp > gss(2)
%48 = x
(15:40) gp > gss(3)
%49 = x^2 + 3
(15:40) gp > gss(4)
%50 = x^2 - 4*x + 8
(15:40) gp > gss(5)
%51 = x^2 - 5
(15:40) gp > gss(6)
%52 = x
(15:40) gp > gss(7)
%53 = x^2 + 7
(15:40) gp > gss(8)
%54 = 3998607*x^2 - 22619537*x + 63977712
(15:40) gp > gss(9)
%55 = x - 3
(15:40) gp > gss(10)
%56 = x
(15:40) gp > gss(11)
%57 = x^2 + 11
(15:40) gp > gss(12)
%58 = 7338631*x^2 - 50843527*x + 176127144
(15:40) gp > gss(13)
%59 = x^2 - 13
(15:40) gp > gss(14)
%60 = x
(15:41) gp > vector(10,n,gss(prime(n)))
%61 = [x, x^2 + 3, x^2 - 5, x^2 + 7, x^2 + 11, x^2 - 13, x^2 - 17, x^2 + 19, x^2 + 23, x^2 - 29]
(15:41) gp > polcoeff(%,0)
  ***   at top-level: polcoeff(%,0)
  ***                 ^-------------
  *** polcoeff: incorrect type in polcoef (t_VEC).
  ***   Break loop: type 'break' to go back to GP prompt
break> break

(15:41) gp > polcoeff(%61,0,x)
  ***   at top-level: polcoeff(%61,0,x)
  ***                 ^-----------------
  *** polcoeff: incorrect type in polcoef (t_VEC).
  ***   Break loop: type 'break' to go back to GP prompt
break> break

(15:42) gp > vector(10,n,polcoeff(gss(prime(n)),0))
%62 = [0, 3, -5, 7, 11, -13, -17, 19, 23, -29]
(15:42) gp > -%
%63 = [0, -3, 5, -7, -11, 13, 17, -19, -23, 29]
(15:42) gp > vector(10,n,p=prime(n);polcoeff(gss(p),0))
%64 = [0, 3, -5, 7, 11, -13, -17, 19, 23, -29]
(15:43) gp > vector(10,n,p=prime(n);polcoeff(gss(p),0)/p/(-1)^((p-1)/2))
%65 = [0, -1, -1, -1, -1, -1, -1, -1, -1, -1]
(15:43) gp > vector(10,n,p=prime(n-);polcoeff(gss(p),0)/p/(-1)^((p-1)/2))
  ***   syntax error, unexpected ')': vector(10,n,p=prime(n-
  ***   );polcoeff(gss(p),0)/p/(
  ***   ^------------------------
(15:44) gp > vector(10,n,p=prime(n);-polcoeff(gss(p),0)/p/(-1)^((p-1)/2))
%66 = [0, 1, 1, 1, 1, 1, 1, 1, 1, 1]
(15:44) gp > GS(p)=z=Mod(x,polcyclo(p));
(15:45) gp > GS(7)
%68 = Mod(x, x^6 + x^5 + x^4 + x^3 + x^2 + x + 1)
(15:45) gp > GS(p)=z=Mod(x,polcyclo(p));sum(k=0,p-1,z^k)
%69 = (p)->z=Mod(x,polcyclo(p));sum(k=0,p-1,z^k)
(15:46) gp > GS(7)
%70 = Mod(0, x^6 + x^5 + x^4 + x^3 + x^2 + x + 1)
(15:46) gp > GS(p)=z=Mod(x,polcyclo(p));sum(k=0,p-1,z^(k^2)
  ***   syntax error, unexpected end of file, expecting )-> or ',' or 
  ***   ')': ...lo(p));sum(k=0,p-1,z^(k^2)
  ***                                    ^-
(15:46) gp > GS(p)=z=Mod(x,polcyclo(p));sum(k=0,p-1,z^(k^2))
%71 = (p)->z=Mod(x,polcyclo(p));sum(k=0,p-1,z^(k^2))
(15:46) gp > GS(7)
%72 = Mod(2*x^4 + 2*x^2 + 2*x + 1, x^6 + x^5 + x^4 + x^3 + x^2 + x + 1)
(15:46) gp > "ζ é um elemento de anel quociente Z[x]/(Φ_p(x)), e soma de Gauss é também um elemento de mesmo anel (apriori)"
%73 = "ζ é um elemento de anel quociente Z[x]/(Φ_p(x)), e soma de Gauss é também um elemento de mesmo anel (apriori)"
(15:48) gp > GS(7)^2
%74 = Mod(-7, x^6 + x^5 + x^4 + x^3 + x^2 + x + 1)
(15:48) gp > vector(20,n,p=prime(n);lift(GS(p)^2))
%75 = [0, -3, 5, -7, -11, 13, 17, -19, -23, 29, -31, 37, 41, -43, -47, 53, -59, 61, -67, -71]
(15:50) gp > lift(%74)
%76 = -7
(15:50) gp > vector(20,n,p=prime(n);lift(GS(p)^2)/p)
%77 = [0, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1, -1, 1, -1, -1]
(15:50) gp > vector(20,n,p=prime(n);lift(GS(p)^2)/p/(-1)^((p-1)/2))
%78 = [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
(15:51) gp > intnum(x=[-oo,1],[oo,1],exp(-x^2))
%79 = 1.7725
(15:51) gp > \p
   realprecision = 19 significant digits (5 digits displayed)
(15:52) gp > \p20
   realprecision = 38 significant digits (20 digits displayed)
(15:52) gp > \p
   realprecision = 38 significant digits (20 digits displayed)
(15:52) gp > intnum(x=[-oo,1],[oo,1],exp(-x^2))^2/Pi
%80 = 1.0000000000000000000