- GiNaC-list - ginac.de

Sqrfree_parfrac Issue GiNaC
by Konstantin Leyde 24 Jul '20

24 Jul '20

Hello, My name is Konstantin Leyde, and I am a master student in physics at the ETH Zurich. For my thesis project, I have been working with the GiNaC library, and I have encountered an issue with the partial fractioning (sqrfree_parfrac) routine: Unless I misunderstood the way the sqrfree_parfrac routine works, it appears that for special cases, the output of the sqrfree_parfrac routine misses a numerical prefactor for one term, and thus returns the wrong answer. Please find attached a pdf, that explains the problem in more detail. In this pdf (as well as separately) we have attached the code that shows the issue. Many thanks in advance, Best wishes, Konstantin

2 1

Release 1.7.10
by Richard B. Kreckel 28 Jun '20

28 Jun '20

Hi, GiNaC 1.7.10 has just been released. The changes are: * Fix collect_common_factors() for hidden zero arguments. * Fix build on modern systems (C++20 and Python3). Take care! -richy. -- Richard B. Kreckel <https://in.terlu.de/~kreckel/>

1 0

C binding for GiNaC
by Xianwen Chen (陈贤文) 22 Jun '20

22 Jun '20

Dear all, I wonder if it is possible to use GiNaC in a C program? Is there a binding for that? I guess there may be one, because GiNaC is used in other software, for example, Octave and Sage Math. But I couldn't find relevant information. Yours sincerely, Xianwen

2 2

problem about collect_common_factors
by Feng Feng 22 Jun '20

22 Jun '20

Hi all, I got a strange problem about collect_common_factors, the problem can be reproduced by the following example: ------ symbol a("a"),b("b"),c("c"); symtab table; table["a"] = a; table["b"] = b; table["c"] = c; parser reader(table); ex res = reader("b*((c+c*(1+4*a+2*b+b^2)-c*(1+2*b+b^2)-(1+4*a)*c))-a*(c+c*(1+4*a+2*b+b^2)-c*(1+2*b+b^2)-(1+4*a)*c)"); cout << res << endl; cout << collect_common_factors(res) << endl; ------ the collect_common_factors will throw an error: division by zero. so is this a bug for collect_common_factors? Best regards! Feng

3 2

Issue with derivative of "abs"
by Pierangelo Masarati 16 Jun '20

16 Jun '20

Dear GiNaC users, we have been using GiNaC for years now, to provide our multibody dynamics solver, MBDyn <http://www.mbdyn.org/>, some symbolic differentiation capabilities (perhaps not the most appropriate use, but some users seem to like it). I'll make the story as short as possible. In one of the most interesting applications, the user needs to: - define the variables that represent strains (and strain rates), - the expression for the configuration-dependent generalized forces in short, a constitutive law. The expression's partial derivatives with respect to the strains and strain rates are then automatically computed by GiNaC, and, when needed, the numerical values are substituted and the expressions are evaluated. The related code snippet looks like > GiNaC::symbol gEps("v"); // parameter symbols > GiNaC::ex gExpr; // expressions > GiNaC::ex gExprDEps; // derivatives > > GiNaC::lst l; > > l.append(gEps); > gExpr = GiNaC::ex("v*abs(v)", l); > > gExprDEps = gExpr.diff(gEps); One user recently tried something like this: var = v expr = v*abs(v) # a turbulent viscous law the result was dexpr/dvar = abs(v)+1/2*v*(v+conjugate(v))*abs(v)^(-1) which fails when v = 0 because of the division by zero. I tried with "realsymbol" (since our values can only be real); things, however, didn't change much: dexpr/dvar = v^2*abs(v)^(-1)+abs(v) I realized that the issue is in the derivative of the abs() function; indeed, the proposed formula is rather general (despite failing when numerically evaluated with 0 as the argument). When the argument is real, a more robust representation would be "d(abs(v))/dv = sign(v)". In the case at hand, "d(v*abs(v))/dv = 2*abs(v)" or, as a second choice, "d(v*abs(v))/dv = abs(v) + v*sign(v)". First, I considered the possibility to implement the "sign()" function in GiNaC, but I am not familiar enough with its internals. If it makes sense, I would definitely recommend its introduction. However, I noticed that the "step()" function is implemented (and step(0) correctly evaluates to 1/2); the "sign()" function can then be obtained as "sign(x) = 2*step(x) - 1" when "x" is real. To this end, I modified GiNaC's code as reported here: > diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp > index d68afbb5..ec74a435 100644 > --- a/ginac/inifcns.cpp > +++ b/ginac/inifcns.cpp > @@ -321,6 +321,8 @@ static ex abs_expand(const ex & arg, unsigned options) > static ex abs_expl_derivative(const ex & arg, const symbol & s) > { > ex diff_arg = arg.diff(s); > + if (arg.info(info_flags::real)) > + return diff_arg*(2*step(arg) - 1); > return > (diff_arg*arg.conjugate()+arg*diff_arg.conjugate())/2/abs(arg); > } which indeed gives the expected result: dexpr/dvar = v*(-1+2*step(v))+abs(v) I wonder whether my fix can be considered correct and, in case, whether it could be included in the distribution. Sincerely, p. -- Pierangelo Masarati Professore Ordinario di Costruzioni e Strutture Aerospaziali Dipartimento di Scienze e Tecnologie Aerospaziali Politecnico di Milano

2 1

Re: [GiNaC-list] Issue with derivative of "abs"
by Vladimir V. Kisil 16 Jun '20

16 Jun '20

Dear Pierangelo, Thanks for the email. It is nice to see another usage of GiNaC! Your patch combines the best from both worlds as was discussed: https://www.ginac.de/pipermail/ginac-devel/2014-April/002105.html So I am voting in favour of this patch to be applied! Meanwhile, it seems the faster solution for your would be to clone abs function with the alteration towards the desired derivative. It is not difficult: 1. Cut & Paste everything for definition abs from GiNaC's inifcns.h and inifcns.cpp to your code. 2. Replace the derivative function by abs_expl_derivative() from your patch. 3. Search & replace abs -> myabs in this portion of code. Then the function myabs() shall be ready to be used. Best wishes, Vladimir -- Vladimir V. Kisil http://www.maths.leeds.ac.uk/~kisilv/ Book: Geometry of Mobius Transformations http://goo.gl/EaG2Vu Software: Geometry of cycles http://moebinv.sourceforge.net/ Jupyter: https://github.com/vvkisil/MoebInv-notebooks >>>>> On Tue, 16 Jun 2020 19:24:32 +0200, Pierangelo Masarati <pierangelo.masar\ ati(a)polimi.it> said: PM> Dear GiNaC users, PM> we have been using GiNaC for years now, to provide our multibody PM> dynamics solver, MBDyn <http://www.mbdyn.org/>, some symbolic PM> differentiation capabilities (perhaps not the most appropriate PM> use, but some users seem to like it). PM> I'll make the story as short as possible. In one of the most PM> interesting applications, the user needs to: PM> - define the variables that represent strains (and strain PM> rates), PM> - the expression for the configuration-dependent generalized PM> forces PM> in short, a constitutive law. The expression's partial PM> derivatives with respect to the strains and strain rates are PM> then automatically computed by GiNaC, and, when needed, the PM> numerical values are substituted and the expressions are PM> evaluated. PM> The related code snippet looks like >> GiNaC::symbol gEps("v"); // parameter symbols >> GiNaC::ex gExpr; // expressions GiNaC::ex gExprDEps; // >> derivatives >> >> GiNaC::lst l; >> >> l.append(gEps); gExpr = GiNaC::ex("v*abs(v)", l); >> >> gExprDEps = gExpr.diff(gEps); PM> One user recently tried something like this: PM> var = v PM> expr = v*abs(v) # a turbulent viscous law PM> the result was PM> dexpr/dvar = abs(v)+1/2*v*(v+conjugate(v))*abs(v)^(-1) PM> which fails when v = 0 because of the division by zero. I tried PM> with "realsymbol" (since our values can only be real); things, PM> however, didn't change much: PM> dexpr/dvar = v^2*abs(v)^(-1)+abs(v) PM> I realized that the issue is in the derivative of the abs() PM> function; indeed, the proposed formula is rather general PM> (despite failing when numerically evaluated with 0 as the PM> argument). When the argument is real, a more robust PM> representation would be "d(abs(v))/dv = sign(v)". In the case PM> at hand, "d(v*abs(v))/dv = 2*abs(v)" or, as a second choice, PM> "d(v*abs(v))/dv = abs(v) + v*sign(v)". PM> First, I considered the possibility to implement the "sign()" PM> function in GiNaC, but I am not familiar enough with its PM> internals. If it makes sense, I would definitely recommend its PM> introduction. PM> However, I noticed that the "step()" function is implemented PM> (and step(0) correctly evaluates to 1/2); the "sign()" function PM> can then be obtained as "sign(x) = 2*step(x) - 1" when "x" is PM> real. To this end, I modified GiNaC's code as reported here: >> diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index >> d68afbb5..ec74a435 100644 --- a/ginac/inifcns.cpp +++ >> b/ginac/inifcns.cpp @@ -321,6 +321,8 @@ static ex >> abs_expand(const ex & arg, unsigned options) static ex >> abs_expl_derivative(const ex & arg, const symbol & s) { >> ex diff_arg = arg.diff(s); + if >> (arg.info(info_flags::real)) + return >> diff_arg*(2*step(arg) - 1); return >> (diff_arg*arg.conjugate()+arg*diff_arg.conjugate())/2/abs(arg); >> } PM> which indeed gives the expected result: PM> dexpr/dvar = v*(-1+2*step(v))+abs(v) PM> I wonder whether my fix can be considered correct and, in case, PM> whether it could be included in the distribution. PM> Sincerely, p. PM> -- Pierangelo Masarati Professore Ordinario di Costruzioni e PM> Strutture Aerospaziali Dipartimento di Scienze e Tecnologie PM> Aerospaziali Politecnico di Milano PM> _______________________________________________ GiNaC-list PM> mailing list GiNaC-list(a)ginac.de PM> https://www.ginac.de/mailman/listinfo/ginac-list

1 0

Release 1.7.9
by Richard B. Kreckel 11 Apr '20

11 Apr '20

Hi, GiNaC 1.7.9 has just been released. The changes are: * Fix unarchiving of overloaded functions. * Fix MinGW64 build. Happy Easter! -richy. -- Richard B. Kreckel <https://in.terlu.de/~kreckel/>

1 0

archive/unarchive on functions with a variable number of arguments
by Feng Feng 08 Apr '20

08 Apr '20

Hi all, I came across a problem about the archive and unarchive on functions with a variable number of arguments. for example, the following c++ code produces the output: zeta(10,x) :> 1/93555*Pi^10. it seems that the unarchive just return zeta(10). Best regards! Feng //--- symbol x("x"); ex zx = zeta(10,x); archive ar_out; ar_out.archive_ex(zx, "zx"); ofstream out("tmp.gar"); out << ar_out; out.close(); lst syms = {x}; archive ar_in; ifstream in("tmp.gar"); in >> ar_in; in.close(); cout << zx << " :> " << ar_in.unarchive_ex(syms,"zx") << endl; //---

3 4

Re: [GiNaC-list] series of x^n around x=0
by Feng Feng 23 Dec '19

23 Dec '19

Of course, the issue here is mathematics, not GiNaC. I totally agree with that. :) And thanks a lot for pointing out the method to find the leading power. Since there is also log(x)-terms in my expression, and leading power seems not be able to capture these log terms. Currently, I am trying to walk around by finding the all pattern: x^n, since n is rational number, and multiply n by lcm of all those denominators, and transform all rational powers to integer powers … BTW, when the input contains x^n (n as an integer) and log(x), it seems that series works fine around x=0. Thanks again! Best regards! Feng 原始邮件发件人: Vladimir V. Kisil<kisilv(a)maths.leeds.ac.uk> 收件人: Feng Feng<F.Feng(a)outlook.com> 抄送: GiNaC discussion list<ginac-list(a)ginac.de>; Vladimir V. Kisil<V.Kisil(a)leeds.ac.uk> 发送时间: 2019年12月22日(周日) 06:17 主题: Re: [GiNaC-list] series of x^n around x=0 >>>>> On Sat, 21 Dec 2019 11:14:58 +0000, Feng Feng <F.Feng(a)outlook.com<mailto:F.Feng@outlook.com>> said: FF> Thanks very much for the reply.! And yes, it is not FF> well-defined in mathematical sense to taylor expand x^(3/2) FF> around x=0. What I exactly want is to get the asyptotic FF> behaviour of a function f(x) around x=0. For example, f(x) = FF> x^(3/2) (1 + x + x^2 + …), the asyptotic behaviour at leading FF> order (LO) and next-to-leading order (NLO) are x^(3/2) and FF> x^(5/3), respectively. FF> So I wonder there is a way to get the asyptotic behaviour of a FF> function f(x), for example, Input: x^(3/2) * (1+x+x^2+x^3), and FF> request the terms at LO and NLO Output: x^(3/2) + x^(5/3) I am not sure that is NLO, but the leading term can be obtained as follows: ex e=pow(x, numeric(3,2))*(1+x+x^2+x^3); cout << (log(e).diff(x)*x).expand().subs(x==0); // -> 3/2 So now e can be multiplied by x^(2/3) to remove the singularity. See for explanations: https://en.wikipedia.org/wiki/Logarithmic_derivative So, the issue here is mathematics, not GiNaC. Once you know the right method, it shall be easier to implement in the code. -- Vladimir V. Kisil http://www.maths.leeds.ac.uk/~kisilv/ Book: Geometry of Mobius Transformations http://goo.gl/EaG2Vu Software: Geometry of cycles http://moebinv.sourceforge.net/ Jupyter (Colab): https://github.com/vvkisil/MoebInv-notebooks Jupyter (CodeOcean): https://codeocean.com/capsule/7952650/tree

1 0

Re: [GiNaC-list] series of x^n around x=0
by Feng Feng 21 Dec '19

21 Dec '19

Thanks very much for the reply.! And yes, it is not well-defined in mathematical sense to taylor expand x^(3/2) around x=0. What I exactly want is to get the asyptotic behaviour of a function f(x) around x=0. For example, f(x) = x^(3/2) (1 + x + x^2 + …), the asyptotic behaviour at leading order (LO) and next-to-leading order (NLO) are x^(3/2) and x^(5/3), respectively. So I wonder there is a way to get the asyptotic behaviour of a function f(x), for example, Input: x^(3/2) * (1+x+x^2+x^3), and request the terms at LO and NLO Output: x^(3/2) + x^(5/3) Thanks again! Best regards! Feng

2 1