January 2002 - GiNaC-list

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

Request for a sample source file to create 8x8 matrix of real numbers.
by Ananda Murthy R S 19 Jan '02

19 Jan '02

Hello: I hereby request anybody to give me sample code to create an 8x8 matrix of real numbers using GiNaC and to print it. It would help me a lot. Thanks in advance for taking the trouble of helping me. Sincerely, Anand

2 1

Invalid conversions and error messages.
by Bob McElrath 19 Jan '02

19 Jan '02

I wrote a program that GiNaC doesn't like, and it told me so. (I'm becoming somewhat of an expert at this) ;) #include <ginac/ginac.h> #include <iostream.h> using namespace GiNaC; using namespace std; int main() { symbol a("a"); ex blah = a; numeric blahn = ex_to<numeric>(blah); // The line below generates the error: // Internal error: statement in file ./float/misc/cl_F_digits.cc, line 25 has been reached!! // Please send the authors of the program a description how you produced this error! cout << "blahn = " << blahn << endl; // If the cout is commented out, the line below generates the error: // Internal error: statement in file ./real/elem/cl_R_minusp.cc, line 30 has been reached!! // Please send the authors of the program a description how you produced this error! if(blahn.is_negative()) { cout << "blahn is negative." << endl; } } Looks to me like a cln error. But GiNaC should probably throw an exception for the invalid conversion of "a" to a numeric. Cheers, -- Bob Bob McElrath (rsmcelrath(a)students.wisc.edu) Univ. of Wisconsin at Madison, Department of Physics

2 1

Feature request: clearing lsts
by Roberto Bagnara 19 Jan '02

19 Jan '02

I did not find in the documentation a way to erase all the elements of an lst. More formally, what is needed is the functionality provided by, e.g., void clear(GiNaC::lst& l) { for (unsigned n = l.nops(); n-- > 0; ) l.remove_first(); assert(l.nops() == 0); } Something like that turns out to be very useful in several situations. For an example, consider the following: if (match(e, pattern_1, substitution)) ... else if (clear(substitution), match(e, pattern_2, substitution)) ... else if (clear(substitution), match(e, pattern_3, substitution)) ... Please, consider the addition of such a method to GiNaC::lst. You may not want to call it clear() though, as this is perhaps too STL-ish (someone may be induced to think that lst is an STL container and we do not want that ;-) All the best Roberto -- Prof. Roberto Bagnara Computer Science Group Department of Mathematics, University of Parma, Italy http://www.cs.unipr.it/~bagnara/ mailto:bagnara@cs.unipr.it

2 1