a = ffgen([3,5],'a);
b = ffgen([3,10],'b);
m = ffembed(a, b);

A = ffmap(m, a)
minpoly(A) == minpoly(a)
A.mod == b.mod

rm = ffinvmap(m);
ffmap(rm,b^9+b^7+b^6+2*b)
ffmap(rm,b)

P = x^2+b*x+1; polisirreducible(P)
[r,n] = ffextend(b, P, 'c);
r
Pc = ffmap(n, P)
subst(Pc, x, r)
ffgen(r)

rn = ffinvmap(n);
nm = ffcompomap(n,m);
ffmap(n,ffmap(m,a)) == ffmap(nm, a)
ffcompomap(rm, rn) == ffinvmap(nm)

a = ffgen([3,5],'a);
f = fffrobenius(a);
ffmap(f,a) == a^3
g = fffrobenius(a, 5);
ffmap(g,a) == a
h = fffrobenius(a, 2);
h == ffcompomap(f,f)

forfactored(n=1,6,print(n))

sum(i=2^30,2^30+10^6,eulerphi(i))
##
my(s=0);forfactored(n=2^30,2^30+10^6, \
            s+=eulerphi(n));s
##

fordivfactored(4!,d,print(d))

my(s=0);fordivfactored(2^128-1,d,\
          s+=moebius(d)^2);s
##
my(s=0);fordiv(2^128-1,d,s+=moebius(d)^2);s
##

forsquarefree(n=1,10,print(n))

forperm(3, p, print(p))

forsubset(3,s,print(s))

forsubset([4,2],s,print(s))

n=36;
M1 = [21,7,7;23,2,7;20,21,31];
matdet(Mod(M1,n))
*** matdet: impossible inverse in Fl_inv:Mod(21,36).
matdetmod(M1,n)
matdet(M1)%n
M2 = matrix(400,400,i,j,if(j<i,1,i==j,prime(i+2)%
n,0));
matdet(M2)%n;
time = 3,472 ms.
matdetmod(M2,n);
time = 44 ms.

Mod(M1,n)^(-1)
matinvmod(M1,n)
[19 22  7]
[15 23 10]
[31 19  5]
M1^(-1)%n == %
N = M2^(-1)%n;
time = 8,945 ms.
N == matinvmod(M2,n)
time = 44 ms.

M3 = [20,7,31;30,25,9;20,21,31];
matimagemod(M3,n)
[6 2 5]
[0 2 1]
[0 0 1]
matkermod(M3,n)
[18  6]
[ 0 18]
[ 0  6]

E=ellinit([-25,1]);
ellratpoints(E,10)

ellratpoints(E,[10^6,1])

hyperellratpoints(x^6+x+1,100) \\ y^2 = x^6+x+1
(19/20)^6+(19/20)+1-(13109/8000)^2
hyperellratpoints(x^5-5*x^3+4*x+1,[1000,1])

E   = ellinit([0,1,1,-2,0]);
N   = ellglobalred(E)[1]
tor = elltors(E) \\ trivial
lfunorderzero(E)
G = ellgenerators(E); \\ with elldata
G = [[-1,1],[0,0]];  \\ without elldata

tam = elltamagawa(E)
reg = matdet(ellheightmatrix(E,G));
bsd = (E.omega[1]*tam)*reg
ellbsd(E)*reg
L1  = lfun(E,1,2)/2!

F=1/(1+x^2);sumnumrat(F,0)
(Pi*I*cotan(I*Pi)+1)/2
sumaltrat(F,0)
(Pi/cos(I*Pi)+1)/2
F=1+1/x^2;prodnumrat(F,1)
exp(sumpos(a=1,log(1+1/a^2)))

prodeulerrat(1+p^-2)
zeta(2)/zeta(4)
sumeulerrat(1/p^2)

zetahurwitz(Pi,Pi)
zetahurwitz(2,1) - zeta(2)
zetahurwitz(1,1) - Euler

E=ellinit([0,-1,1,-10,-20]);
L=lfuntwist(E,Mod(2,5));
lfunan(E,10)
lfunan(Mod(2,5),10)
lfunan(L,10)

nf=nfinit(polcyclo(5,'a));
E2=ellinit(E[1..5],nf);
localbitprec(64); lfun(E2,2)
L2=lfuntwist(E,Mod(4,5));
lfun(E,2)*lfun(L2,2)*norm(lfun(L,2))

p=Vecsmall([4,5,1,2,3]);
permorder(p)
permsign(p)
p^5

G=galoisinit(x^6+108);
[T,o]=galoischartable(G);
T~

binomial(7)
V1=vector(5001,i,binomial(5000,i-1));
##
V=binomial(5000);
##
V==V1

isprime(3^500+196)
##
C=ecpp(3^500+196);
##
ecppisvalid(C)

vecprod([2,3,5,7])
zeta(1/2+x+O(x^4))
laurentseries(x->zeta(x),4)
exponent(4.^111)
exponent(Pi^2)
printp([a,b;c,d])
bestapprnf(ellj((1+sqrt(-23))/2),quadhilbert(-23))