7 ms·

Nemo: computer algebra package for Julia

53 points by winestock 11y ago

IndianAstronaut 11y ago
Since Julia is about speed, how does this compare to Sympy speed?
- zokier 11y ago
  They have a benchmarks page where it compares quite favorably to SageMath. I don't know but I'd kinda assume that Sage uses Sympy under the hood.
  - fdej 11y ago
    Sage doesn't really use SymPy a whole lot. Most of the core algebra in Sage is based on C libraries like MPIR, Flint, Singular, Pari... 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).
- fdej 11y ago
  Not necessarily the most representative benchmark, but: >>> 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) Edit: oops, I forgot a variable which made the benchmark totally unfair!
- srean 11y ago
  > Since Julia is about speed 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.
  - IndianAstronaut 11y ago
    I honestly see little reason for its existence besides speed. R and Python are very well established.
vortico 11y ago
This is a fantastic collection of software into what seems to be a well organized and well documented package.
Cyph0n 11y ago
Holy shit! Those benchmarks look absolutely incredible, outperforming Magma by a factor of 5 on average. Incredible work. This is the kind of work that makes Julia really look like the next step for scientific computing. http://nemocas.org/benchmarks.html http://nemocas.org/benchmarks.html
- deleted 11y ago
  [deleted]
- clusterfoo 11y ago
  It's really turning out to be a beautiful general-purpose language as well. 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 :)
  - nickpsecurity 11y ago
    That's why I keep referencing it here, esp in Go or Rust discussions. Really wish the other languages would've been more like it. I hope it keeps getting better and with more use cases.
- rp21a 11y ago
  I'd like to know the CPU to compare. For example, the first one is 5 seconds on a Core i5 4570 3.2 GHz in Maple. 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.
  - wbhart 11y ago
    It was on a single core of an Opteron K10 6174 at 2.2GHz. I'm not sure what you mean by KS on the input. Nemo doesn't use KS on that benchmark (except possibly in Flint for the lowest level). This is comparable to what the other systems are doing, except Pari which uses a recursive sparse representation. 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.
wbhart 11y ago
There's a little bit more information in Fredrik J's blog post: http://fredrikj.net/blog/2015/09/finding-nemo/ http://fredrikj.net/blog/2015/09/finding-nemo/
jpfr 11y ago
The Axiom [1] developers created a special purpose dependently typed language to capture mathematical abstractions: SPAD/Aldor [2]. Maxima uses Lisp. [3] These packages contain many man-years of work. 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. [1] http://www.axiom-developer.org/ http://www.axiom-developer.org/ [2] http://www.aldor.org/ http://www.aldor.org/ [3] http://maxima.sourceforge.net/ http://maxima.sourceforge.net/
- nickpsecurity 11y ago
  I mostly just follow articles on Julia for now. Aside from it having LISP-like macros, where's a page that says it's a LISP underneath? And I didn't see dependent types mentioned in the page below despite some things sounding similar. You got a link that straight up shows dependent types in Julia or how to emulate them? http://docs.julialang.org/en/release-0.3/manual/types/ http://docs.julialang.org/en/release-0.3/manual/types/
  - jpfr 11y ago
    > Where's a page that says it's a LISP underneath? 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/metaprogramming/ http://julia.readthedocs.org/en/latest/manual/metaprogrammin... [2] https://github.com/JuliaLang/julia/tree/master/src/flisp https://github.com/JuliaLang/julia/tree/master/src/flisp [2] https://github.com/JuliaLang/julia/issues/6113 https://github.com/JuliaLang/julia/issues/6113
mrcactu5 11y ago
If it is built on PARI how does it improve upon it? Personally I want to know how to define algebraic number classes. I attempted and failed write my own algebraic number types in Python.
- wbhart 11y ago
  Pari is not used for elements of algebraic number fields. We use Antic for that. 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.
- wbhart 11y ago
  As for your second question, do you mean how do you define algebraic number fields in Nemo, or in Pari/GP? For Nemo we document it in the Antic fields section of the manual.

53 points by winestock 11y ago

IndianAstronaut 11y ago

Since Julia is about speed, how does this compare to Sympy speed?

zokier 11y ago
They have a benchmarks page where it compares quite favorably to SageMath. I don't know but I'd kinda assume that Sage uses Sympy under the hood.
- fdej 11y ago
  Sage doesn't really use SymPy a whole lot. Most of the core algebra in Sage is based on C libraries like MPIR, Flint, Singular, Pari... 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).
fdej 11y ago
Not necessarily the most representative benchmark, but: >>> 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) Edit: oops, I forgot a variable which made the benchmark totally unfair!
srean 11y ago
> Since Julia is about speed 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.
- IndianAstronaut 11y ago
  I honestly see little reason for its existence besides speed. R and Python are very well established.

vortico 11y ago

This is a fantastic collection of software into what seems to be a well organized and well documented package.

Cyph0n 11y ago

Holy shit! Those benchmarks look absolutely incredible, outperforming Magma by a factor of 5 on average. Incredible work. This is the kind of work that makes Julia really look like the next step for scientific computing. http://nemocas.org/benchmarks.html http://nemocas.org/benchmarks.html

deleted 11y ago
[deleted]
clusterfoo 11y ago
It's really turning out to be a beautiful general-purpose language as well. 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 :)
- nickpsecurity 11y ago
  That's why I keep referencing it here, esp in Go or Rust discussions. Really wish the other languages would've been more like it. I hope it keeps getting better and with more use cases.
rp21a 11y ago
I'd like to know the CPU to compare. For example, the first one is 5 seconds on a Core i5 4570 3.2 GHz in Maple. 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.
- wbhart 11y ago
  It was on a single core of an Opteron K10 6174 at 2.2GHz. I'm not sure what you mean by KS on the input. Nemo doesn't use KS on that benchmark (except possibly in Flint for the lowest level). This is comparable to what the other systems are doing, except Pari which uses a recursive sparse representation. 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.

wbhart 11y ago

There's a little bit more information in Fredrik J's blog post: http://fredrikj.net/blog/2015/09/finding-nemo/ http://fredrikj.net/blog/2015/09/finding-nemo/

jpfr 11y ago

The Axiom [1] developers created a special purpose dependently typed language to capture mathematical abstractions: SPAD/Aldor [2]. Maxima uses Lisp. [3] These packages contain many man-years of work. 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. [1] http://www.axiom-developer.org/ http://www.axiom-developer.org/ [2] http://www.aldor.org/ http://www.aldor.org/ [3] http://maxima.sourceforge.net/ http://maxima.sourceforge.net/

nickpsecurity 11y ago
I mostly just follow articles on Julia for now. Aside from it having LISP-like macros, where's a page that says it's a LISP underneath? And I didn't see dependent types mentioned in the page below despite some things sounding similar. You got a link that straight up shows dependent types in Julia or how to emulate them? http://docs.julialang.org/en/release-0.3/manual/types/ http://docs.julialang.org/en/release-0.3/manual/types/
- jpfr 11y ago
  > Where's a page that says it's a LISP underneath? 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/metaprogramming/ http://julia.readthedocs.org/en/latest/manual/metaprogrammin... [2] https://github.com/JuliaLang/julia/tree/master/src/flisp https://github.com/JuliaLang/julia/tree/master/src/flisp [2] https://github.com/JuliaLang/julia/issues/6113 https://github.com/JuliaLang/julia/issues/6113

mrcactu5 11y ago

If it is built on PARI how does it improve upon it? Personally I want to know how to define algebraic number classes. I attempted and failed write my own algebraic number types in Python.

wbhart 11y ago
Pari is not used for elements of algebraic number fields. We use Antic for that. 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.
wbhart 11y ago
As for your second question, do you mean how do you define algebraic number fields in Nemo, or in Pari/GP? For Nemo we document it in the Antic fields section of the manual.