nf = nfinit(1+x^2);
nfeltispower(nf, -4, 4, &y)
y

\z
    doctest = 1
   ? nf = nfinit(1+x^2);
   ? nfeltispower(nf, -4, 4, &y)
   %2 = 1
   ? y
   %3 = [-1, -1]~
\z

$ cat <<EOF > qf.gp
a=2;
\qf
b=3;
EOF
$ gp qf.gp
a
b
\qf

N = 2^64; bitset(~N,0); N - 2^64
bitclear(~N, 0); N - 2^64
bitflip(~N,3); N - 2^64

sieve_slow(N) =
 {
   my(V = 1<<(N+1) - 1);
   V -= 1; V -= 2;
   for (p = 2, sqrtint(N),
      if (bittest(V,p),
        forstep (j = 2*p, N, p,
          if(bittest(V,j), V -= 2^j))));
   vecextract([0..N], V);
 }
sieve_slow(2^20) == primes([1,2^20])
##

sieve_fast(N) =
 { 
   my(V = 1<<(N+2) - 1);
   bitclear(~V, 0); bitclear(~V, 1);
   for (p = 2, sqrtint(N),
      if (bittest(V,p),
        forstep (j = 2*p, N, p, bitclear(~V,j))));
   vecextract([0..N], V-(1<<(N+1)));
 }
sieve_fast(2^20) == primes([1,2^20])
##

bnf=bnfinit(x^2+3321607); bnf.cyc
pr2=idealprimedec(bnf,2)[1];
c=bnfisprincipal(bnf,pr2,0)
matmodhnf(3*c,bnf.cyc)
bnfisprincipal(bnf,idealpow(bnf,pr2,3),0)

polclass(-1243,111)
polclass(-1243,117)
polclass(-1243,119)
polclass(-1243,129)

E = ellinit("59097a4");
c = ellmaninconstant(E,1)
d = ellmoddegree(E)
d*c^2

y = 'y;
[e,P] = ellfromeqn(6*y^2 - (x^4-1), [7,20])
E = ellinit(e);ellisoncurve(E,P)
[e,P] = ellfromeqn(6*y^2 - (x^4-1), [x,y])

E1 = ellinit([2,5]);
[ellcharpoly(E1,3),ellcharpoly(E1,101)]
Ep = ellinit([2,5],3);
ellcharpoly(Ep)

nf = nfinit(a^2+2); E = ellinit([1,1+a],nf);
pr = idealprimedec(nf,19);
C1 = ellcharpoly(E,pr[1])
C2 = ellcharpoly(E,pr[2])
C1*C2
lfuneuler(E,19)

genus2charpoly([x^6+x+7,0],23)
genus2charpoly([x^6+x+7,0],29)

nf = nfinit(y^2-2);
al = alginit(nf, [1+y, -7]);
al2 = alginit(nf, [-1, 2*y-1]);
al3 = alginit(nf, [1-y, -7]);
algisisom(al,al2)
algisisom(al,al3)

pr = idealprimedec(nf,2)[1];
algisisom(al,al3,pr)
algisisom(al,al3,1)

al = alginit(nfinit(y^2-2),[-2,-3-y]);
u = [0..7]~; \\ generates a field in al
fu = alginvol(al,u); \\ quaternion conjugation
g = algskolemnoether(al,u,fu);
algisinv(al,g)
algmul(al,g,u) == algmul(al,fu,g)\\ g*u == fu*g

theta(Pi/4,I,0)
theta(Pi/4,I,[0,0])
theta(Pi/4,I,[0,1])
theta(Pi/4,I,[1,0])
theta(Pi/4,I,[1,1])

thetanull(I,0)
theta(0,I,0)
derivnum(z=0,theta(z,I,[1,1]))
thetanull(I,[1,1])

real(elljacobi(1,5))
elljacobi(x+O(x^7),0)

f=fileextern("paste -s","w");
for(i=1,20,filewrite(f,i));
fileclose(f);