This site is supported by donations to The OEIS Foundation.
User:Peter Luschny/SageForOEIS
Tools, tips and tricks for using
SAGE in the context of OEIS
Get the software
Browse the docs
search docs, reference, oeis module
Trimming a sequence.
Historically the OEIS entries had 3 term lines with a total number of characters (including signs). Today the length of 260 characters is recommended for the 'DATA' field in a submission to OEIS. The function below trims a given list to this length.
def OEIStrim(L) :
n = 0
leng = 0
T = []
for l in L :
s = 1 if l < 0 else 0
leng = leng + s + len(str(l))
if leng > 260 : break
leng += 2
T.append(l)
n += 1
print n, "terms"
return T
Sample usage: OEIStrim([(-n)^n for n in (0..100)])
Using the OEIS interface
In Sage you can use the builtin 'sloane.Axxxxxx' functions. For example a fast implementation of A102467 is
def is_A102467(n) :
return bool(sloane.A001221(n) <> 1 and sloane.A001222(n) <> 2)
def A102467_list(n) :
return [k for k in (1..n) if is_A102467(k)]
The functions of the 'sloane'-sequences that Sage computes are documented here. If you are curious how these functions are implemented look at github here (Jaap Spies wrote most sequences).
Memoization
In Sage memoization is activated by applying the memoization decorator 'CachedFunction' to the function. Our first example is lame: the factorial function.
@CachedFunction
def f(x):
return x*f(x-1) if x > 1 else 1
for i in (0..9) : print f(i)
Now if you compute, say, f(1000), you not only have executed 1000 multiplications but also build up a cache holding 1000 more or less big integers. If you do not subsequently use these values in further calculations this is a huge waste of resources. Therefore this is a potential misleading example (taken from the wiki). Techniques of dynamic programming are no substitutes for good algorithms. If you want to know how to compute the factorial function efficiently you might look here.
If only a loop over the values is required the economic way is to use a generator (this does not require an internal remember table).
def F(max) :
n = 0
f = 1
while n <= max:
yield f
n += 1
f *= n
for f in F(9) : print f
A much more meaningful example for the use of memoization then the factorial function is provided by the exponential transform:
def exptrans(s) :
@CachedFunction
def g(n) :
return add(binomial(n-1,k-1)*s(k)*g(n-k) for k in (1..n)) if n > 0 else 1
return g
The exponential transform maps a sequence s to a sequence t. For example let s be the constant sequence 1.
def s(n) :
return 1
t = exptrans(s)
[t(n) for n in (0..9)]
out> [1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147]
Or let s be the constant sequence -1.
def s(n) :
return -1
t = exptrans(s)
[t(n) for n in (0..9)]
out> [1, -1, 0, 1, 1, -2, -9, -9, 50, 267]
A more succinct form to express this is the use of the lambda notation, which creates a callable Python function that sends m to 1.
A000110 = exptrans(lambda m : 1) [A000110(n) for n in (0..9)]
out> [1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147]
Or even shorter:
[exptrans(lambda m : 1)(n) for n in (0..9)]
More examples of the exponential transform can be found here .
Flattening an iterable
An often needed utility function is 'flatten'.
def flatten(I) :
L = []
for i in I :
if hasattr(i, "__iter__") :
L.extend(flatten(i))
else :
L.append(i)
return L
A sample use is to flatten a ruler given by its type (see Wichmann rulers).
Scientific format
A printing routine for arbitrary precision (MPFR) reals.
def SciFormat(f, n) :
s = f.str()
t = ""
count = 0
w = true
for c in s :
if c <> 'e' and w :
if count < n+2 :
t += c
count += 1
else :
t += c
w = false
return t
# For example:
b = 2.2379923576571269975e4767529
for i in (0..20): print SciFormat(b, i)
A sample use can be seen in Asymptotics of Bernoulli Numbers where only valid digits of an asymptotic approximation are displayed.
Benchmarking
Let us compute the Bell numbers with their generating function ...
def A000110_list(n) :
R.<t> = PowerSeriesRing(QQ, default_prec=n+1)
G = exp(exp(t)-1)
return G.ogf().padded_list()
... and with the implementation given in the database, A000110 :
def A000110_list(m):
A = [0 for i in range(0, m)]
A[0] = 1
R = [1, 1]
for n in range(1, m):
A[n] = A[0]
for k in range(n, 0, -1):
A[k-1] += A[k]
R.append(A[0])
return R
The timing:
time a = A000110_list(1000)
Time: CPU 4.31 s, Wall: 5.56 s (in the case of generating function)
Time: CPU 0.99 s, Wall: 1.35 s (in the case of the OEIS implementation)
Thus the OEIS implementation is 4 times faster than the implementation
based on the generating function ;-)
Fredrik Johansson implemented Bell numbers in mpmath. He writes on his blog: "I did a bit of benchmarking for computing large Bell numbers Bn, and this approach appears to be much faster than what other systems use. " So let us look at his implementation.
from mpmath import *
def bellexact(n,x=1) :
mp.prec = 20
size = int(log(bell(n,x),2))+10
mp.prec = size
return int(bell(n,x)+0.5)
def A000110_list(n) :
return [bellexact(k) for k in (0..n)]
time a = A000110_list(1000)
Time: CPU 19.69 s, Wall: 27.72 s
Of course Johansson's algorithm was optimized for single evaluations and not for computing a list of values. However, the example shows how easily Phython libraries can be used within Sage as in this case mpmath.
Still we have another way to compute the Bell numbers with Sage: expnums! (The paper of E. T. Bell had the title "Exponential numbers".)
time a = expnums(1001,1)
Time: CPU 0.10 s, Wall: 0.19 s
This Sage class (author Nick Alexander) is the clear winner. On the other hand our implementation is a plain vanilla Python script whereas the 'expnums' function is compiled code based on the ultra-fast MPIR (formerly GMP) library.
There is yet another function witch computes Bell numbers, in the class sage.combinat.combinat, called bell_number which wraps GAP's Bell.
def A000110_list(n) :
return [bell_number(k) for k in (0..n)]
time a = A000110_list(1000)
This time we have to wait for more than three minutes! The stopwatch says
Time: CPU 5.25 s, Wall: 182.39 s
Well, compared to that our little 10-liner which took about one second shines ;-)
Anyway, these examples show that Sage is a big system that provides a multitude of ways to perform a task. And secondly that every major sequence in OEIS should have two implementations: one to compute an arbitrary value of the sequence and another to compute the list of values [a0, a1, a2, ..., an] efficiently.
Using the database locally
The function SloaneEncyclopedia.install() installs a compact copy of the database locally. The files are stripped.gz and names.gz. It is sufficient to call this function only once in a while to update it. To do so use it in the form SloaneEncyclopedia.install( overwrite=True ). More information here.
SloaneEncyclopedia.install()
For various reasons I found it advantageous to wrap the library function SloaneEncyclopedia.find() like this:
def OEIS_ID(T, num) :
# T is the seq given as a list to identify.
# num is the number of possible matches to be returned.
if len(T) < 3:
return []
S = T[1:min(len(T),20)]
R = SloaneEncyclopedia.find(S, num)
if R == [] :
U = filter(lambda z: z != 0 , S)
del U[0]
R = SloaneEncyclopedia.find(U, num)
A = [r[0] for r in R]
s = []
for a in A :
sa = str(a)
L = "A000000"[0:7-len(sa)]
s.append(L + sa)
return s
Example call:
OEIS_ID([1,2,3,4,5,6,7,8,9,10], 4)
out> ['A000027', 'A000401', 'A000926', 'A001477']
For instance I made good use of this function while computing a list of positive definite quadratic forms.
Using Pari
Using published Pari scripts in Sage is easy. For example consider computing the number of subgroups of the Abelian group (C_2)^n.
in> gp.read(get_remote_file('http://www.cse.sc.edu/~maxal/gpscripts/nsg.gp'))
out> Attempting to load remote file:
http://www.cse.sc.edu/~maxal/gpscripts/nsg.gp
Loading: [.]
in> [gp('numsubgrp(2,vector(%d,i,1))' %(n)) for n in (0..9)]
out> [1, 2, 5, 16, 67, 374, 2825, 29212, 417199, 8283458]
This is A006116. You do not need to install Pari to do this. Pari installs automatically with Sage and works together with Sage without further ado.
We could also write a little wrapper function for this:
def A006116(n) : """ Implementation of the formula for the number of subgroups of an abelian group given the paper G.A Miller "On the Subgroups of an Abelian Group", The Annals of Mathematics, 2nd Ser., Vol. 6, No. 1 (Oct., 1904), pp. 1-6 (c) 2007 by Max Alekseyev """ gp.read(get_remote_file('http://home.gwu.edu/~maxal/gpscripts/nsg.gp')) return gp('numsubgrp(2,vector(%d,i,1))' %(n))
A disadvantage of this method is that you have to be online to use this function. However you can push things a little further and incorporate the code directly into a wrapper — provided you do not infringe a copyright. The next example shows how to do.
in> gp.read(get_remote_file('http://luschny.de/hacking/gp/A094348.gp'))
out> Attempting to load remote file: http://luschny.de/hacking/gp/A094348.gp Loading: [.]
in> gp('A094348')
out> (lim)->local(n,i,R,A,len,count,change,high);R=vector(500);A=vector(50\ );A[1]=1;A[2]=2;A[3]=4;A[4]=6;A[5]=12;count=5;high=0;n=12;while(n<=li\ m,d=divisors(n);len=length(d);change=0;for(i=1,min(len,high),if(R[i]>\ d[i],R[i]=d[i];change=1));if(len>high,for(i=high+1,len,R[i]=d[i]);hig\ h=len);if(change,count++;A[count]=n);n+=12;);vector(count,i,A[i])
Now you can copy the code into a wrapper. Just replace the "(lim)->" at the head of the code by "A094348(lim)=".
def A094348(lim) : gp('A094348(lim)=local(n,i,R,A,len,count,change,high);\ R=vector(500);A=vector(50);A[1]=1;A[2]=2;A[3]=4;A[4]=6;A[5]=12;\ count=5;high=0;n=12;while(n<=lim,d=divisors(n);len=length(d);\ change=0;for(i=1,min(len,high),if(R[i]>d[i],R[i]=d[i];change=1));\ if(len>high,for(i=high+1,len,R[i]=d[i]);high=len);\ if(change,count++;A[count]=n);n+=12;);vector(count,i,A[i])') return gp('A094348(%d)'%(lim))
Calling A094348(7207200) returns after a while
[1, 2, 4, 6, 12, 24, 36, 48, 60, 72, 120, 180, 240, 360, 420, 720, 840, 1260, 1680, 2520, 5040, 7560, 10080, 15120, 20160, 25200, 27720, 30240, 45360, 50400, 55440, 83160, 110880, 166320, 221760, 277200, 332640, 360360, 498960, 554400, 665280, 720720, 1081080, 1441440, 2162160, 2882880, 3603600, 4324320, 6486480, 7207200]
Strictly speaking the wrapper is not necessary. After 'gp.read(...)' you can for example call gp('A094348(2000)') directly. But it is good software practice to wrap the call. Sage speaks Python and all the layers Sage rests on should be encapsulated. So when you decide one day to implement A094348 in Python or base it on a different library there will be only one place where you have to change things and not throughout your code.
On a side note: Matthew Vandermast's A094348 is a very interesting sequence. See the comments of David Wasserman. A more sophisticated Pari script can be found on seqcomp and also a discussion A good algorithm for A094348?.
Compute in the cloud, publish on the web.
Go to SMC (Sage Math Cloud) create an account (which is free) and you can instantly start computing with Sage without any installation of software. Many different formats are preinstalled on SMC. So you can not only work with Sage notebooks but for instance also with IPython-notebooks which I prefer lately. To make Sage available in an IPython-notebook you just have to include the command "%load_ext sage" once at the start of the IP-notebook. For example my blog post on perfect and optimal rulers is nothing but a Sage notebook exported to Html and the corresponding IP-notebook can be browsed here and downloaded here.
More information and help can be found here sagemath-cloud-help and here: cloud-forum.
