Interesting examples:
- Matplotlib doesn't support Python 3.x at all yet.
- Mercurial doesn't support Python 3.x, and they don't have any plans to do so. (Which is another point in favor switching Sage to GIT.)
This is my blog about things related to Sage.
Totals grouped by language (dominant language first): python: 530370 (96.41%) ansic: 14538 (2.64%) cpp: 5188 (0.94%)
SLOC Directory SLOC-by-Language (Sorted) 88903 rings python=87720,cpp=1183 72913 combinat python=71629,cpp=1284 47747 schemes python=46255,cpp=1492 39815 graphs python=28377,ansic=11438 31540 matrix python=31540 31019 modular python=31012,ansic=7 24475 libs python=21171,ansic=2845,cpp=459 20517 misc python=20383,ansic=134 18006 interfaces python=18006 17577 geometry python=16936,cpp=641 12775 categories python=12775 12093 server python=12093 11971 groups python=11971 11961 plot python=11961 10686 crypto python=10686 9920 modules python=9920 8389 symbolic python=8260,cpp=129 8150 algebras python=8150 7260 ext python=7198,ansic=62 7093 structure python=7093 6364 coding python=6364 5670 functions python=5670 5249 homology python=5249 4798 numerical python=4798 4323 quadratic_forms python=4323 3919 gsl python=3919 3911 calculus python=3911 3879 sandpiles python=3879 3003 sets python=3003 2647 databases python=2647 2074 logic python=2074 1736 finance python=1736 1608 games python=1608 1465 monoids python=1465 1435 tests python=1383,ansic=52 1370 stats python=1370 971 interacts python=971 959 tensor python=959 906 lfunctions python=906 308 parallel python=308 275 probability python=275 219 media python=219 197 top_dir python=197
#!/usr/bin/env python
import os, shutil
for dirpath, dirnames, filenames in os.walk('.'):
for f in filenames:
if f.endswith('.pyx') or f.endswith('.pxi'):
print f
shutil.move(os.path.join(dirpath, f),
os.path.join(dirpath, os.path.splitext(f)[0] + '.py'))
sage: crt(2, 1, 3, 5) # Chinese Remainder Theorem
11
sage: crt? # ? = documentation and examples
Returns a solution to a Chinese Remainder Theorem...
...
sage: crt?? # ?? = source code
def crt(...):
...
g, alpha, beta = XGCD(m, n)
q, r = (b - a).quo_rem(g)
if r != 0:
raise ValueError("No solution ...")
return (a + q*alpha*m) % lcm(m, n)
Moreover, by browsing the Mercurial repository, you can see exactly who wrote or modified any particular line of code in the Sage library, when they did it, and why. Everything included in Sage is free and open source, and it will foreover remain that way.def python_sum2(n):
s = int(0)
for i in xrange(1, n+1):
s += i*i
return s
Then enter the following in another cell:%cython
def cython_sum2(long n):
cdef long i, s = 0
for i in range(1, n+1):
s += i*i
return s
The second implementation, despite looking nearly identical, is nearly a hundred times faster than the first one (your timings may vary).sage: timeit('python_sum2(2*10^6)')
5 loops, best of 3: 154 ms per loop
sage: timeit('cython_sum2(2*10^6)')
125 loops, best of 3: 1.76 ms per loop
sage: 154/1.76
87.5
sage: var('k, n')
sage: factor(sum(k^2, k, 1, n))
1/6*(n + 1)*(2*n + 1)*n
And now our simpler fast implementation is:def sum2(n):
return n*(2*n+1)*(n+1)/6
Just as above, we can also use the Cython compiler:%cython
def c_sum2(long n):
return n*(2*n+1)*(n+1)/6
Comparing times, we see that Cython is 10 times faster:sage: n = 2*10^6
sage: timeit('sum2(n)')
625 loops, best of 3: 1.41 microseconds per loop
sage: timeit('c_sum2(n)')
625 loops, best of 3: 0.145 microseconds per loop
sage: 1.41/.145
9.72413793103448
In this case, the enhanced speed comes at a cost, in that the answer is wrong when the input is large enough to cause an overflow:sage: c_sum2(2*10^6) # WARNING: overflow -407788678951258603Cython is very powerful, but to fully benefit from it, one must understand machine level arithmetic data types, such as long, int, float, etc. With Sage you have that option.
sage: R.We compare this with Magma (V2.17-4), which has a more ad hoc coercion system:= PolynomialRing(ZZ) sage: f = x + 1/2; f x + 1/2 sage: parent(f) Univariate Polynomial Ring in x over Rational Field
> R:= PolynomialRing(IntegerRing()); > x + 1/2 ^ Runtime error in '+': Bad argument types Argument types given: RngUPolElt[RngInt], FldRatElt
sage: f(x,y) = sin(x - y) * y * cos(x)
sage: plot3d(f, (x,-3,3), (y,-3,3), color='red')
%magma R<x> := PolynomialRing(Integers()); print Parent(x); /// Univariate Polynomial Ring in x over Integer Ring
In Sage:
R.<x> = ZZ[] parent(x) /// Univariate Polynomial Ring in x over Integer Ring
x.parent() /// Univariate Polynomial Ring in x over Integer Ring
isinstance(ZZ, Parent) /// True
isinstance(2, Parent) /// False
Automatic Coercions:
"The primary goal of coercion is to be able to transparently do arithmetic, comparisons, etc. between elements of distinct parents."
When I used to try to get people to use Magma, perhaps the number one complaint I heard about Magma was that doing arithmetic with objects having distinct parents was difficult and frustrating.
For the first year, in Sage, there was a very simple coercion system:
That seriously sucked. E.g.,
Mod(2,7) + 6
was completely different than
6 + Mod(2,7)!
The first was Mod(1,7), and the second was the integer 8. This makes understanding code difficult and unpredictable.
So I rewrote coercion to be a bit better (this was a really painful rewrite that I mostly did myself over several hard months of work):
Then we decided that there is a canonical homomorphism Z --> Z/7Z, but there is not one Z/7Z --> Z since there is no ring homomorphism in this direction, hence the above makes sense in either order.
One implication of this new model was that parent objects have to be immutable, i.e., you can't fundamentally change them after you make them. This is why in Sage you must specify the name of the generator of a polynomial ring at creation time, and can't change it. In Magma, it is typical to specify the name only later if you want.
Objects must be immutable because the canonical maps between them depend on the objects themselves, and we don't want them to just change left and right at runtime.
%magma R := PolynomialRing(RationalField(), 2); f := R.1^3 + 3*R.2^3 - 4/5; print f; /// $.1^3 + 3*$.2^3 - 4/5 [ $.1, $.2 ]
%magma AssignNames(~R, ["x", "y"]); print f; [R.1, R.2] /// x^3 + 3*y^3 - 4/5 [x, y]
%magma AssignNames(~R, ["z", "w"]); print f; /// z^3 + 3*w^3 - 4/5
R = PolynomialRing(QQ) /// TypeError: You must specify the names of the variables.
R.<x,y> = PolynomialRing(QQ) f = x^3 + 3*y^3 - 4/5; f /// x^3 + 3*y^3 - 4/5
Note: In Sage, you can can use a with block to temporarily change the names if you really need to for some reason. This is allowed since at the end of the with block the names are guaranteed to be changed back.
with localvars(R, ['z','w']):
print f
print "back?", f
///
z^3 + 3*w^3 - 4/5
back? x^3 + 3*y^3 - 4/5
But this new model had a major problem too, e.g., if x in Z[x] then "x + 1/2" would FAILS! This is because 1/2 does not coerce into Z[x] (the parent of x), and x does not coerce into Q (the parent of 1/2).
Maybe the implementors of Magma have the answers? Evidently not.
%magma R<x> := PolynomialRing(Integers()); x + 1/2; /// Runtime error in '+': Bad argument types Argument types given: RngUPolElt[RngInt], FldRatElt
Robert Bradshaw did though, and now it is in Sage:
R.<x> = ZZ[] x + 1/2 /// x + 1/2
His new design is (for the most part) what Sage actually uses now.
He launched an effort in 2008 (see the Dev Days 1 Wiki) to implement a rewrite of the coercion model to his new design. This ended up swallowing up half the development effort at the workshop, and was a massive amount of work, since every parent structure and element had to have some modifications made to it.
This meant people changing a lot of code all over Sage that they didn't necessarily understand, and crossing their fingers that the doctest test suite would catch their mistakes. This was SCARY. After much work, none of this went into Sage. It was just way too risky. This failure temporarily (!) burned out some developers.
Robert Bradshaw, on the other hand, persisted and came up with a new approach that involved migrating Sage code gradually. I.e., he made it so that the old coercion model was still fully supported simultaneously with the new one, then he migrated a couple of parent structures, and got the code into Sage. I'm sure not everything is migrated, even today. There are two points to what he did:
The coercion model is explained here: http://sagemath.org/doc/reference/coercion.html