rburing's profile - activity

Looking at the source of desolve the difference is in how the variable tempic is set, namely the offending expression is (dvar==ics[1])._maxima_():

sage: y == e
e == y(t)

versus

sage: y == 1*e
y(t) == e

So it seems that Sage cannot be trusted to preserve the ordering LHS = RHS of an equation. This is definitely a bug which should be reported.

Please quote the respective lines which are causing the conflict (to make it easier for people to help you).

Maybe label the "unlabeled" by '' (the empty string) instead of None?

This is not my forte but it seems this labeling is a convention of Sage. Changing it would entail changing the implementation of BooleanFunction. Do you really want to change B and/or P? It is much easier to e.g. define P1 = P.subs(dict(zip(P.parent().gens(), reversed(P.parent().gens())))).

There is an option for MATLAB but not BASIC. However I imagine it's not too difficult to convert output in one of the supported languages into BASIC.

You're welcome. I wouldn't try to do that myself. I recommend (also in my answer, slightly edited since the original) to convert the expression to a sympy expression and to use their codegen which can generate C code.

Ah, I didn't think about the interfaces to other programs. What you can do is use your original variable names, and substitute them by the ones I defined at the very last step before output. Something like expr.subs({var('a_%d' % j) : SR.symbol('A[%d]' % j) for j in range(4)}).

Indeed var('A[0]') is not allowed, because it requires the name to be a valid Python identifier (probably because it wants to be able to inject it into the namespace). However we don't care about injecting into the namespace, because this is effectively already done through the name of the list.

Reading the source of var we find it calls SR.var which calls SR.symbol. This is the one we want:

sage: A = [SR.symbol('A[%d]' % j) for j in range(4)]
sage: A[3]^2
A[3]^2

Note that not all expressions are immediately valid C-expressions in this way, e.g. you still have to transform a^b to pow(a, b) or something similar. For this purpose you can e.g. convert the expression to a sympy expression with sympy.sympify(...) and use their C-code generator.

That is due to the implementation of fraction field elements. To render the expression in the desired way, see my answer.

My system has three equations and three unknowns, which I want to transform into a series of three 1-D root finding problems.

For this to be possible the system must have a finite number of solutions, but your system has an infinite number of solutions.

is there a test to see if my system of equations is pathological in some way?

You can compute I.dimension(). (By the way, this uses the Gröbner basis.) For 0-dimensional ideals in $\mathbb{Q}[x_1,\ldots,x_n]$ (like the Scholarpedia example) there are finitely many solutions over $\mathbb{C}$, and the output will be what you expect.

You can also call I.variety(ring=R) to solve the system of equations over R (e.g. CC or RR for approximate solutions, AA or QQbar for exact solutions), which will give an error if the ideal is not 0-dimensional.

This is a rational function, and an 'ordering of terms' is not part of the data, so Sage should not be expected to preserve it. Nevertheless the desired output can be obtained; see my answer.

The object you want to construct is a rational function (quotient of polynomials), and we can define it as such:

R.<x> = QQ[]
s = (x^2 + 2*x + 1) + 1/(x^2 + 2*x + 1)

We can obtain the two parts you want e.g. as follows (or you could have defined them individually):

sage: q, r = s.numerator().quo_rem(s.denominator())
sage: parts = (q, r/s.denominator())
sage: parts
(x^2 + 2*x + 1, 1/(x^2 + 2*x + 1))

Now we can generate the LaTeX for each part individually and join them by a plus symbol:

sage: latex_code = latex(parts[0]) + '+' + latex(parts[1])
sage: latex_code
x^{2} + 2 x + 1 + \frac{1}{x^{2} + 2 x + 1}
sage: show(latex_code)

Showing the LaTeX gives the desired output: $x^{2} + 2 x + 1 + \frac{1}{x^{2} + 2 x + 1}$.

Note: Buchberger's algorithm (to find a Gröbner basis) for linear forms is just Gaussian elimination.

My naive idea was to take all subchains of all maximal chains. This works, but it is inefficient because you get duplicates.

I would wrap the last command in sum(..., []) to get a single list of chains.

This is not my field so I can ask a silly question: Is every saturated chain a subchain of a maximal chain?

One approach is to use symbolic expressions instead of functions, e.g.

var('t')
r = vector([1, t])

The addition example becomes slightly verbose (but certainly unambiguous):

r(t=1) + r(t=2)

An apparent downside is that this blurs the line between variables and symbolic constants which might also appear:

var('t,c')
s = vector([1, c*t])

However, many operations which are really operations on functions are implemented in Sage as operations on symbolic expressions, where the dependent variables are passed alongside the symbolic expression.

For example, to calculate the Jacobian:

sage: jacobian(s, [t])
[0]
[c]

It seems Sage has trouble parsing the output from Maxima. Can you/someone try solving the same system in Maxima (with desolve I presume) and adding the output?

How about wrapping commands in sage_eval(...)? That should give you access to the preparser from Python.

Use Ric[0,0].expr() to convert Ric[0,0] to the desired symbolic expression:

desolve(Ric[0,0].expr(), h, ivar=t, contrib_ode=True)

Here Sage is letting PARI do the hard work, creating a "Buchmann's number field" (bnf), which is a data structure containing the defining data of the number field and invariants such as the class group structure. The section Class group, units, and the GRH of the page I linked earlier states that

Some of the functions starting with bnf are implementations of the sub-exponential algorithms for finding class and unit groups under GRH, due to Hafner-McCurley, Buchmann and Cohen-Diaz-Olivier.

Some more names are dropped further along the page.

In the source code of PARI one sees that bnfinit corresponds to the C function bnfinit0 which calls Buchall_param, which seems to be doing much of the work. The source file credits McCurley and Buchmann. At this point it is probably it is easier to find a reference than it is to read Buchall_param.

Note that Sage's parameter proof tells PARI to run bnfcertify. From the same section quoted before:

Warning. Make sure you understand the above! By default, most of the bnf routines depend on the correctness of the GRH. In particular, any of the class number, class group structure, class group generators, regulator and fundamental units may be wrong, independently of each other. Any result computed from such a bnf may be wrong. The only guarantee is that the units given generate a subgroup of finite index in the full unit group. You must use bnfcertify to certify the computations unconditionally.

The source for bnfcertify0 refers to Zimmert's bound $B$ and consists of two phases. The first phase for the class group calls bnftestprimes and according to the User's Guide

Running this function successfully proves that the classes of all prime ideals of norm $\leq B$ belong to the subgroup of the class group generated by the factorbase used to compute the bnf (equal to the class group under GRH). If the condition is not true, then (GRH is false and) the function will run forever.

The second phase tests the units and relations, but I am not sure how.

That is because you didn't fix the case l == 1, which I did in my answer. To be explicit:

yield alphabet

is not the same as

for a in alphabet:
    yield [a]

Note the original function words is recursive: it calls itself.

Your new function non_rep_words should also call itself, rather than calling the old words:

def non_rep_words(alphabet,l):
    if l == 0:
        yield []
    elif l == 1:
        for a in alphabet:
            yield [a]
    else:
        for word in non_rep_words(alphabet, l-1):
            for a in alphabet:
                if word[-1] != a:
                    yield word+[a]

Here I also corrected the case l == 1.

Note that you could avoid this case by adding an extra clause to your if statement:

def non_rep_words(alphabet,l):
    if l == 0:
        yield []
    else:
        for word in non_rep_words(alphabet, l-1):
            for a in alphabet:
                if len(word) == 0 or word[-1] != a:
                    yield word+[a]

However this will perform many pointless comparisons to 0, so it is probably slower.