- GiNaC-list - ginac.de

Problem with subs
by Tatiana Zolo 28 Feb '02

28 Feb '02

We do not understand how subs works in some cases: we tried the following ginsh session > p=x+1; 1+x > subs(p,x==k-1); k > subs(p,x==x-1); -1+x > subs(p,x==x-2); -1+x > subs(p,x==x^2-1); x^2 > subs(p,x==x^2); 1+x^2 so, apparently, subs is sometimes wrong. We couldn't find anything in the documentation to explain this fact. Thank you very much. Best wishes Alessandro Zaccagnini Tatiana Zolo purrs developers, at the University of Parma, Italy

2 1

Problems with subs
by Tatiana Zolo 28 Feb '02

28 Feb '02

We do not understand how subs works in some cases: we tried the following ginsh session > p=x+1; 1+x > subs(p,x==k-1); k > subs(p,x==x-1); -1+x > subs(p,x==x-2); -1+x > subs(p,x==x^2-1); x^2 > subs(p,x==x^2); 1+x^2 so, apparently, subs is sometimes wrong. We couldn't find anything in the documentation to explain this fact. Thank you very much. Best wishes Alessandro Zaccagnini Tatiana Zolo purrs developers, at the University of Parma, Italy

1 0

Problem with expand()
by Tatiana Zolo 28 Feb '02

28 Feb '02

We have a few questions about GinaC's behaviour. For our purposes, it would be useful to have a slightly more general form of expand that takes care of non-ambiguous iterated exponentials, as in the case (2^x)^3 which we would like to see in the form 2^(3*x) or even 8^x We understand that there are cases when this should not be done automatically, but we believe that this is an instance when there is no harm in simplifying since 2 is a positive real number and there is no question of complex numbers being involved. A "real" version of expand would be very useful for us. Do you plan to add these features to a future release of GiNaC? Thank you very much. Best wishes Alessandro Zaccagnini Tatiana Zolo purrs developers, at the University of Parma, Italy

1 0

collect TeX output
by Christophe TROESTLER 15 Feb '02

15 Feb '02

Hi, It seems that GiNaC, when (for example) collecting expressions and then printing them to LaTeX adds some braces in additions to parentheses. For example symbol x("x"), y("y"), t("t"); ex e = x*t + y*t; collect(e,t).print(GiNaC::print_latex()); ouputs t {(y+x)} Is there any way (an option?) to remove the extra braces? Indeed, they force TeX NOT to BREAK the expression (i.e., it must fit on a single line) which is problematic when the expression is long. Thanks, ChriS

4 4

GiNaC for Win32 ?
by Herve Soulard 28 Jan '02

28 Jan '02

Hi there, I'm currently looking for a free (for now) C/C++ CAS library for the Win32 platform. GiNaC's functionnalities seem OK for me, but I'm wondering if a Win32 version exists? Do you know if I can find one ? Thanks very much for any help, Hervé Soulard.

4 3

GiNaC 1.0.5 released
by Christian Bauer 27 Jan '02

27 Jan '02

Hi! You probably won't see it on national TV but it's new, improved, and finger-licking good. :-) GiNaC 1.0.5 is binary compatible to 1.0.x and has - a slightly more versatile degree()/coeff()/collect() facility - a bugfix for the output of imaginary numbers As usual, you can get it from ftp://ftpthep.physik.uni-mainz.de/pub/GiNaC Bye, Christian -- / Coding on PowerPC and proud of it \/ http://www.uni-mainz.de/~bauec002/

1 0

collect
by tatiana 25 Jan '02

25 Jan '02

Hi, I don't understand some answer about the method collect(). If it is able to collect with 'sin($1)' the polynomial 'a' why not the polynomials 'b' or 'c' with '2^$1' or '$1^n'? > a=d*p*sin(x)+p*sin(x)+q*sin(y)+q*d*sin(y)+sin(y)+d*sin(y)+q*sin(x); q*sin(y)+sin(x)*p+d*sin(y)+sin(y)+d*q*sin(y)+sin(x)*q+sin(x)*d*p > collect(a,find(a,sin($1))); sin(y)*(1+d+q+d*q)+sin(x)*(q+p+d*p) > b=q*2^y+2^x*p+d*2^y+2^y+d*q*2^y+2^x*q+2^x*d*p; q*2^y+d*2^y+2^x*d*p+2^x*p+2^y+d*q*2^y+2^x*q > collect(b,find(b,2^($1))); q*2^y+d*2^y+2^x*d*p+2^x*p+2^y+d*q*2^y+2^x*q > c=n*2^n+n-3+3^n-n^2*2^n; -3+2^n*n-2^n*n^2+3^n+n > collect(c,find(c,$1^n)); -3+2^n*n-2^n*n^2+3^n+n Thanks in advance, Tatiana Zolo

2 1

1.0.4 unleashed
by Richard B. Kreckel 24 Jan '02

24 Jan '02

Like all other 1.0.n releases, this one claims binary compatibility to 1.0.0. Cebix has opened a branch in CVS now where ABI-breaking patches may go in (but is too shy to make the announcement). Anyway, the news for 1.0.4 are: * Speedup in expand(). * Faster Bernoulli numbers (thanks to Markus Nullmeier). * Some minor bugfixes and documentation updates. .rpm's already on FTP, .deb's to follow soon thru all woody-mirrors. Enjoy -richy.

1 0

Calculation of determinant takes extremely long time.
by Ananda Murthy R S 22 Jan '02

22 Jan '02

Hello: Below I have given a program which tries to find the determinant of a matrix having real numbers and and symbols. The program compiles properly. But when I run the program it takes extremely long time to give the output. I have waited for nearly 4 hours before pressing Ctrl+C. How can I find determinants of such matrices of larger sizes? What method does GiNaC employ to calculate the determinant of a matrix having symbols? What method does GiNaC employ to calculate the characteristic polynomial of a matrix? // This is Example-1 given in the paper ``Numerical Computation of // Decentralized Fixed Modes'' by R.V.Patel and P.Misra, Automatica, // Vol.27, No.2, pp.375-382, 1991. #include <iostream> #include <ginac/ginac.h> using namespace GiNaC; using namespace std; int main() { symbol F11("F11"), F12("F12"), F13("F13"), F21("F21"), F22("F22"), F23("F23"), F31("F31"), F32("F32"), F33("F33"), F44("F44"), F45("F45"), F46("F46"), F54("F54"), F55("F55"), F56("F56"), F64("F64"), F65("F65"), F66("F66"), F77("F77"), F78("F78"), F79("F79"), F87("F87"), F88("F88"), F89("F89"), F97("F97"), F98("F98"), F99("F99"), s("s"); ex w1[8][8] = { {-2.902280128, 4.902280128, 3.698135576, -1.091710896, -1.339597088, -3.833588112, 0.1354525359, 0.1354525359}, {14.57598539, 8.663631092, 4.471485536, -6.163631092, -1.031929416, 1.692145555, -11.85191042, 0.5280631913}, {-13.01717766, 5.242845433, 6.808655713, 0.3658471056, -1.250335699, -9.401666250, 5.593010537, -0.1838540200}, {-9.358517860, 9.358517860, 9.233430521, -1.249825322, -1.102425449, -10.21076863, 0.9773381096, 0.9773381096}, {10.90157850, 10.33803799, 6.577811264, -6.527468755, -1.463336353, -1.277505940, -10.98858465, 1.391388959}, { 7.438961098, 6.026323164, 4.365195826, -4.228199859, -1.659477876, -1.136995967, -5.916479189, 0.6866298675}, { 4.590977560, 7.475549489, 4.982190956, -3.664980257, -2.282092104, -2.544374131, -5.126096156, 0.4770129007}, {16.85292654, 5.785447163, 2.204831037, -5.785447163, -1.798583916, 3.580616126, -11.47372649, 1.906247120}}; ex w2[8][9]={ { -6.288091946, -8.384122594, -10.48015324, 24.06096115, 8.869414499, 7.808913485, -2.620038311, -3.144045973, -13.10019155}, { 6.989292416, 9.319056555, 11.64882069, 16.40558839, 19.28474045, 8.940485304, 2.912205174, 3.494646208, 14.56102587}, {-13.54775782, -18.06367710, -22.57959637, 25.20750633, 18.08277192, 31.47867615, -5.644899094, -6.773878912, -28.22449547}, { -5.198332933, -6.931110578, -8.663888222, 23.84960895, 19.94822168, 17.03165000, -2.165972056, -2.599166467, -10.82986028}, { -7.648224421,-10.19763256, -12.74704070, 10.15737156, 9.939200440, 4.700713320, -3.186760175, -3.824112210, -15.93380088}, { -5.899533404, -7.866044538, -9.832555673, 19.48795695, 26.58611759, 31.63366095, -2.458138918, -2.949766702, -12.29069459}, {-6.288091946, -8.384122594, -10.48015324, 29.74830670, 25.04543540, 17.08815424, -2.620038311, -3.144045973, -13.10019155}, { 1.360132475, 1.813509966, 2.266887458, 17.73934984, 17.02037591, 14.47725433, 0.5667218645, 0.6800662374, 2.833609322}}; ex w3[9][8]={ {-2.632399367, -0.8816532905, 0.2216458770, 0.8816532905, 0.3060021972, -1.103299168, 1.223098002, 0.5276480743}, { 0.6965370444, -0.6965370444, -0.6386855981, 0.0878917469, 0.0698972579, 0.6507314098, -0.0120458116, -0.0120458116}, { 0.1251321667, 0.7413375337, 3.578384236, -0.2660527380, -0.0994345890, -3.687097343, -3.700738679, 2.940754836}, { 0.9728813679, -0.6824007787, 0.8321248154, 0.0737554812, 0.0515436172, -0.8059427381, -1.726710642, 1.440410961}, { 0.6338502954, 1.253223829, 1.874876846, -0.7779391327, -0.2226901107, -1.471703558, -0.7446180000, 1.113991924}, { 0.0725179078, -0.0725179078, -0.4003742177, 0.6048139599, -0.1581541911, 0.5700763364, -0.1697021187, -0.1697021187}, { 0.0952676686, -0.0687101141, 0.0066780601, 0.0905768308, -0.1072959550, 0.1535121891, 0.1560969942, 0.0472149362}, {-0.2252498152, 0.0992129486, 0.1715238903, 0.1417124980, -0.0477814471, -0.0948070092, 0.1258808196, 0.0636080200}, {-0.2314675213, 0.1207746598, 0.2092921292, 0.0112031403, -0.0280553044, -0.0858390289, 0.3357822914, -0.0594501770}}; ex w4[9][9]={ {F11, F12, F13, 0, 0, 0, 0, 0, 0 }, {F21, F22, F23, 0, 0, 0, 0, 0, 0 }, {F31, F32, F33, 0, 0, 0, 0, 0, 0 }, { 0, 0, 0, F44, F45, F46, 0, 0, 0 }, { 0, 0, 0, F54, F55, F56, 0, 0, 0 }, { 0, 0, 0, F64, F65, F66, 0, 0, 0 }, { 0, 0, 0, 0, 0, 0, F77, F78, F79}, { 0, 0, 0, 0, 0, 0, F87, F88, F89}, { 0, 0, 0, 0, 0, 0, F97, F98, F99}}; ex w5[8][8]={ {1,0,0,0,0,0,0,0}, {0,1,0,0,0,0,0,0}, {0,0,1,0,0,0,0,0}, {0,0,0,1,0,0,0,0}, {0,0,0,0,1,0,0,0}, {0,0,0,0,0,1,0,0}, {0,0,0,0,0,0,1,0}, {0,0,0,0,0,0,0,1}}; matrix A(8,8); matrix B(8,9); matrix C(9,8); matrix F(9,9); matrix I(8,8); //This loop fills A matrix with values. for(unsigned r=0;r<8;r++) { for(unsigned c=0;c<8;c++) { A.set(r,c,w1[r][c]); } } //This loop fills B matrix with values. for(unsigned r=0;r<8;r++) { for(unsigned c=0;c<9;c++) { B.set(r,c,w2[r][c]); } } // This loop fills C matrix with values. for(unsigned r=0;r<9;r++) { for(unsigned c=0;c<8;c++) { C.set(r,c,w3[r][c]); } } // This loop fills F matrix with values. for(unsigned r=0;r<9;r++) { for(unsigned c=0;c<9;c++) { F.set(r,c,w4[r][c]); } } // This loop fills I matrix with values. for(unsigned r=0;r<8;r++) { for(unsigned c=0;c<8;c++) { I.set(r,c,w5[r][c]); } } matrix AplusBFC = A.add(B.mul(F.mul(C))); matrix Af = I.mul_scalar(s).sub(AplusBFC); ex e1=determinant(Af); ex e2=charpoly(AplusBFC,s); cout << "R V Patel et al EXAMPLE-1:" << endl; cout << "MATRIX [A] = " << endl; cout << A << endl << endl; cout << "MATRIX [B] = " << endl; cout << B << endl << endl; cout << "MATRIX [C] = " << endl; cout << C << endl << endl; cout << "MATRIX [F] = " << endl; cout << F << endl << endl; cout << "MATRIX [AplusBFC] = " << endl; cout << AplusBFC << endl << endl; cout << "MATRIX [Af] = " << endl; cout << Af << endl << endl; cout << "Determinant |sI-AplusBFC| = " << endl; cout << e1 << endl; cout << "Characteristic polynomial of Af = " << endl; cout << e2 << endl << endl << endl; return 0; } Thanks for helping me. Sincerely, Anand

2 1

Roots of unity
by Bob McElrath 21 Jan '02

21 Jan '02

ginsh: > evalf((-1)^(1/3)); 0.5+0.86602540378443864673*I While this answer is technically correct: exp(I*pi/3) = -1 It's not what I was expecting, and the cause of some headache today. I was trying to write a little routine like RootOf to solve a cubic polynomial, and was getting three complex answers. (at least one of them must be real if you start with a real polynomial) How should I go about enforcing that (-1)^(1/3) = -1? I'm interested in polynomials over the reals, so the complex solutions don't interest me. This behaviour of evalf is preventing me from extracting the real solution. Here's my little function if anyone is interested: // Find the roots of a polynomial. // The polynomial can be at most 4th order. lst RootOf(const ex& poly, const ex& x) { ex p = poly; ex D,Q,R,S,T; p.expand(); ex a0 = poly.coeff(x,0); ex a1 = poly.coeff(x,1); ex a2 = poly.coeff(x,2); ex a3 = poly.coeff(x,3); ex a4 = poly.coeff(x,4); switch(p.degree(x)) { case 0: return lst(); case 1: return lst(-a0/a1); case 2: return lst((-a1 + sqrt(a1*a1-4*a2*a0))/(2*a2), (-a1 - sqrt(a1*a1-4*a2*a0))/(2*a2)); case 3: a0 /= a3; a1 /= a3; a2 /= a3; Q = (3*a1-a2*a2)/9; R = (9*a2*a1-27*a0-2*a2*a2*a2)/54; D = pow(Q,3) + pow(R,2); S = pow(R+sqrt(D),numeric(1,3)); T = pow(R-sqrt(D),numeric(1,3)); // These formulas come from: http://mathworld.wolfram.com/CubicEquation.html return lst( -numeric(1,3)*a2+S+T, -numeric(1,3)*a2-numeric(1,2)*(S+T)+numeric(1,2)*I*sqrt(ex(3))*(S-T), -numeric(1,3)*a2-numeric(1,2)*(S+T)-numeric(1,2)*I*sqrt(ex(3))*(S-T) ); //case 4: // http://mathworld.wolfram.com/QuarticEquation.html // It seems this is defined ambiguously. eq. 30 wants "a real root" but there // may be as many as 3! Can there 12 roots!??! default: throw(std::domain_error("RootOf: Cannot take the root of a polynomial with degree > 4")); } return lst(); } Cheers, -- Bob Bob McElrath (rsmcelrath(a)students.wisc.edu) Univ. of Wisconsin at Madison, Department of Physics

4 17