That’s all for this week. Cheers !!
Happy Coding.
After completing the work on gf_factor
, I moved on to implement GathenShoup’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.
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
Continue nonlinsolve :
PR #11111
Added some XFAIL testcases 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
).
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)
Out[2]:
⎧x ⎫
ℂ ∩ ⎨─  n ∊ ℂ⎬
⎩n ⎭
In [3]: img = ImageSet(Lambda(n, x/n), S.Integers)
In [4]: Intersection(img, S.Complexes)
Out[4]:
⎧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))
Out[7]:
⎧ 2 ⎫
[0, 2] ∩ ⎨n  n ∊ ℤ⎬
⎩ ⎭
continue..
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!
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 n . This 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.
Next Week
That’s all for this week. Cheers !!
Happy Coding.
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.
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 planetsympy/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.
> 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
.
elimination_technique_2
, which is a variant of the elimination procedures.
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 :)
अलविदा
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() Out[5]: cos(x⋅y) + 1 ───────────── y 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**2x)).series() Out[23]: 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() Out[28]: 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 x
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 Gfunction. 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()) Out[17]: x ℯ In [19]: hyperexpand(expr_to_holonomic(log(x)).to_meijerg()).simplify() Out[19]: log(x)
These are at #11360.
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 3^{rd} 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.
just updates that dictionary. The method boundary_conditions would return the dictionary itself.
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 upto date with the timeline.
Since, the end of week 5, the progress is as follows:
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.
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);
}
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.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 nonexpr polys, all methods which are
// common for all nonexpr polys go here eg. degree, eval etc.
template <typename Cont, typename C>
class UNonExprPoly : UPolyBase<Cont>
// the specialized nonexpr 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
Tata!
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.
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.


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 2π
is a period of sin(x)
.
But sin(x)
has many other periods, such as 4π
, 6π
, and so on.
However, sin(x)
has no (positive) period shorter than 2π
.
If
p
is a period off(x)
, andH
is any function, thenp
is a period ofH(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 2π
, 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 2π
. But their sum is the 0
function, which has every positive number p
as a period!
If
p
andq
are periods off(x)
andg(x)
respectively, then any common multiple ofp
andq
is a period ofH(f(x),g(x))
for any functionH(u,v)
, in particular whenH
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 2π
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.
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.
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.
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
and
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.
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
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
Meanwhile :
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.
continue..
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.
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)
eg.
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
update_gen_set(z)
gen_set = (x**(1/2), z)
// modification
update_gen_set(x**(1/3))
gen_set = (x**(1/6))
// no change
update_gen_set(x)
gen_set = (x**(1/2))
find_gen_pow
:
// "find_gen_pow"
gen_set = (x/2)
// addition
update_gen_set(x**(1/2))
gen_set = (x/2, x**(1/2))
// modification
update_gen_set(x/3)
gen_set = (x/6)
// no change
update_gen_set(x)
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!
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.
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));
Laters!
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.
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 midterm 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
>>> 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
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.
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.
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.
Examples:
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
.(x1)^2/(x1)
simplifies to x1
.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 setD \{a}
. The pointa
is a removable discontinuity if there exists a neighborhood ofa
on whichf
is bounded.
Mathematically,
lim (z a)*f(z)
z>a+
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.
First of all the good news, I have passed the midterm 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.
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 topdown approch instead
of the usual bottomup 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.
This week, I had been working on distinct degree factorization in this PR.
Distinct degree factorization splits a squarefree 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 distinctdegree factorization, and break it into polynomials with equal degrees.
Currently I am implementing CantorZassenhaus’s algorithm.
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.inv DenseMatrix.det DenseMatrix.transpose DenseMatrix.submatrix(row_start, col_start, row_end, col_end, row_step, col_step) DenseMatrix.rows DenseMatrix.cols + operator for Matrix and Scalar addition * operator for Matrix and Scalar addition DenseMatrix.to_s DenseMatrix.inspect
Soon to be wrapped are:
DenseMatrix.LU 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 timeline.
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 nonnegative 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.
Cheers!
Continue nonlinsolve :
PR #11111
_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 :
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
.
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 :
~sympy.solvers.diophantine.diophantine
and other helper functions of the Diophantine module.a1*x1 + a2*x2 + ... + an*xn = b
.
ax^2 + bxy + cy^2 + dx + ey + f = 0
ax^2 + by^2 + cz^2 + dxy + eyz + fzx = 0
a_{1} * x_{1}^2 + a_{2} * x_{2}^2 + .... + a_{n} * x_{n}^2 = a_{n+1} * x_{n+1}^2
x_{1}^2 + x_{2}^2 + ... + x_{n}^2 = k
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).
In last 2 (Extended Pythagorean equation and General sum of squares) we need to do permute_signs
to get all the soln.
E.g.
```
>>> 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.
e.g.
```
>>> 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.
classify_diop
can returns these diop_type
:
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 :
checksol
defined in sympy/solvers/solvers.py
and opened a PR that fixes
the small issue, #11339continue…
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.
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 nonnegative power, which also has a nonzero 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);
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.
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.Tata!
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.
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.
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
.
 (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 timeout 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)
isT
, then the period ofg(a*x)
isT/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.
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 polynomialequation 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
Out[70]:
⎧ √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
Out[72]:
⎧ √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
nonlinsolve
then check these gist nonlinsolve_till_24june_2016,
nonlinsolve_till_21jun2016.py.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 :
continue..
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.