July 20, 2016


After completing the work on gf_factor, I moved on to implement Gathen-Shoup’s factorization algorithm. Like Zassenhaus’s agorithm, it is also a probabilistic algorithm.
The paper is available here.

First question: Why this algorithm ?
> Because, it is kind of faster than zassenhaus’s algorithm.

Second question: What is “kind of” here ?
> Well, it is faster when a specific condition satisfies.

Its asymptotic runtime is O(n**2 + n log q).(log n)**2.loglog n, where n is degree of polynomial and q is the field characteristics.
In “Soft O” notation, it is O~(n**2 + n log q). While the cantor zassenhaus’s algorithm has O~(n**2 log q) asymptotic runtime.
So when log q approaches n, the difference is remarkable.

Algorithm Asymptotic Runtime
Shoup O~(n**2)
Zassenhaus O~(n**3)

It also works in three parts, firstly square free factorization, then distinct degree and finally equal degree factorization.
I have completed the implementation of this algorithm on shoup_factorization branch. And also, I changed the container (which stores factors) type to set.
Will send a PR soon.

July 18, 2016

Continue - Diophantine in Solveset :

PR 11234

  • For general Pythagorean diop_type (Diophantine eq type), it seems diophantine always returns parameterized solution so I did some changes in the PR. commit

  • You can refer this comment.

Continue Simplified Trig soln

PR #11188

  • After some changes, the PR is ready for review.

Continue nonlinsolve :

PR #11111

  • Added some XFAIL test-cases of system of Trigonometric equations. solveset trig solver (solve_trig) is not smart enough(solveset returns ConditionSet, where soln can be simply inverse trig functions using _invert or inverse Trigonometric functions). So solveset returns ConditionSet that means substitution is not getting soln.

  • It is better to replace trigonometric functions or other Function with symbols

(e.g. sin(x) –> u, sin(y)–> v, f(x)–> f_x, g(x) –> g_x)

and then solve for the symbols. After getting solution from nonlinsolve user can invert or do solveset

(e.g. solveset(Eq(sin(x), soln_u), x, domain) to get value of x).

  • Ready for review.

Meanwhile :

  • We already know that solveset need improved invert_real , invert_complex and Imageset Intersections.

  • previous work is in this PR 10971. Trying to improve them.

Some cases is here :

# In 2, 4, 5 intersection is not needed.

In [1]: img = ImageSet(Lambda(n, x/n), S.Complexes)

In [2]: Intersection(img, S.Complexes)
    ⎧x        ⎫
ℂ ∩ ⎨─ | n ∊ ℂ⎬
    ⎩n        ⎭

In [3]: img = ImageSet(Lambda(n, x/n), S.Integers)

In [4]: Intersection(img, S.Complexes)
⎧x        ⎫    
⎨─ | n ∊ ℤ⎬ ∩ ℂ
⎩n        ⎭    

In [5]: Intersection(ImageSet(Lambda(n, 2*n*I*pi), S.Integers), S.Complexes)
Out[5]: {2⋅ⅈ⋅π⋅n | n ∊ ℤ} ∩ ℂ

# ImageSet Intersection is not implemented when inverter returns multiple values.
# here ans should be {0, 1}
In [6]: img1 = ImageSet(Lambda(n, n**2), S.Integers)

In [7]: Intersection(img1, Interval(0,2))
         ⎧ 2        ⎫
[0, 2] ∩ ⎨n  | n ∊ ℤ⎬
         ⎩          ⎭


July 16, 2016

When this week started, most of the work up to week 7 were done. In the first couple of days this week, code for parsing was merged, and I started looking into wrapping lambdification of SymEngine Basic expressions.

This was a seemingly easy task, with two options ahead of me. One was to directly wrap C++’s lambdify_double method, while the other was to create Ruby lambdas directly.

After considering various aspects, and with feedback from other contributors, I decided to go ahead with writing them directly in Ruby.

It is now completed, and undergoing review.

The code implemented allows lamdifying SymEngine::Basic expressions with free symbols.

i.e. let, f = x + y, where x and y are SymEngine::Symbols,

a lambda is expected, such that it can be called, f_lambda.call(3, 4) which would return 7.

For expressions with single or multiple symbols, it works as following:

f = x + y
f_lambda = SymEngine::lambdify(f, [x, y])

For expressions with a single symbol, it can be used through the to_proc method, using a & call.

f_diam = 2 * SymEngine::PI * r
radii = [2, 3, 4]
diams = radii.map(&f_diam)


With another day left already, I am looking at Exception Handling, which I expect to be complicated in implementing.

See you next week!



July 15, 2016

Hello, guys. Welcome back. It’s been eight weeks into the coding period and this week I enjoyed a lot while working. I learned to build a list without using a loop, in case, it can be done more concisely with a list comprehension such as:

Instead of doing :

In [ ]: b = []
In [ ]: for x in a:
       b.append(10 * x)

we can do :

In [ ]: b = [10 * x for x in a]

I came to know about stuff , like using zip for transposing a matrix and dividing a list into groups of nThis week I had my weekly meeting with Jason and Sartaj on 13th of this month. They were in Austin, Texas attending Scipy 2016, the 15th annual Scientific Computing with Python conference.

So Far

  •  PR 11178  and PR 11237 got merged. So, now the master branch is updated with a module of singularity functions which can handle almost all the mathematical operations.
  • I have used the Singularity Functions  to continue developing the beam bending module at PR 11374. I have added the following methods:
    • load_as_SingularityFunction : This is a private method which represents a PointLoad and DistributedLoad object into a Singularity Function object and stores in the load beam object attribute.
    • load_distribution : This is a public method which outputs a load distribution curve of a beam.
    • shear_force : This is also a public method. It at first checks whether the moment list in boundary conditions dictionary have elements. If it has then it call the bending_moment method and then differentiate the output w.r.t to the free symbol else it just integrates the output of load_distribution method.
    • bending_moment : The initial part of this method is similar to the method shear_force. If the test passes then it integrates the load_distibution twice using contants of integrations and later the contants are solved else it just integrates the output of shear_force method.
    • slope : It outputs the slope of a beam by solving the constants.
    • deflection : It outputs the elastic curve of the beam.

Next Week

  • Modify the Beam Module.
  • Add Sphinx documentations.

That’s all for this week. Cheers !!

Happy Coding.

July 13, 2016


Previous week, I had started the work on Equal degree factorization here. Here I was storing the factors in a vector and then sorting the factors. Certik pointed out that its better to use some other data structure instead of vector. So I have changed the container type to set.
Then after both Distinct degree factorization and Equal degree factorization being implemented, I started to work on factoring a polynomial in finite field, this needed integrating square free factorization with these two. I worked on gf_factor() method. In this method we take a polynomial in a given field as input, and return all the factors and their respective powers, and polynomial’s leading coefficient as output.

GaloisFieldDict::gf_factor() const
  integer_class lc;
  GaloisFieldDict monic;
  gf_monic(lc, outArg(monic));
  if (monic.degree() < 1)
      return std::make_pair(lc, factors);
  std::vector<std::pair<GaloisFieldDict, integer_class>> sqf_list
      = monic.gf_sqf_list();
  for (auto a : sqf_list) {
      auto temp = (a.first).gf_zassenhaus();
      for (auto f : temp)
          factors.insert({f, a.second});
  return factors;

For a given polynomial firstly we get its monic representation and it leading coefficient. Then we find the Square free factors of the monic representation. And on each of the square free factor we run the zassenhaus’s algorithm.
I have been working on this branch, will create a PR after the Equal degree factorization PR gets merged.
Then I have started working on Shoup’s algorithm for polynomial factorization, will post about it in the coming weeks.

July 10, 2016

Hi everyone.

Here's what we have been doing for 7th week of GSoC.

Kalevi mentioned about the $LaTex$ not getting rendered on Planet Sympy website, thought it works fine on my website. Following the conversation Aaron opened issue planet-sympy/issues/45, though it hasn't been fixed.

This week we completed PR #11295 on Reidemeister Schreier algorithm. Remember the blog issue I asked on our gitter channel. Kalevi suggested to be more on the 'descriptive' side.

What do we do this week?

Implemented the Reidemeister Schreier algorithm (shortly RS algorithm) and Titze Transformations (shortly TT) in PR #11295 . Most of the pseudo code for RS algorithm is there in the Handbook [1], perhaps majority of the time was spent working with TT. There isn't any pseudo code available for making that work, but we gathered detailed information from [1] combined with Havas Paper [2]. What issues were there? In the both [1] and [2], there a few assumptions made, like
> shall assume that all relators are always cyclically reduced; that is, that
whenever a relator is changed, it is immediately replaced by its cyclic reduction.

I opened the issue #11352 regarding a typo I left in coset_enumeration_c. Though I didn't expect a PR, @kritkaran94 fixed it. Perhaps a good test case which makes use of different value of max_stack_size was suggested by Kalevi. One thing came back to me, "everthing is an object in Python", the fact that initially CosetTableDefaultMaxLimit was a module level variable initially and we changed it to CosetTable.CosetTableDefaultMaxLimit, module is an object so is a class.

What hasn't been done?

  • No testing for a few of techniques, for example currently no tests exist for elimination_technique_2, which is a variant of the elimination procedures.
  • TT doesn't produce the best possible result.

Though doing this seemed easy at first. But I didn't wanted to apply identity_cyclic_reduction every change as there can be a few instance where it not necessary to do this, perhaps because of the property of words.

The good thing about this week was that I could now understand the limitations that will remain after my GSoC project. The scope of Computational Group Theory is more than what I expected. Apart from that I will be leaving for my college this week stars from 15th of July, perhaps I don't know how things will shift regarding time for GSoC. Let's be realistic :)



  • 1. Derek F. Holt, Bettina Eick, Bettina, Eamonn A. O'Brien, "Handbook of computational group theory", Discrete Mathematics and its Applications (Boca Raton). Chapman & Hall/CRC, Boca Raton, FL, 2005. ISBN 1-5848-372-3 .
  • 2. George Havas, "Reidemeister-Schreier program"

I started working on the support for additional symbolic parameters in the module. So I added an extra argument in the method expr_to_holonomic for using custom domains. This will be helpful in integrations. For instance:

In [5]: expr_to_holonomic(sin(x*y), x=x, domain=QQ[y]).integrate(x).to_expr()
-cos(x⋅y) + 1
In [15]: expr_to_holonomic(log(x*y), x=x, domain=QQ[y]).integrate((x, 1, x)).to_expr()
Out[15]: x⋅log(x) + x⋅log(y) - x - log(y) + 1

Code upto here was merged at #11330.

After that I implemented the Frobenius method to support functions where the series may have negative or fractional exponents. Although there are limitations in the algorithm. First one being that we don’t have an algorithm to compute the general solution. There doesn’t seem to be a direct algorithm if roots of the indicial equation differ by an integer. Second one is the initial conditions. Initial conditions of the function (most of the times they won’t even exist) can’t be used if the exponents are fractional or negative.

In [23]: expr_to_holonomic(sqrt(x**2-x)).series()
            3/2       5/2       7/2         9/2         11/2
        C₀⋅x      C₀⋅x      C₀⋅x      5⋅C₀⋅x      7⋅C₀⋅x        ⎛ 6⎞
C₀⋅√x - ─────── - ─────── - ─────── - ───────── - ────────── + O⎝x ⎠
           2         8         16        128         256            

In [28]: expr_to_holonomic(cos(x)**2/x**2, initcond=False).series()
         2         4                       3         5
     C₂⋅x    2⋅C₂⋅x    C₁   2⋅C₁⋅x   2⋅C₁⋅x    4⋅C₁⋅x    C₀    ⎛ 6⎞
C₂ - ───── + ─────── + ── - ────── + ─────── - ─────── + ── + O⎝x ⎠
       3        45     x      3         15       315      2

I will be working more on this and hope to add more functionality in this method.

Right now I am working on an algorithm to convert a given Holonomic Function to Meijer G-function. This also doesn’t look so straightforward. Kalevi helped me with the theory and we wrote a basic that works for some cases. We hope to add more things in it.

In [17]: hyperexpand(expr_to_holonomic(exp(x)).to_meijerg())
In [19]: hyperexpand(expr_to_holonomic(log(x)).to_meijerg()).simplify()
Out[19]: log(x)

These are at #11360.

July 09, 2016

Hi there! It’s been seven weeks into the coding period and the second half has started now. I had a meeting with Jason on 3rd of July. We discussed some modifications that were needed to be done in the existing PRs. Further, we had a conversation on the beam bending module. We ended up with some good ideas regarding the implementation of beam object.

So Far

  • In PR  PR 11178  and PR 11237, I have added the suggested modifications.
  • In PR  11266, I have done some modifications in the DistributedLoad class. Initially, it had “start”,” end” and “value” as its attribute but Jason suggested me to add an attribute which would denote the order of the load. This attribute would help while representing a load in the form of Singularity Functions. Such as:
    • Order = 0 will denote Step Function,
    • Order = 1 will denote Ramp Function,
    • Order = 2 will denote Parabolic Ramp Function and
    • so on …
  • In the same PR, I have added some more method for taking boundary conditions explicitly for different cases such as for deflection, slope and moment as inputs. I made a private attribute _boundary_conditions and initiated as an empty dictionary of lists with keywords  deflection, slope and moment. Later on each of the methods, namely,
    • apply_boundary_conditions
    • apply_moment_boundary_conditions
    • apply_slope_boundary_conditions
    • apply_deflection_boundary_conditions

    just updates that dictionary. The method boundary_conditions would return the dictionary itself.

Next Week

  • I will try to get  PR 11178  and PR 11237 merged.
  • Add apply_load method to the Beam class.
  • Successfully convert all the load inputs through the  apply_load method into Singularity Functions.
  • Write methods which would output load curve, shear curve and moment curve.

That’s all for this week. Cheers !!

Happy Coding.

Week 6 was quite a rush, with many other things happening, and the work getting lagged behind. But with week 7, I was able to put much extra time into the project and get more or less up-to date with the time-line.

Since, the end of week 5, the progress is as follows:

  • Completed and merged Matrix CWrappers in PR #992.
  • Completed and merged Parser CWrapper in PR #1029.
  • Completed and merged PR #1027 which had wrappers for various functions for Number class, needed for the tests in Ruby to function properly.
  • PR #56 in the Ruby Wrapper for Matrices is almost ready, with several review feedback to be incorporated.
  • PR #60 in the Ruby Wrapper for Parser is ready to be merged.
  • PR #59 contains the updates for the IRuby Notebooks, which will be used by @certik for SciPy 2016’s SymEngine talk.
  • PR #55 with the Ruby wrappers for evalf was finally merged after quite a time going into accmmodating floating point tests.
  • Issue #57 was closed in the Ruby Wrapper.

In the proposal, I had promised an auxiliary feature of LaTeX friendly outputs for SymEngine Expresions, given I am on time. But considering I am barely hanging on to the time line, I decided to push this to the end of the time line.

So I will immediately start working on the lambdify functions, and hopefully start exception catching in the same week. This will give me time to come back to the LaTeX friendly expressions in the last weeks.

See you!


This last week was divided into two parts. The first half was spent in fixing miscellaneous issue with the last PR, regarding the basic to polynomial converisons. The second half was spent in thinking of how SymEngine should store rational polynomials (polynomials with rational coefficients) and how infact should they be implemented.

Leftovers from Basic Conversions

As mentioned, the first half of my week was spent finalizing the basic to polynomial conversions PR.

  • Using divnum and mulnum methods instead of using rcp_cast<Number> on mul and div respectively. These made the code much more readable too.

  • Use static_cast<const Rational &>(*rat_rcp).get_den() instead of rcp_static_cast<const Rational>(rat_rcp)->get_den(). The logic behind this is, that in the second case, another RCP is being constructed, and the get_den function is called on this new RCP. This is unecessary and can be avoided using the former approach.

  • At many places within the visitor pattern, I was calling the external function basic_to_poly. This was not wrong, but caused necessary overhead and recursion, which is not usually desired. Isuru suggested that I can reuse the visitor and successfully avoid the recursion. Here’s an example

dict_type _basic_to_poly(basic, gen)
    BasicToPolyV v;
    v.apply(basic, gen);

// dict is the answer we will be returning
void BasicToPolyV::bvisit(const Add &x)
    for (auto const &it : x.dict_)
        // unnecessary overhead
        dict += _basic_to_upoly(it, gen);

becomes this

// dict is the answer we will be returning
void BasicToPolyV::bvisit(const Add &x)
    dict_type res;
    for (auto const &it : x.dict_)
        // no overhead
        res += BasicToPolyV::apply(it, gen);
    dict = std::move(res);
  • There was also a discussion on how the visitor for Pow in find_gen_pow was implemented. Isuru pointed out, how it would give wrong answers in the cases he put forward. But I had thought of other cases (like 2**(2**(x+1))) which would would then not work, if the code was changed. These were special type of cases where the expression could be ‘simplified’. There’s some discussion on the issue here.

Starting Rational Polynomials

I began thinking of how the polynomials with rational coefficients class should be implemented. In Flint, the fmpq_poly class is implemented as a fmpz_poly class with another just another integer acting as the common denominatior. This kind of implementation saves up on a lot of time! (think of canonicalization of coefficients, after multiplication) But this kind of implementation would require a lot of extra code be written, and we could not reuse any of the code we have already written. So, after discussion with Isuru, he suggested that I can go ahead with the trivial implementation of using a map from uint to rational_class.

class URatDict : ODictWrapper<uint, rational_class>

Add in some constructors and that’s it! The internal container for rational polynomials in SymEngine is ready, because of the way we had strcutured the earlier code.

Now the question was how the end classes should be implemented (i.e. URatPoly) If you think about it, the methods of both the rational and integer polynomials will look exactly the same. from_vec, container_from_dict, eval, compare will all look the same with minor differences. The minor differences are only what type of container is being used (URatDict or UIntDict?) and what coefficient class is being used (integer_class or rational_class?) Infact, I checked this for our other polynomial types Flint and Piranha and it holds true over there too!

So, what came to mind was that all the methods had to be reused! Thus I formed a USymEnginePoly class which acts as the ‘SymEngine’ implementation of polynomials. You can pass it a BaseType from which it will deduce what type of polynomial would be constructed. Ideally, all you need to do now for implenting end polynomial classes is choose a implementation type (SymEngine or Piranha or Flint?) and a polynomial type (Integer or Rational?)

Here’s some pseudo code for the class structure aimed at

// the base
template <typename Container>
class UPolyBase : Basic

// super class for all non-expr polys, all methods which are
// common for all non-expr polys go here eg. degree, eval etc.
template <typename Cont, typename C>
class UNonExprPoly : UPolyBase<Cont>

// the specialized non-expr classes
template <typename D>
class UIntPolyBase : UNonExprPoly<D, integer_class>
// and
template <typename D>
class URatPolyBase : UNonExprPoly<D, rational_class>

// a specific implementation
// similar classes like UPiranhaPoly/UFlintPoly
template <template <typename X> class BaseType, typename D>
class USymEnginePoly : BaseType<D>

// end classes (few examples)
class UIntPoly : USymEnginePoly<UIntDict, UIntPolyBase>
class URatPoly : USymEnginePoly<URatDict, URatPolyBase>
class UIntPolyPir : UPiranhaPoly<piranha::poly, UIntPolyBase>
class URatPolyFlint : UFlintPoly<flint::poly, URatPolyBase>

Here’s a bonus picture!

I hope to get all the rational polynomials ready and tested in the next two to three days. The existing work can be seen in #1028


July 08, 2016

This week the focus was on the support code for the featherstone method. I finished adding examples to the docstrings of every function I made. I then wrote up test code for all of the new functions primarily focusing on expected outputs but included some expected error messages for one of the functions. Lastly I have coded up the functions themselves. This work can be seen in PR #11331.

I continued following the PR I reviewed last week and give suggestions as he worked on it. The PR is now in my opinion ready to be merged and is a beneficial addition to the sympy codebase.

The last thing I did this week was have a meeting about the presentation that I will be aiding in on Monday. After the meeting I have spent some time looking over my portions of the presentation and making sure I am prepared to speak.

Future Directions

Next week is the SciPy conference were I will be aiding in PyDy’s tutorial. Also I will be meeting with my mentor in person and during our time there I suspect we will work on a variety of things from Featherstone’s articulated body method to the base class.

PR’s and Issues

  • (Open) Added docstrings to ast.py PR #11333
  • (Open) [WIP] Featherstones EOM support PR #11331

July 07, 2016

This week I updated my PR#11277 to find the period of a general function.


In the past few weeks, I dedicated a lot of my time reading about the property of periodicity of a function. Earlier, I had implemented a trivial(and restricted) functionality for this task. This motivated me to study this topic as I planned to generalise the function.

Here are my notes on periodicity which were the literature reference for the development of the method:

  • Note that is a period of sin(x). But sin(x) has many other periods, such as , , and so on. However, sin(x) has no (positive) period shorter than .

  • If p is a period of f(x), and H is any function, then p is a period of H(f(x)).

  • For sums and products, the general situation is complicated.

    Let p be a period of f(x) and let q be a period of g(x). Suppose that there are positive integers a and b such that ap=bq=r.
    Then r is a period of f(x)+g(x), and also of f(x)g(x).
    However, the point to note here is that r need not be the shortest period of f(x)+g(x) or f(x)g(x).

    For example: The shortest period of sin(x) is , while the shortest period of (sinx)**2 is π.

    Another example: Let f(x)=sin(x), and g(x)=−sin(x). Each function has smallest period . But their sum is the 0-function, which has every positive number p as a period!

  • If p and q are periods of f(x) and g(x) respectively, then any common multiple of p and q is a period of H(f(x),g(x)) for any function H(u,v), in particular when H is addition, multiplication or division. However, it need not be the smallest period.

  • The sum of two periodic functions need not be periodic.

    For example: Let f(x)=sin(x)+cos(2πx). The function is not periodic.
    The problem is that 1 and are incommensurable. There do not exist positive integers a and b such that (a)(1)=(b)(2π).


I am abstracting the details of implementation so as not to make the post even further boring.

During the period of development, I faced few issues and had a lot of queries to make.

  1. The new implementation returns a value which might not be the fundamental period of the given function. The previous implementation, though limited, returned the fundamental period of the given function.

  2. The ability to find the LCM of irrationals. We will be dealing with the iconic π(and its multiples) in many of our cases(as is evident from the example above). Currently, we donot have the functionality to find the LCM of irrational numbers. A method needs to be developed to handle this issue.

  3. Issue with automatic simplification while verifying the result.

After Thoughts

I am looking forward to addressing all these issues in tonight’s meeting. Apart from that, implementing this was a lot of fun. I got to learn about inheritance and abstraction while implementing instance methods for periodic functions.

Hopefully, all my effort doesn’t go in vain.

Till next time !


After Distinct degree factorization being merged in this PR, I started working on Equal degree factorization.
Following up from previous example, we had:

x**5 + 10


x**10 + x**5 + 1

as distinct degree factors.

Now we run Equal degree factorization on both of the above given polynomial, from x**10 + x**5 + 1

We get:

x**2 + x + 1
x**2 + 3x + 9
x**2 + 4x + 5
x**2 + 5x + 3
x**2 + 9x + 4

And for:

x**5 + 10

We get:

x + 2
x + 6
x + 7
x + 8
x + 10

See that the first one gave degree two factors and the second one gave degree one factors.

I have implemented the algorithm in this PR.
Combining the distinct degree and equal degree factorization, we get the factors of a polynomial in finite field. For factorizing a polynomial in integral fielf we need to convert it to some finite field polynomial, factor it and then lift it back to integral field. This is Hensel’s Lifting, I will write about it in the coming blog posts.

I am sorry that I couldn’t work this week as I was at Robo Cup. Looking forward to a good week ahead.

July 04, 2016

ImageSet.put_values() :

PR #11343

  • After the discussion Harsh told that it is better to use like imageset.lamda(values_for_lambda_var) directly, also don’t make lambda variables public , it should be local and before doing imageset.lamda(values_for_lambda_var) this one need to check whether values are in ` base_set` or not.

  • You can see the previous code here : gist

  • I updated the docs stating how to put certain values in ImageSet lambda variables.(reverted my changes and edited the ImageSet docs in the PR.)

Continue nonlinsolve :

PR #11111

  • How Intersections and Complements are handled :

see these examples:

In [ ]: intr = Intersection(FiniteSet(x), Interval(1,10))

In [ ]: comp = Complement(FiniteSet(x), Interval(1,10))

In [ ]: intr
Out[ ]: [1, 10] ∩ {x}

In [ ]: comp
Out[ ]: {x} \ [1, 10]

In [ ]: type( Intersection(comp, Interval(1,11)))
Out[ ]: sympy.sets.sets.Complement

In [ ]: type(Complement(intr, Interval(1,2)))
Out[ ]: sympy.sets.sets.Complement

So first handling the Complements and then intersection will be checked in solveset soln.

  • nonlinsolve can handle simple trigonometric system of equations but when complex equations is used then it will return ConditionSet since solveset trig solver is not smart enough right now.

  • Trying to break the code into the functions to make the code better.

Continue Simplified Trig soln

PR #11188

  • I was getting problem in ImageSet union when you run test. Sometimes it pass all the cases but not always. I found the problem (most probably because order of args in union.reduce is not always same for all the python version I am using FiniteSet so I hope it is good way to handle this.) Now it passes all the checks all time.

Meanwhile :

  • I found that Mod is not defined for complex numbers. e.g.
In [1]: g = Mod(-log(3), 2*I*pi)

In [2]: g
Out[2]: Mod(-log(3), 2⋅ⅈ⋅π)

In [3]: simplify(g)
Out[3]: Mod(-log(3), 2⋅ⅈ⋅π)

In [4]: simplify(g)
Out[4]: Mod(-log(3), 2⋅ⅈ⋅π)

In [5]: g2 = Mod(-log(3), 2*pi)

In [6]: simplify(g2)
Out[6]: -log(3) + 2⋅π

In [7]: simplify(Mod(I,I))
Out[7]: 0

In [8]: simplify(Mod(2*I,I))
Out[8]: 0

In [9]: simplify(Mod(2*I,3*I))
Out[9]: 5⋅ⅈ

In [10]: simplify(Mod(2*I,1+3*I))
Out[10]: Mod(2⋅ⅈ, 1 + 3⋅ⅈ)

Need to implement Mod for complex number as well. There is concept og Gaussian Integers

Some resources I found is this : link1, link2, link3 link4

  • One can see the issue 11391 for detail discussion.

  • Tried to fix the bug of is_zero_dimensional in this PR #11371.


July 03, 2016


This week I continued my work from last week, but in a cleaner and structured fashion. The PR with last week’s work was closed in favour of a new PR, with better and scalable code. I’ll briefly explain how the functionality is implemented and the main ideas behind it. All the work is in #1003.

Finding Generators

You can refer to last week’s post to know what a generator for a polynomial is at a preliminary level. An expression may have any number of generators, but our job is to get the least number of generators which will be able to construct the given expression. The only rule is that a generator cannot be a Number (ie Rational or Integer). So sqrt(3) and pi can be generators but 2 or 1/2 can’t.

Initially, the approach I was trying out was not general (as I catered to only the Univariate Int case), and expanding it to other polynomials would prove to be immensely difficult. Isuru suggested that the function should return all the possible generators. This will be the most general case for this function, and we can adapt to specific cases based on how many / which generators the function returns. I will try and summarize how the function works at a high level.

If the expression passed to find_gen is a

  • Add and Mul : We can just call find_gen recursively on each of the expressions being added/multiplied. x + y or x*y will return both x and y as the generators. This follows from the fact that polynomials can be multiplied or added up to give the resultant polynomial.

  • Number : Do nothing, as numbers can never act as generators of polynomials.

  • Pow : Few cases arise in if the expression is a Pow. If the exponent is a positive integer, it suffices to find the generators of only the base of the expression. If it is a negative integer, we update the current generator set with base**(-1). For eg.

find_gen((x + y)**4) = find_gen(x + y) = (x, y)
find_gen((x**-2)) = (x**-1)

If the exponent is not an integer, the situation becomes a little complicated. It would seem intuitive that the generators in this case would be the base powered to the generators of the exponent. Like so,

find_gen(base**exp) = base**find_gen(exp)
find_gen(2**(x+y)) = 2**find_gen(x+y) = (2**x, 2**y)

This would seem to be working, but actually the same find_gen function cannot be used to get the generators of the exponent expression. The find_gen works in a different way once we are “inside” an exponent. Take for example :

// Incorrect
find_gen(2**(x + (1/2)))
= 2**find_gen(x + (1/2))
= 2**x

// Correct
find_gen(2**(x + (1/2)))
= 2**find_gen_pow(x + (1/2))
= (2**x, 2**(1/2))

find_gen_pow is another function which returns generators keeping in mind that the current expression it is dealing with is actually the exponent of another expression. So, it’s behaviour varies from the simple find_gen. Here is another example :

// Incorrect
find_gen(2**(x**2 + 3*x + (1/2))
= 2**find_gen(x**2 + 3*x + (1/2))
= 2**x

// Correct
find_gen(2**(x**2 + 3*x + (1/2))
= 2**find_gen_pow((x**2 + 3*x + (1/2))
= (2**x, 2**(x**2), 2**(1/2))
  • Basic : For all other structures like Symbol and Function, they themselves act as generators of the polynomial. We just need to update the generator set.

It is to be kept in mind that whenever we obtain a new potential generator, we update the current generator set. This may lead to modification of an already existing generator or add in a new one. This method thus takes care of some cases, not done by SymPy. A similar updation rule is followed for find_gen_pow where this is done on the coefficients of each expression instead of their powers.

find_gen :

gen_set = (x**(1/2))

// addition
gen_set = (x**(1/2), z)

// modification
gen_set = (x**(1/6))

// no change
gen_set = (x**(1/2))

find_gen_pow :

// "find_gen_pow"
gen_set = (x/2)

// addition
gen_set = (x/2, x**(1/2))

// modification
gen_set = (x/6)

// no change
gen_set = (x/2)

I have a short description on the storage of these generators in the PR description and subsequent comments. Please do have a look!

Polynomial Conversions

The next half of the week’s work involved actually converting a Basic into a polynomial. The API is too look as follows

RCP<Basic> UIntPoly::from_basic(Basic x);
RCP<Basic> UIntPoly::from_basic(Basic x, Basic gen);

The user can either supply a generator or one will automatically be deduced and used. The functions will throw of exactly one generator is not found, as we are dealing with univariate polynomials. Initially, I planned to write two different functions for UIntPoly and UExprPoly conversions, but soon I realized that most of the code was going to be identical. So, I formed a template function base with most of the common code, while the specific visitors inherit from this template base class with special methods.

The actual conversion was not too difficult, it just had to be broken down into cases like the find_gen problem. I leave it to the reader to try and theorize how this can be done. For eg. Imagine an Add, you can call b2poly on each expression added and add up the resulting polynomials. The code right now is scalable, and if we add in more polynomial classes like URatPoly there will be no issues in adding them in.

Miscellaneous Work

  • Templatized pow_upoly so that it can be used by any general polynomial class. Also removed some redundant code in the ExpandVisitor related to said function. Can be seen in #1010

  • I was just testing out how the parser is working out with the basic to polynomial conversions. It is working very seamlessly, constructing polynomials has never been easier! It’s as simple as

s = "2*(x+1)**10 + 3*(x+2)**5";
poly = UIntPoly::from_basic(parse(s));

s = "((x+1)**5)*(x+2)*(2*x + 1)**3";
poly = UIntPoly::from_basic(parse(s));


July 02, 2016

This week was majorly focused on fixing issues. Thanks to Ondrej, Aaron and Kalevi for testing the module and pointing them out. Some of the issues required a new functionality to be added and some were just plain technical bugs.

I started by writing a method which allows the user to change the point x0 where the initial conditions are stored. Firstly it tries converting to SymPy and then convert back to holonomic i.e. from_sympy(self.to_sympy(), x0=b). If the process fails somewhere then the method numerically integrates the differential equation to the point b. This method was then used in adding support for initial conditions at different points in addition and multiplication. Issues are  #11292 and #11293.

After that I implemented differentiation for a holonomic function. The algorithm is described here. A limitation is that sometimes it’s not able to compute sufficient initial conditions (mainly if the point x0 is singular).

There was also a bug in to_sequence() method. Apparently, whenever there weren’t sufficient initial conditions for the recurrence relation,  to_sympy() would fail. So I changed it to make to_recurrence() return unknown symbols C_j , where C_j representing the coefficient of x^j in the power series. So now even if we don’t have enough initial conditions, to_sympy() will return an answer with arbitrary symbols in it. This is also useful in power series expansion.

Earlier the method to_sympy() would only work if initial conditions are stored at 0 which is a big limitation. Thanks to Kalevi for the solution we now can find the hyper() representation of a holonomic function for any point x0 .(of course only if it exists.) This is further used by to_sympy() to convert to expressions. There is also something interesting Kalevi said “The best points to search for a hypergeometric series are the regular singular points”. I observed this turned out to be true a lot of times. We should keep this in mind when using to_sympy().

Later I added functionality to convert an algebraic function of the form p^(m/n) , for a polynomial p . Thanks to Aaron and Ondrej for the solution. Also added custom Exception to use in the holonomic module. Some small bugs #11319#11316 and 11318 were also fixed.

Currently I am trying to make from_sympy() work if additional symbols are given in the expression and also support other types of fields (floats, rationals (default right now), integers  ). This involves extending the ground domain of the Polynomials used internally. Hopefully this should be done in the next couple of days.

The unmerged fixes are at #11330.


July 01, 2016

Hi there! It’s been six weeks into the coding period, and it marks the half of GSoC. The Midterm evaluations have been done now and I have passed the mid-term evaluations. I am very much thankful to my mentors. This week I was having some physical illness due to this reason I was not able to attend the meeting with my mentors. This affected my workflow too. I was not able to work for almost 3 days. I will try to compensate the time lag by working some extra hours in the upcoming days.

So Far

  • Until last week, I have completed the implementation of PointLoad and  DistributedLoad. This week I started implementing the class for the principle object of the Beam Bending problems. I have implemented the Beam class. The work is not completed yet. I have initiated the class with length, Young’s Modulus and Moment of Inertia. I have to discuss with the mentor about the mutability of this class.
  • BoundaryConditions :- It is a public method of the beam class which takes the boundary conditions of the beam bending problem as input. I found out that there are only three types of boundary conditions : Deflection, Shear and Moment. So, I made them input in the form of keyword arguments and for each key, the value would be a list of tuples where each of the tuples would contain location and value in this format : (location, value). This method returns the boundary conditions in the form of a python dictionary. For example:  Let a problem have boundary conditions as: Moment(x = 0) = 1, Deflection(x = 0) = 0, Deflection(x = 4) = 0. Then the API is:
  • >>> from sympy.physics.mechanics.beam import Beam
    >>> from sympy import Symbol
    >>> E = Symbol('E')
    >>> I = Symbol('I')
    >>> b = Beam(4, E, I)
    >>> b.BoundaryConditions(moment = [(0, 1)], deflection = [(0, 0), (4, 0)])
    {'deflection': [(0, 0), (4, 0)], 'moment': [(0, 1)]}

Next Week

  • I will try to get  PR 11178  and PR 11237 merged.
  • Modify the Beam class with new methods.

That’s all for this week. Cheers !!

Happy Coding.

The main theme of this week is Featherstone’s method. I have finished reading all of the text book that I feel I need to in order to finish my project. After reading I realize that I have been improper about addressing my project. Instead of saying I am introducing Featherstone’s method to SymPy, it would be more accurate to say that I am introducing one of Featherstone’s methods. The book introduced two equations of motion generation methods for open loop kinematic trees and one method for a closed loop kinematic tree (I stopped reading after chapter 8 and so there may have been even more methods). For my project I have decided to focus on the articulated body method of equation of motion generation for kinematic trees. This method is presented as being more efficient than the composite body method and the closed loop method seems rather complicated.

With this in mind I began digging deeper into the articulated body method and better learning how it works. With this mindset I went over the three passes that the method uses and looked for places where code would be needed that isn’t specifically part of the method. I compiled a list of these functions and have written docstrings and presented them in PR #11331. The support code for the method includes operations for spatial (6D) vectors and a function and library for extracing relevant joint information.

This week I reviewed PR #11333. The pull request adds docstrings to method that did not have them previously which is a big plus but the docstrings that were added were minimal and vague. I asked that the contributer add more information to the docstrings and he said he will get to it.

Future Directions

Next week I plan on furthering my work on the articulated body method. I hope to have the support functions completely written up and to begin writing the equation of motion generator itself. These plans may be set aside, however, as my most active mentor will be coming back from traveling next week and so work may resume on the base class.

PR’s and Issues

  • (Open) Added docstrings to ast.py PR #11333
  • (Open) [WIP] Featherstones EOM support PR #11331

This week I worked on the topic of singularities.


A singularity is in general a point at which a given mathematical object is not defined.


  • 1/x has a singularity at x = 0 as it seems to reach infinity.
  • |x| (Absolute) has a singularity at x = 0 since it is not differentiable at that point.
  • √x (Square root) has a singularity at x = 0 since it doesnot admit a tangent there.

In real analysis, singularities are either discontinuities or discontinuities of the derivative (sometimes also discontinuities of higher order derivatives).

A removable singularity of a function is a point at which the function is undefined, but it is possible to redefine the function at that point in such a way that the resulting function is regular in a neighbourhood of that point.

For instance,

f(z) = sin(z)/z has a removable singularity at z = 0.
This singularity can be removed by defining f(0) := 1, which is the limit of f as z tends to 0.

In SymPy, we try to automatically simplify some of the expressions before passing it for further processing.

For example,

  • x/x simplifies to 1.
  • (x-1)^2/(x-1) simplifies to x-1.

At times, this can lead to incorrect behaviour. Suppose, we are interested to find the singularities ( and/or domain ) of the function f(x) = x/x.

In []: f = x/x
In []: singularities(f, x)
Out[]: S.EmptySet # no singularities
In []: continuous_domain(f, x, S.Reals) # finding the valid domain of f on the real line
Out[]: S.Reals # should be (-oo, 0) U (0, oo)

After some discussion about this issue of automatic simplification, we came to the conclusion that removable singularities should be treated as part of the domain of the function. This will help in justifying the behaviour of expressions which are automatically simplified.

So, I studied about how to classify a singularity as a removable one. For this, we use the Reimann’s theorem:

Let f a function defined on the set D \{a}. The point a is a removable discontinuity if there exists a neighborhood of a on which f is bounded.


lim (z -a)*f(z)

I have implemented this theorem in the form of a function but I face certain issues with the values returned by solveset due to check_domain function (which checks the bound of a function at a particular value). This remains to be discussed in the meeting.

After Thoughts

I have lots to discuss with my mentors on these issues.
After merging the PR#11277, I think solve_decomposition method will be good to go.

June 30, 2016

First of all the good news, I have passed the mid-term evaluations. I received the feedback from Kalevi, "I have had some problems with the blog posts.". Though Kalevi, mentioned to not take the post seriously. Well, not to say, that is definitely true. I don't claim that I will be able to remedy this, but lets see what happens.

Kalevi, had sent me a mail, mentioning "next things to do", it included "Cosets in Permuations group" and "Finitely Presented Abelian groups". (its been quite sometime since that mail). It was more of my decision to implement the topic of "presentation of subroups". To be true, I particularly chose this topic, since that makes the output of out work done till now, to look beautiful :) . Earlier I had a constant doubts when I was making proposal about what it exactly means by "presentation of subgroups". So while I read the examples from the Handbook, it became clear to me.

As to the progress in the topic, here's what we have done, I opened the PR on sunday. This week we started with the implementation of subgroup presentations. There are two methods to compute "presentation of a subgroup", first being called "computing presentations on schreier generators" also called the Reidemeister Schreier procedure. Since we have already implemented the "Todd Coxeter" procedure, which is used as the first step in Reidemeister Schreier for coset enumeration. It uses a routine which is called "define_schreier_generators", Second being "computing presentation on the generators of subgroup". Well the pseudo code is there for both the two implementations. The code of

The important part of finding the "presentation of subgroup" being simplifying the "presentations", i.e reducing the three factors |X| (no. of subgroup generators), |R| (no. of subgroup relators) and l (total length of subgroup relators). Though there doesn't seem to any pseudo code available anywhere for the implementation (this makes it quite interesting to implement). But the paper by Havas [Reidemeister Schreier program] (http://staff.itee.uq.edu.au/havas/1974h.pdf) neatly explains the things.

Last 2 days were spent on trying to make "elimination_technique_1" work. Its still not fully functional, since it still returns some dependent relators as a part of the presentation. As for an example, currently Reidemeister Schreier procedure presentation works this way.

>>> F, x, y = free_group("x, y")
>>> f = FpGroup(F, [x**2*y**2, y**-1*x*y*x**-3])
>>> H = [x]
>>> reidemeister_presentation(f, H)
([x_1], [x_1**-4, x_1**-8])

Two things need to be fixed here. First being, removal of -ve powers of single syllable relators, $x_1^{-4}$ -> $x_1^4$, similarly for the other relator. Second thing being, since $x_1^4 = 0$ implies that $x_1^8 = 0$ (dependent relator) so the relator should be removed. Perhaps the issue with this is, where should the code for checking this dependence be introduced? Kalevi seems a bit busy these days with other project. That is not the issue. I believe there may be something more to the simplification in the "Handbook", (have been referring to [2] for the elimination of generators, relators and simplification of presentations).

Yesterday, I Thought of focussing on a new topic "simplification of presentations", the method is short descripted in [2], but that is sufficient enough. It starts by checking whether any relator is of the form $gen^n$, where $n \in \mathbb{N}$. If any such relators are found then all other relators are processed for strings in the $gen$ known order. All string length exceeding length $n/2$ are replaced by their shortest equivalent. Further, for even `n` strings gen$^{-n/2}$ replaced by gen$^{n/2}$. I have written its code, perhaps it only requires a few more tests to be added to it.

I think I will be busy with "presentation of subgroups", this week and next week as well.


  • 1. Derek F. Holt, Bettina Eick, Bettina, Eamonn A. O'Brien, "Handbook of computational group theory", Discrete Mathematics and its Applications (Boca Raton). Chapman & Hall/CRC, Boca Raton, FL, 2005. ISBN 1-5848-372-3 .
  • 2. George Havas, "Reidemeister-Schreier program"

Hi everyone. Here's a brief summary of 5th week of my GSoC.

In the 5th week we worked upon the implementation of Low Index Subgroups algorithm. The basic problem dealt with by this algorithm is:

find all subgroups H such that the index $|G:H|$ is at most a specified number $N > 0$.

The term "low" is used to indicate the impracticalness of solving the problem for large $N$ values. The Handbook [1] contains the necessary pseudo code for its implementation, in which it has been explained in a top-down approch instead of the usual bottom-up approach. Complexity of the algorithm: exponential.

Status update: completion of Low Index Subgroup algorithm. The main function for calling the algorithm is low_index_subgroups with signature (G, N, Y=[]), here $G$ is the group for which we want to find the subgroups for, $N$ being a positive integer specifying the upper bound on the index value of subgroup generated, and $Y$ is optional argument specifying.

This routine also (coset_enumeration_c to needed this) makes use of calling cyclic conjugates of cyclic reductions of a set of relators, so I made a method in CosetTable named conjugates, which does this task now. We have a list $S$, which stores the coset tables for the founds subgroups. The algorithm later calls a routine called descendant_subgroups. One of the things being that in low_index_subgroups the descendant_subgroups routine can be called on any of the generators of inverse of generator of group $G$, though this isn't clear in the Handbook, but I some intuition to it.

descendant_subgroups($C$, $\{R_{1x}^C \}$, $R_2$, $N$)

try_descendant($C$, $\{R_{1x}^C\}$, $R_2$, $N$, $\alpha$, $x$, $\beta$)

Here is an example of Low Index Subgroup algorithm in work, Yay :) !!

>>> from sympy.combinatorics.free_group import free_group
>>> from sympy.combinatorics.fp_groups import FpGroup, low_index_subgroups
>>> F, x, y = free_group("x, y")
>>> f = FpGroup(F, [x**2, y**3, (x*y)**4])
>>> L = low_index_subgroups(f, 10)
>>> for i in L:
...     print(i.table)

[[0, 0, 0, 0]]
[[0, 0, 1, 2], [1, 1, 2, 0], [3, 3, 0, 1], [2, 2, 3, 3]]
[[0, 0, 1, 2], [2, 2, 2, 0], [1, 1, 0, 1]]
[[0, 0, 1, 2], [3, 3, 2, 0], [4, 4, 0, 1], [1, 1, 5, 4], [2, 2, 3, 5], [5, 5, 4, 3]]
[[1, 1, 0, 0], [0, 0, 1, 1]]
[[1, 1, 0, 0], [0, 0, 2, 3], [4, 4, 3, 1], [5, 5, 1, 2], [2, 2, 5, 6], [3, 3, 6, 4], [7, 7, 4, 5], [6, 6, 7, 7]]
[[1, 1, 1, 2], [0, 0, 2, 0], [3, 3, 0, 1], [2, 2, 4, 5], [5, 5, 5, 3], [4, 4, 3, 4]]
[[1, 1, 2, 3], [0, 0, 4, 5], [5, 5, 3, 0], [4, 4, 0, 2], [3, 3, 5, 1], [2, 2, 1, 4]]

Things I would love to have for this algorithm is to give, a presentation in the form Magma gives the subgroup returned in the form of generators of subgroup. That is coding the Reidemeister Schreier and its variant algorithms.

June 27, 2016


This week, I had been working on distinct degree factorization in this PR.

Distinct degree factorization splits a square-free polynomial into a product of polynomials whose irreducible factors all have the same degree.
Suppose we have a polynomial:

x**15 - 1

When we do distinct degree factorization of it, we get:

x**5 + 10
x**10 + x**5 + 1

Both the results consist of product of some polynomials with equal degree.

x**5 + 10

It is product of polynomials of degree 1, while:

x**10 + x**5 + 1

is product of polynomials of degree 2.

The factorization is achieved by exploiting the fact that:


So, for every degree i less than half of polynomial’s degree, we find the gcd of above computed polynomial and the given polynomial. If it’s not one, we add it to existing list of factors.

Then we have to find its equal degree factorization, I have implemented it and I am testing it with different cases. It will run on every input provided from distinct-degree factorization, and break it into polynomials with equal degrees.
Currently I am implementing Cantor-Zassenhaus’s algorithm.

June 26, 2016

We’re halfway there!

The past two weeks were supposed to be for implementing Matrices. But as I mentioned in the previous post, things got complicated with evalf methods, and I could not start working on the Matrices. This week PR #987 was merged up for evalf CWrappers. (The spec for the Ruby Wrapper needs some changes, which I have kept aside for now)

In SymEngine, matrices are of two forms, the DenseMatrix and the SparseMatrix. I have implemented the following constructors for them:

mat = SymEngine::DenseMatrix.new()
mat = SymEngine::DenseMatrix.new(SymEngine::DenseMatrix)
mat = SymEngine::DenseMatrix.new(Array)
mat = SymEngine::DenseMatrix.new(no_rows, no_columns)

mat = SymEngine::SparseMatrix.new()
mat = SymEngine::SparseMatrix.new(no_rows, no_columns)

Apart from that, for the DenseMatrix, the following methods were wrapped:

DenseMatrix.get(row, col)
DenseMatrix.set(row, col, SymEngine::Basic)
DenseMatrix.submatrix(row_start, col_start, row_end, col_end, row_step, col_step)
+ operator for Matrix and Scalar addition
* operator for Matrix and Scalar addition

Soon to be wrapped are:

DenseMatrix.LDLDenseMatrix.LU_solve # (to solve Ax = B type equations with LU decomposition)
DenseMatrix.FFLU # (fraction free LU)
DenseMatrix.FFLDU # (fraction free LDU)

The work can be views at PR #992 at symengine repo and PR #56 at symengine.rb repo.

This week was really nice for me in getting significant parts of the Matrices done,and with the luck continuing I will be able to wrap the necessary methods for SparseMatrix, get done with the tests and finish the Matrix part within the coming week, putting me back in track with the time-line.

Wish me luck!

Earlier this week, I continued my work on the method to convert a given Holonomic Function to Hypergeometric.

There were some major bugs that produced wrong results because of the fact that the recurrence relation doesn’t holds for all non-negative n. So I changed it to also store the value, say n0, so the recurrence holds for all n >=n0. I changed the method .to_hyper() accordingly and it works pretty fine now.

I also added the method to convert to symbolic expressions to_sympy() . It uses hyperexpand , however will only work for cases when a hyper representation is possible.

These were implemented at [#11246].

After that I have been majorly working on some issues Ondrej pointed out. These are #11292#11291#11287#11285#11284. The fixes are implemented at [#11297].

This week I will be fixing some more bugs, the ones listed at [Label: Holonomic Functions] and also a bug discussed with Kalevi on Gitter.


Continue nonlinsolve :

PR #11111

  • Removed the dependency of old solvers _invert. It was basically helpful for removing sqrt or Rational power. Now using unrad and little changes in substitution function can handle the system havning sqrt or cube root.

Continue - Diophantine in Solveset :

PR 11234

  • Some points we discussed :

    1. Solveset should return all the solution, so we need to check for which kind of eq. diophantine returns all the solution. We can get the equation type using classify_diop defined in diophantine.py.

    2. If there is the case where diophantine doesn’t return all the solutions then just return ConditionSet or try to make all solution (e.g. if just need to permute sign in values then use permute_signs).

  • After reading the diophantine Documentation and diophantine doctests and examples, I found that :

    1. Currently, following five types of Diophantine equations can be solved using ~sympy.solvers.diophantine.diophantine and other helper functions of the Diophantine module.
    • Linear Diophantine equations:

    a1*x1 + a2*x2 + ... + an*xn = b.

    • General binary quadratic equation:

    ax^2 + bxy + cy^2 + dx + ey + f = 0

    • Homogeneous ternary quadratic equation:

    ax^2 + by^2 + cz^2 + dxy + eyz + fzx = 0

    • Extended Pythagorean equation:

    a_{1} * x_{1}^2 + a_{2} * x_{2}^2 + .... + a_{n} * x_{n}^2 = a_{n+1} * x_{n+1}^2

    • General sum of squares:

    x_{1}^2 + x_{2}^2 + ... + x_{n}^2 = k

    1. If I am correct then Diophantine returns all the soln if eq is Linear Diophantine equations, General binary quadratic equation and Homogeneous ternary quadratic equation. I tried some Linear Diophantine equations and cross checked with online Diophantine solvers(one of the solver is here).

    2. In last 2 (Extended Pythagorean equation and General sum of squares) we need to do permute_signs to get all the soln.


          >>> from sympy.utilities.iterables import permute_signs
          >>> list(permute_signs((1, 12)))
          [(1, 12), (-1, 12), (1, -12), (-1, -12)]

    In general if all variables have even powers then we should do permute_signs. Diophantine solves the factors of the eq so if we have eq like this x**2 - y**2, diophantine solves factors x + y and x - y.

        >>> from sympy.solvers.diophantine import diophantine
        >>> from sympy.abc import x, y, z
        >>> diophantine(x**2 - y**2)
        set([(-t_0, -t_0), (t_0, -t_0)])

    so we should check even powers and permute sign.

    other examples :


    >>> from sympy.solvers.diophantine import diop_general_sum_of_squares
    >>> from sympy.abc import a, b, c, d, e, f
    >>> diop_general_sum_of_squares(a**2 + b**2 + c**2 + d**2 + e**2 - 2345)
    set([(15, 22, 22, 24, 24)])
    >>> from sympy.solvers.diophantine import diop_general_pythagorean
    >>> diop_general_pythagorean(a**2 + b**2 + c**2 - d**2)
    (m1**2 + m2**2 - m3**2, 2*m1*m3, 2*m2*m3, m1**2 + m2**2 + m3**2)
    >>> diop_general_pythagorean(9*a**2 - 4*b**2 + 16*c**2 + 25*d**2 + e**2)
    (10*m1**2  + 10*m2**2  + 10*m3**2 - 10*m4**2, 15*m1**2  + 15*m2**2  
      + 15*m3**2  + 15*m4**2, 15*m1*m4, 12*m2*m4, 60*m3*m4)
    >>> from sympy.solvers.diophantine import diop_general_sum_of_even_powers
    >>> diop_general_sum_of_even_powers(a**4 + b**4 - (2**4 + 3**4))
    set([(2, 3)])


    In above these types of cases we need permute_signs. If we check these solution you can see that permute_signs is needed when solutions is not parameterized.



      >>> from sympy.solvers.diophantine import diophantine
      >>> from sympy.abc import x, y, z
      >>> diophantine(x**2 - y**2)
      set([(-t_0, -t_0), (t_0, -t_0)])


    solution set([(-t_0, -t_0), (t_0, -t_0)]) is same as set([(-t_0, -t_0), (t_0, t_0), (-t_0, t_0), (t_0, -t_0)]).(because t_0 can take any integer value.)

    I discussed these things with @thilinarmtb (He have worked on diophantine. Blog link : https://thilinaatsympy.wordpress.com/). Main points are :

    Only the linear solver is incomplete, the algorithm in Diophantine module should be fixed (`permute_signs` or something like that won't help). We can use `permute_signs` when we have even powers in all variables. We should update Diophantine module to use permute sign. But we should not returns `ConditionSet` for linear diophantine eq. because for linear diophantine eq `diophantine()`returns parameterized solution which is complete most of the time.
    1. classify_diop can returns these diop_type :

      • linear
      • univariate
      • binary_quadratic
      • inhomogeneous_ternary_quadratic
      • homogeneous_ternary_quadratic_normal
      • homogeneous_ternary_quadratic
      • inhomogeneous_general_quadratic
      • inhomogeneous_general_quadratic
      • homogeneous_general_quadratic
      • general_sum_of_squares
      • general_pythagorean
      • cubic_thue
      • general_sum_of_even_powers

    If the equation type is none of these then solveset_integers should returns ConditionSet.Because currently diophantine can handle these kinds of eq only.

  • PR Update diophantine to get some missing solution: 11334

Continue simplified Trigonometric eq solution :

PR #11188

  • Shifted the _union_simplify function code to _union in fancyset/ImageSet with some changes.

  • Still facing problem (It sometimes pass all cases sometimes fail).

Meanwhile :

  • I found a problem in checksol defined in sympy/solvers/solvers.py and opened a PR that fixes the small issue, #11339



All of my week was mostly consumed by asking questions like “Should this throw?” and “What are it’s generators?”. My work for this week was to make a mechanism for converting arbitrary expression into polynomials. Turns out, it isn’t as trivial as it initially sounded. I will discuss my progress below. This week has been a lot more of logic involved, and thinking about how to go ahead instead of just development, so I do not have a lot to write about.

The Idea

A multivariate polynomial is made of generators. A univariate polynomial is made of a single generator. Each term in a polynomial is a monomial which is multiplication of it’s generators to a non-negative power, which also has a non-zero coefficient multiplied to it. Here are some polynomials and their generators :

x**2 + x**5          		->      x

sqrt(3) + sqrt(3)**3        ->		sqrt(3)

(1/x**2) + (1/x**5) 		-> 		(1/x)

x**2 + y**2 + x				-> 		x, y

2**x + 1					-> 		2**x

I hope you get the idea. The main task is given a Basic, we need to construct a polynomial with appropriate generator. The task is broken down into two parts. First off, we try and extract a generator from the given expression. Secondly, use the found (or user given) generator to construct the polynomial.

Ideally the API should look as follows:

RCP<Basic> UIntPoly::from_basic(Basic x);
RCP<Basic> UIntPoly::from_basic(Basic x, Basic gen);

My Progress

I have only worked on the univariate polynomial class as of now, and was to write the from_basic method just for the univariate classes for now. Although I just had to deal with the specific case of univariate integer polynomials, I haven’t successfully been able to get the job done just yet.

The work started with converting Symbol var to Basic var. Initially, the generators of the polynomial class could only be symbols, however this needed to be changed, as we’ve seen generators can be anything! This change did not take a lot of effort.

I started off with a find_gen function which given a Basic will try and find a unique generator for the polynomial to be constructed. Ideally, if no generator/ more than one generator is found it will throw a runtime error. I covered most of the cases but a few still remain.

Also, an initial implementation of the to_poly functionality has been done, complete with the required API. Some cases related to conversion still need to be handled. I do not want to go into the logic involved into each conversion as it is tedious and not very structured. All the work done till now can be seen in #998.

I hope with all the information I have gained over this week, I will wind up the Basic to UPoly conversions for both the integer coefficients as well as the general expression coefficients by the end of next week. I also intend to make the logic more concrete, so that the same API can be used in the multivariate case.

Miscellaneous Work

  • I made some changes in the expand function related to polynomials. They mostly have to deal with removal of code, which has already been implemented. It is also a partial fix for #886, work is in #999.


June 24, 2016

Hi there! It’s been five weeks and midterm evaluations are going on. I have completed the implementations of Singularity Function module and headed towards the Implementation of Beam Bending Module.

So far

  • In PR 11266 , I have implemented three major classes :
    • DistributedLoad : Have the attributes start, end and value. This object can be used to generate a distributed load. It is a mutable object. The start and  end attribute should be the instance of mechanics Point. The value can take any function as input. This Load can be represented as a Singularity Function (i have not implemented this functionality yet) and used for solving the beam bending problems.
    • PointLoad : There are two types : Moment and Point Load. Moments are basically the couple of forces and are represented by the Doublet function whereas Point Load is the force which acts on a single unique point of the Beam and is represented as a DiracDelta Function.
    • Beam : It is the principle object of the Beam Bending problems. I have started working on it.

Next Week

Thats all for this week. Cheers !!

Happy Coding.


Well I started this week off by getting sick and as such productivity took a little bit of a hit. The majority of this week was spent reading Featherstone’s text book. The example documentation showcasing the base class API still hasn’t been reviewed and so that part of the project will just have to be set aside until later. Overall the project will not suffer, however, because of my progress in learning Featherstone’s method.

I’ve done a couple of other things this week as well. In a discussion with Jason it was determined that the LagrangesMethod class must have access to the dynamic system’s bodies through the Lagrangian input. Upon research the Lagrangain turned out to not be an object instance but rather a function that simply returned the Lagrangian for input bodies. This meant that LagrangesMethod did not in fact have access to the dynamic system’s bodies. Due to this I decided that an easy way to get LagrangesMethod to have body information would be to add an optional keyword argument for it. This was LagrangesMethod can have a more similar API to KanesMethod. This change can be found in PR #11263.

This week I reviewed PR #10856 which claimed to fix Issue #10855. Upon review it seemed that the “fix” was to just not run tests that were failing. When researched it looks like a whole module has not been updated for Python 3.X and is failing its relative imports. When run in Python 2.X it’s still not working either but rather is throwing up many KeyError flags. I think this has not been caught sooner due to the module being a component directly dealing with another project (pyglet) thus the tests are not run by TravisCI.

Lastly there were some test errors in the example documentation for the base class on PyDy. I was not too worried about these because the PR is not currently awaiting merging and is simply a discussion PR. The failing tests, however, were not related to the changes in the PR and so a PyDy member submitted a PR that fixed the tests and asked me to review it. After I looked it over and determined that the fix addressed the issue correctly he merged the PR.

Future Directions

Next week I plan to continue forward with reading Featherstone’s book and, if possible, begin implementing one of the methods outlined in the book. Also I plan on beginning work on mirroring Jason’s overhaul of KanesMethod on LagrangesMethod.

PR’s and Issues

  • (Open) Added support for a bodies attribute to LagrangesMethod PR #11263
  • (Open) Added a depencency on older version of ipywidgets PR #100
  • (Open) Blacklisted pygletplot from doctests when pyglet is installed PR #10856
  • (Open) sympy.doctest(“plotting”) fails in python 3.5 Issue #10855
  • (Merged) Fix multiarray import error on appveyor PR #354

Hi folks !

The past couple of weeks were spent of developing heuristics for determining the fundamental period of a given trigonometric function.

In our higher school, we all must have come across Trigonometric Functions. One of the most striking properties of these functions is their periodicty.
The ability of a function to repeat its values in regular intervals has always caught my imagination.


Well, SymPy ought to have a functionality to determine the period of a function. The instigated me to implement this function now was the build failure in my PR#11141 on solve_decomposition.


function_range(sin(x), x, domain=S.Reals) causes the build to time-out as the number of critical points of sin(x) in the entire real domain is infinite. The same goes for other trigonometric functions as well.

However, if we can set the domain argument to be a finite interval which encompasses the entire behaviour of sin(x) over the entire real domain, our issue can be solved.

This led me to the idea of using the periodicity of the function as its domain.


f = sin(x) + sin(2*x) + cos(3*x)

We know that the period of f is 2⋅π i.e. the LCM of the periods of individual function.

It is known !

Hence, in order to find the period of f, we need the functionality to determine the period simpler trigonometric functions such as sin(x), sin(2*x) and cos(3*x)


If the period of g(x) is T, then the period of g(a*x) is T/a.

Using this property, we can easily compute the periods of sin(2*x) and cos(3*x) with our knowledge of the periodicity of the fundamental trigonometric functions.

June 21, 2016

Continue- nonlinsolve :

PR #11111

  • I found the case where complex solution is missed. When solve_poly_system is solving polynomial equation and that result is used in solving non polynomial-equation to get other soln. May be solve_poly_system return only real solution( so complex solution is not being used for further steps) thats why further step is with this real solution.

  • It seems solve_poly_system is doing it’s job since it is not designed for all solution. But we have substitution method using solveset_real and solveset_complex and retuning all solutionn. In new commits I improved the code and now substitution method can solve all kind of system independently.

  • I come across these kind of situation many times :

In [68]: x, y, z, r, t = symbols('x, y, z, r, t', real = True)

In [69]: soln = solveset(sqrt(r)*Abs(tan(t))/sqrt(tan(t)**2 + 1) + x*tan(t),x,S.Reals)

In [70]: soln
          ⎧      √r⋅│tan(t)│      ⎫
(-∞, ∞) ∩ ⎪───────────────────────⎪
          ⎨   _____________       ⎬
          ⎪  ╱    2               ⎪
          ⎩╲╱  tan (t) + 1 ⋅tan(t)⎭

In [71]: soln = solveset(sqrt(r)*Abs(tan(t))/sqrt(tan(t)**2 + 1) + x*tan(t),x)

In [72]: soln
⎧     -√r⋅│tan(t)│      ⎫
⎨   _____________       ⎬
⎪  ╱    2               ⎪
⎩╲╱  tan (t) + 1 ⋅tan(t)⎭

In other words :

In [1]: Intersection(FiniteSet(-x), S.Reals)
Out[1]: (-∞, ∞) ∩ {x}

In [2]: Intersection(FiniteSet(x), S.Reals)
Out[2]: ℝ ∩ {x}

So I am not able to extract the -x here. Because of this nonlinsolve returning some extra soln. One case is in test file : test_issue_5132.

I opened a PR for this #11280

previous PRs update :

  • #11188 - Simplified Solution for Trigonometric equation : I did some minor changes.Now I hope it is good to go.

  • #11234 - Connecting Diophantine and solveset to get Integer solution: added some testcases. I have defined solveset_integers for this, it take take list of symbols(to get soln in that order).t But Currently solveset takes one symbol. I hope this PR is good to go.

  • 11257 - solveset univariate trig inequality solvers : Actually to simplify the soln I need to use PR #11188, so if it got merged that code will help me in this PR.

Meanwhile :

  • I worked on these previous PRs regarding issues I reported #11239, 11280.


Older blog entries

Planet SymPy is made from the blogs of SymPy's contributors. The opinions it contains are those of the contributor. This site is powered by Rawdog and Rawdog RSS. Feed readers can read Planet SymPy with RSS, FOAF or OPML.