From SageMath to Julia

question

#1

With SageMath I can create a symbolic variable x and then work very
simply with it symbolically. How this is done with Julia is for me a book
with seven seals.

An example would spare me days of despair. Who has time and feels
like to translate the small function below to Julia?

def B(n):
    var('x')
    A, R = [], []
    for m in range(0, n, 1):
        A.append((-x)^m)
        for j in range(m, 0, -1):
            A[j - 1] = j * (A[j - 1] - A[j])
        S = sum((-1)^k*A[0].coefficient(x, k)/(k + 1) for k in (0..m))
        R.append(S)
    return R 

print B(13)

#2

You’ve tried SymEngine.jl or SymPy.jl?


#3

One issue you might have translating this to Julia is that it seems you are using Python’s negative array indexing (your inner loop at the first iteration will be A[-1] = 0 which makes sense in Python as it will wrap around, but will give an error in Julia, which doesn’t have circular indexing behavior. Outside of that using SymEngine or SymPy the code should be as tricky as changing the syntax.


#4

No, intentionally not. I would like to implement this very simple small
function without external tools (without modules, and without Python in
the background) to see how this could be done in Julia from the scratch.

Or is it really so difficult?


#5

Julia isn’t a CAS, so it doesn’t do symbolic calculations by default. You will have to use a library of some kind if you want to do symbolic calculations. Just like Python, Matlab, C++ do. The main ones in Julia are the ones mentioned by @ChrisRackauckas.


#6

SageMath is essentially Python + some external libraries for symbolic and numerical mathematics. I am not sure why you’d allow Python to use libraries here but not Julia… of course no programming language (except something like Mathematica which isn’t general purpose at all) won’t have this built in.


#7

It looks like this is just a polynomial, in which case you can use the Polynomials.jl package. If you use rational coefficients (possibly with BigInt) the result will be exact.

If you want to do it “without any packages” (but why would you?) you can just manipulate the coefficients by hand.


#8

In the end, much (but by no means all) of symbolic manipulation seems to be manipulating polynomials, for which you don’t actually need a real symbolic manipulation engine.


#9

There’s also Nemo, which may be overkill for this application: http://www.nemocas.org/


#10

Nemo definitely won’t be of any help here. There is no symbolic ring type in Nemo (though there may be a rudimentary one day).

What you require here is a symbolics package. Note that Sage isn’t quite Python and preprocesses some input so that it can behave more like a CAS, and behind the scenes it has some symbolics package which handles symbolics; it’s not baked into Python.

General purpose programming languages (like Julia and Python) are actually very poor for symbolics. You really need a specialised language to handle this properly, though you can get some of the way there with a specialised symbolics package. The really hard part is making it very fast.

Maple, for example, stores all its data in a DAG (Directed Acyclic Graph) format, which is completely different to the way general purpose languages handle variables and bindings.

The interesting question is whether it might be possible to introduce a special REPL mode in Julia one day that allows symbolics. I haven’t looked to see if anyone has already worked on such a thing.

If you are really asking how to implement a symbolic ring in Julia, you first need to decide how symbolic expressions will be stored (what data format), and how complex you want the supported expressions to be. Your simple example just requires polynomials, which means you just need terms and factors, essentially. But sophisticated systems like Maple support extremely complex expressions, along with very complex algorithms for expanding and simplifying expressions.

By the way, SymEngine is supposed to be a C++ library. So you could interface to that from Julia if you wanted to (using Cxx.jl). This would be the most logical way to support symbolic computation in Julia, unless you have a team of dozens of coders who are willing to write a Julia-only symbolic engine for you!


#11

Thank you for this excellent answer!

… This would be the most logical way to support symbolic
computation in Julia, unless you have a team of dozens of
coders who are willing to write a Julia-only symbolic engine for you!

Perhaps some day some smart guys will announce:
“Romeo: A Fresh Approach to Symbolic Computing”.

I’m afraid I’ll have to wait until then.


#12

For reference to others reading this thread: there’s already a SymEngine.jl package, as @ChrisRackauckas pointed out above (here’s the announcement post). It looks like @leiteiro was looking specifically for a Julia-only solution, though.


#13

Actually, I don’t think anyone has mentioned Symata.jl in this tread yet. Its development is currently halted until the Julia API stabilizes in 1.0, but perhaps this is more akin to what you’re looking for?

EDIT: nevermind. From the README:

Symata depends on the PyCall package and the Python SymPy module.


#14

… it’s not baked into Python.

It looks like @leiteiro was looking specifically for a Julia-only solution.

Right. “SymPy” is a pure Python library for symbolic mathematics.
So what I am looking for is an analogue package, pure Julia,
say “SymJulia”. Since SymPy is BSD licensed it should be
allowed to be inspired by its code.


#15

Many of you may not be aware that Reduce.jl allows you to use the full REDUCE computer algebra system within julia, in addition it can symbolically manipulate julia Expr objects.

It looks like you are just trying to compute the Bernoulli numbers. This is a bulit-in feature of Reduce.jl

If you use Pkg.checkout("Reduce") to get the master branch of my package right now, you can get this feature and use it like this:

julia> [bernoulli(RExpr("$i")) for i in 0:13]
14-element Array{Reduce.RExpr,1}:
 1            
 ( - 1)/2     
 1/6          
 0            
 ( - 1)/30    
 0            
 1/42         
 0            
 ( - 1)/30    
 0            
 5/66         
 0            
 ( - 691)/2730
 0       

This is the same exact result as you get in sage, but in Julia code. Hope that solves your issue.

print(B(13))
[1, 1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66, 0, -691/2730]

EDIT:
The syntax in Reduce.jl has been updated so it is more conveniet:

B(n::Int) = bernoulli.([0:n-1...])

#16

Many of you may not be aware that Reduce.jl allows you to use the full REDUCE computer
algebra system within julia, in addition it can symbolically manipulate julia Expr objects.

Thanks for pointing out Reduce. But the Bernoulli numbers can be implemented directly
in Julia in many ways, without the need for an additional package.

This is really about something else: I was looking for a basic approach for a rudimentary
implementation of symbolic computations in pure Julia.


#17

I’m actually working on a pure Julia package for some basic symbolic computations and simplification routines, with convenient conversion back into pure Julia functions. The current master hasn’t been updated in a while, and is actually probably broken, but if you go back a few commits, it should work quite consistently. I discovered some issues with the design, so I overhauled most of it, but I plan to push a new working version soon. https://github.com/eveydee/Sylvia.jl