Nemo: computer algebra package for Julia
nemocas.org
nemocas.org
I've been playing it for a few days now, and it's such a joy to work with. So many great forward-thinking libraries: (Reactive, Compose), excellent, readable code quality even down to the language implementation from very smart people...
If the language stack of the future looks something like Rust for systems, Julia for applications, and Elm for UI, I'll be a pretty happy developer :)
The second problem is also not clear. Is there a Kronecker substitution on the input? If so, what are the degrees? I'd like to see the code used.
The code is in the Benchmarks test in test/Benchmark-test.jl.
It's not surprising the first benchmark is faster in Maple, since Maple can make use of a true sparse representation, possibly quite a bit of vectorisation on the processor you have and possibly multiple cores. The benchmark here uses only a dense representation, a single core and there is no explicit vectorisation (possibly none at all).
Pihanha for example will do both of those first examples in a fraction of the time that Nemo will. But again, it uses sparse representation and again can use multiple cores.
We'll do a sparse representation in Nemo later on, perhaps even wrap Pirhana.
The main purpose of the benchmark is actually to show off what the really fast Julia generics do for us, not to actually do this particular benchmark as fast as is humanly possible. In order to do that as fairly as possible, we deliberately use univariate polynomials over other univariate polynomials in all systems (except Pari, as noted), rather than dedicated multivariate polynomial rings.
That first benchmark is pretty dense, so it's not at all unreasonable to use KS.
On the other hand, if we force it to use KS for Z[x][y][z][t] down to Z[x][y], then classical down to Z[x], then Flint for Z[x] itself, it gets marginally faster, with 40.6s.
Clearly neither of these approaches is an optimal strategy, as the Maple timings show.
We have a lot of generic matrix stuff missing at present though, so that's the priority for now.
There's also things like Puiseux series and capped absolute power series and a few things like that which we can't do without.
It would be a hell of a lot of fun to optimise the hell out of both those polynomial benchmarks, but we have to be realistic and go for breadth first, rather than depth, I think.
Julia is JIT-compiled, a Lisp (under the hood) and (somewhat) dependently typed.
Is that enough to port (or even transpile) modules from the big open source CAS without an entire rewrite? Is there enough similarity between CAS to make "foreign" modules even a remote possibility?
Otherwise, Nemo will likely not achieve a great unification of math packages, since the required effort goes way beyond the resources of a small dedicated group.
Every expression is transformed internally to an AST representation that can be seen as a lisp-style s-expression. And this gets transformed, JIT-compiled, and so on. [1]
Furthermore, the Julia internals contain an entire lisp implementation (femtolisp [2]). That has influenced the metaprogramming capabilities.
> And I didn't see dependent types mentioned in the page below despite some things sounding similar.
Integer matrices can be represented as Array{Int64,2}. So the type of the matrix depends on the Int64 datatype and the dimensionality 2. But that's not necessarily what type theorists understand to be a dependent type. There was an epic discussion that led to the phrasing "dependent" being removed from the Julia docs. [2]
[1] http://julia.readthedocs.org/en/latest/manual/metaprogrammin...
[2] https://github.com/JuliaLang/julia/tree/master/src/flisp
When we first wrote Nemo in Julia, we used a form of emulated dependent typing. But the way we did this led to performance problems that we just hadn't thought about.
For example, if we did lots of operations over Z/pZ for many primes p, e.g. in multimodular algorithms, it meant that the Jit compiler had to recompile everything for every prime p, which was awfully slow.
Julia is absolutely not meant to be used that way, and it took us a long time to see that. After redesigning Nemo to use Julia much more like the way it was designed to be used, we actually ended up speeding everything up.
Reason being that the only emulation of dependent types in Julia being way to slow would be an extra argument against any claims for Julia having dependent types vs real, dependently-typed languages that don't punish you for using them. I'd rather be clear before going that far, though.
Of course you could emulate them some other way in Julia, but they are deliberately not baked into the language.
In the end we used the model that some computer algebra systems, like Sage and Magma use, of having parent objects for the rings in which the element objects live. I'll leave it up to the experts to figure out if this is dependent typing or not. It obviously works very well either way.
Thanks for the answer and links!
Edit: Just read the debate. I agree with 'tel' that the phrase dependent types has such a specific meaning, capabilities, and whole literature devoted to it that nobody should use it unless it's consistent. The Google search point was especially smart. Per that discussion, not my own analysis of Julia, Julia doesn't have dependent types & it's good they took that phrase off the site. Its different enough to justify a custom term to promote its usefulness while not misleading.
>>> from sympy import *
>>> from time import clock
>>> var('x y z t')
(x, y, z, t)
>>> f = poly(1 + x + y + z + t)
>>> p = f ** 20
>>> t1=clock(); q=p*(p+1); clock()-t1
53.828634
julia> using Nemo
julia> Rx,x=PolynomialRing(ZZ,"x");
julia> Ry,y=PolynomialRing(Rx,"y");
julia> Rz,z=PolynomialRing(Ry,"z");
julia> Rt,t=PolynomialRing(Rz,"t");
julia> p = (1+x+y+z+t)^20;
julia> @time p*(p+1);
1.593379 seconds (129.24 k allocations: 6.758 MB, 0.82% gc time)
Is it ? At least that's not the part that excites me most, although its a nice by product to have. Its a nice flexible language with good syntax and a good set of features. I still dislike 1 based indexing though.
Nemo uses some of the same libraries under the hood, but there's much less overhead. This is in part due to Julia being faster than Python (although Sage works around that to some extent by using Cython), and in part due to Nemo using tighter algorithms/code.
Note that Nemo is currently more geared towards purely algebraic computation, though. It doesn't currently have S-expressions with support for symbolic integration and the like (Sage mostly uses Pynac and Maxima for this).
Personally I want to know how to define algebraic number classes. I attempted and failed write my own algebraic number types in Python.
I should point out we actually spent some time speeding up that particular benchmark. The others are much more raw.
Having said that, we have an even faster way of doing it which will be in Nemo 0.4.
It's worth noting though that Pari doesn't have Jit compilation, so it would still be possible to beat something written in Pari/GP for example.
One of the big things with complex functionality like that in Pari is that years of mathematical knowledge have gone into it. One can't come along with some tricks like Jit compilation and expect to beat it with some simple scripts. It actually took months of work on the Antic library to get faster algebraic number field element arithmetic than what was already in Pari, and that's trivial functionality compared to the rest of Pari.
The winning strategy isn't actually to reimplement the whole of Pari (or Singular, or Gap, etc.) but to use them for complex functionality that they are actually good at, and then either improve Pari itself, or implement something even better over time, say by extending Antic.