V = binomial(10)
foreach(V,n,print1(sigma(n)," "))
1 18 78 360 576 728 576 360 78 18 1
my(s=0);parforeach(V,n,sigma(n),S,s+=S);s

gcdext(135,95)
[,v,d]=gcdext(135,95); [v,d]
my([u,,d]=gcdext(135,95)); [u,d]

setdebug()
setdebug("qflll",5);
nfinit(x^8+1);
setdebug("qflll")
default(debug,0)

permcycles(Vecsmall([2,7,1,8,4,5,9,10,3,6]))
permcycles(Vecsmall([3,1,4,5,9,2,6,8,7]))

[M,C] = halfgcd(23,59)
M*[23,59]~

P = truncate(sqrt(1+x+O(x^4))); Q = x^4;
[M,C] = halfgcd(P, Q)
M*[P,Q]~

bnf = bnfinit(a^2+47);
bnr = bnrinit(bnf,77); bnr.cyc
bnr3 = bnrinit(bnf,77,,3); bnr3.cyc
bnrclassfield(bnr3)

p = nextprime(2^64);
\\ bnr = bnrinit(bnf,p); very very slow
bnr = bnrinit(bnf,p,, (2*3*5*7)^10); \\ fast
bnr.cyc
F=bnrclassfield(bnr,3)

\\ R = rnfconductor(bnf,F[1]); \\very slow
R = rnfconductor(bnf,F[1],1); \\fast
[cnd,bnr2,subg] = R; cnd
bnr2.cyc

[cnd,cndf] = rnfconductor(bnf,F[1],2);
cnd
cndf

P = x^5+20*x+16;
polgalois(P)
G = galoissplittinginit(P);
G.pol == nfsplitting(P)
galoisidentify(G)
galoisfixedfield(G,[G.group[2],G.group[6]],1)

E = ellinit([2,3]);
[N,k,vga] = lfunparams(E)

L = lfunqf(matdiagonal([1,2,3,4]));
Ld = lfundual(L);
eps = lfunrootres(L)[3]
lfunlambda(L,Pi)/lfunlambda(Ld,2-Pi)

lfuneuler(Mod(2,5),3)
lfuneuler(Mod(2,5),5)
L=lfungenus2([x^2 + x,x^3 + x^2 + 1]);
lfuneuler(L,11)
lfuneuler(L,13)

polylogmult([2,2,2],[1,-1,1])

v = [3,5,2,2];
vd = zetamultdual(v)
zetamultdual(vd)
zetamult(v)
zetamult(vd)

bz=besseljzero(Pi,1)
besselj(Pi,bz)
bz=besselyzero(2,10)
bessely(2,bz)

harmonic(10)
sum(i=1,10,1/i)
harmonic(10,3)
sum(i=1,10,1/i^3)

harm(n,k) = 
{
  my(s=k-1);
  (-1)^s*derivnum(x=0,psi(n+1+x)-psi(1+x),s)/s!;
}
harm(100000,3)
harmonic(100000,3)*1.

S = sum(i=0,100,binomial(2*i,i)/(i+1)*T^(i+1))\
    + O(T^101);
seralgdep(S,3,3)
S^2-S+T
serdiffdep(S,3,3)
(4*T-1)*S'-2*S==-1

D=4*nextprime(2^40);
norm(quadunit(D))
quadunitnorm(D)

quadunitindex(5,13)
w^7

polsubcyclo(10^7+19,7)
P=polsubcyclofast(10^7+19,7)[1]
polredabs(P)

snfrank([4,4,2], 2)
snfrank([4,4,2], 4)
snfrank([4,4,2], 8)
snfrank([4,4,2], 0)

poltomonic(9*x^2 - 1/2)
U = poltomonic(9*x^2 - 1/2, &L)
L
U / subst(9*x^2 - 1/2, x, x/L)

eulerreal(100)
eulerfrac(100)
Qfb(1,2,3).disc
solve(x=-oo,oo,exp(x)-3*x)