mm-314 === Subject: Re: Number of roots for polynomial ?> If you have a polynomial equation (degree between 6 and 14), is it possible> to: > - Determine if there is at least one real root, or> - Determine the number of real roots, > without actually determining all the roots ?Yes. The answers of Herman Rubin and Fernando Gonzales only addressthe case where the polynomial f has real coeffiecients. Here is acompletly different way. Let f be any analytic function, notnecessarily a polynomial, and let C be any simple closed curve in thecomplex plane that contains no root of f on the curve itself. Thenint(f'(z)/f(z), z= C) = 2*Pi*i*n where n is the number of roots of finside C. For a polynomial f, we find an upper bound M on the modulusof the roots. Then find (if possible) a lower bound m on the absolutevalue of the imaginary parts of the non-real roots. Then let C be therectangle centered at the origin with width 2M and height 2m. Performthe integration numerically. High accuracy in the integration is notrequired because the answer is guaranteed to be integer after divisionby 2*Pi*i. So just use enough accuracy so that this can round to aninteger after this last === you have a polynomial equation (degree between 6 and 14), is it possible >>to: >>- Determine if there is at least one real root, or >>- Determine the number of real roots, >>without actually determining all the roots ? > Yes. The answers of Herman Rubin and Fernando Gonzales only address > the case where the polynomial f has real coeffiecients. Here is a > completly different way. Let f be any analytic function, not > necessarily a polynomial, and let C be any simple closed curve in the > complex plane that contains no root of f on the curve itself. Then > int(f'(z)/f(z), z= C) = 2*Pi*i*n where n is the number of roots of f > inside C. For a polynomial f, we find an upper bound M on the modulus > of the roots. Then find (if possible) a lower bound m on the absolute > value of the imaginary parts of the non-real roots. Then let C be the > rectangle centered at the origin with width 2M and height 2m. Perform > the integration numerically. High accuracy in the integration is not > required because the answer is guaranteed to be integer after division > by 2*Pi*i. So just use enough accuracy so that this can round to an > integer after this last step.1) what is the provenance/name of this method? Is it an obviousextension of the aforementioned Sturm sequences (using contourintegration)? Resources/texts?2) Are there symbolic methods for manipulating such integrals? (withspecified contours like yours, rectangles or === Let f be any analytic function, not> necessarily a polynomial, and let C be any simple closed curve in the> complex plane that contains no root of f on the curve itself. Then> int(f'(z)/f(z), z= C) = 2*Pi*i*n where n is the number of roots of f> inside C. For a polynomial f, we find an upper bound M on the modulus> of the roots. Then find (if possible) a lower bound m on the absolute> value of the imaginary parts of the non-real roots. Then let C be the> rectangle centered at the origin with width 2M and height 2m. Perform> the integration numerically. High accuracy in the integration is not> required because the answer is guaranteed to be integer after division> by 2*Pi*i. So just use enough accuracy so that this can round to an> integer after this last step.> 1) what is the provenance/name of this method?The fact that Int(f'(z)/f(z), z= C) = 2*Pi*i*n is called the ArgumentPrinciple. You can find it in textbooks on Complex Analysis, forexample _Classical Complex Analysis_ by Liang-Shin Hahn and BernardEpstein (Jones and Barlett, 1996). The idea of applying numericalintegration to it and isolating the real axis is my own, though Iguess others have come with it also.> Is it an obvious extension of the aforementioned Sturm sequences (using > contour integration)?If f is a polynomial, and you completely factor it over the complexnumbers, then consider what the partial fraction decomposition off'(z)/f(z) would be. Of course we can't get the factors is general. Maybe an analogy between this partial fraction idea and Sturmsequences could be made. > 2) Are there symbolic methods for manipulating such integrals? (with> specified contours like yours, rectangles or circles)I don't think so. A symbolic method would need the partial fractiondecomposition. But note that this method gives an exact answerdespite being numerical. The answer is guaranteed to be an integer. But some symbolic computation helps with setting up the numericintegrals. Here's an example in Maple where I count the roots of arandom degree 19 polynomial inside the square with diagonal -1-I to1+I.> f:= unapply(randpoly(z, degree= 19), z); 19 15 13 12 7 4 f := z -> 58 z - 94 z - 68 z + 14 z - 35 z - 14 z> p:= unapply(D(f)(z)/f(z), z): > Digits:= 14:> evalf(> Int(p(t-I), t= -1..1, method= _CCquad)> +Int(p(1+I*t)*I, t= -1..1, method= _CCquad) > -Int(p(t+I), t= -1..1, method= _CCquad) > -Int(p(-1+t*I)*I, t= -1..1, method= _CCquad)> ); 0. + 94.247779607694 I> evalf(%/2/Pi/I); === polynomial ?>Mitch Harris harris@tcs.inf.tu-dresden.de >If you have a polynomial equation (degree between 6 and 14), is it >possible>>to:>>- Determine if there is at least one real root, or>>- Determine the number of real roots,>>without actually determining all the roots ?Yes. The answers of Herman Rubin and Fernando Gonzales only address> the case where the polynomial f has real coeffiecients. Here is a> completly different way. Let f be any analytic function, not> necessarily a polynomial, and let C be any simple closed curve in the> complex plane that contains no root of f on the curve itself. Then> int(f'(z)/f(z), z= C) = 2*Pi*i*n where n is the number of roots of f> inside C. For a polynomial f, we find an upper bound M on the modulus> of the roots. Then find (if possible) a lower bound m on the absolute> value of the imaginary parts of the non-real roots. Then let C be the> rectangle centered at the origin with width 2M and height 2m. Perform> the integration numerically. High accuracy in the integration is not> required because the answer is guaranteed to be integer after division> by 2*Pi*i. So just use enough accuracy so that this can round to an> integer after this last step.>1) what is the provenance/name of this method? Is it an obvious>extension of the aforementioned Sturm sequences (using contour>integration)? Resources/texts?>2) Are there symbolic methods for manipulating such integrals? (with>specified contours like yours, rectangles or circles)I don't know about a name, but this was close to the method used by te Riele inhis analysis of the Riemann zeta function.--irascible since === knowing any further information from the polynomial, I think you canonly say:Complex roots come in pairs, so if your polynomial has odd degree, then ithas to have at least one real root. If the polynomial degree is even, thenit may not have real roots.Fernando.ckruger (degree between 6 and 14), is itpossible> to:> - Determine if there is at least one real root, or> - Determine the number of real roots,> without actually determining all the roots ?> === polynomial ?>If you have a polynomial equation (degree between 6 and 14), is it possible>to:>- Determine if there is at least one real root, or>- Determine the number of real roots,>without actually determining all the roots ?This is old stuff, and used to be routinely taught incourses called Theory of Equations. One can use thisfor more than polynomials.Suppose one has continuous functions f_0, ..., f_k on aninterval (a, b) on the extended real line, with the propertythat f_0 has only simple roots (for polynomials, the resultmodifies if there are multiple roots), at roots of f_0, thederivative exists and f_1 has the same sign as thederivative of f_0, and for all j > 1, f_j(x) = f_{j-1}(x)*g_j(x) - f_{j-1}(x).Furthermore, let f_k have constant sign on (a,b). Fora <= x <=b, let s(x) be the number of sign changes inthe sequence of f's, possibly taking limits at endpoints.The number of roots in the interval is then s(a) - s(b).For polynomials, one can use f_1(x) to be the derivativeof f_0, and keep the degrees decreasing until a constantis reached. This is the original Sturm sequence.-- This address is for information only. I do not claim that these viewsare those of the Statistics Department or of Purdue University.Herman Rubin, Department of Statistics, Purdue Universityhrubin@stat.purdue.edu Phone: (765)494-6054 FAX: === Problem is that I haveto do it for the general case (a7*x^7+a6*x^6....), and the whole thingquickly becomes unmanagable. I gave up when I got to the third long divisionstep.Let me explain more, maybe there is another way to do it:The parameters ai of the equation must be selec so that the equation isalways negative and touches 0 at at least one point, which is known. It maytouch at more points (unknown), but it must never cross 0 and becomespossitive (0 to +oo). All but one of the parameters can be calcula fromother criteria (known points or slopes), but currently I manually adjust thelast parameter so that it never crosses 0 (I view it graphically whilemaking the change). I wish to eliminate this manual step, and somehow comeup with an equation that relates this last parameter to the others, based onthe criteria that the equation is always negative or (degree between 6 and 14), is itpossible>to:>- Determine if there is at least one real root, or>- Determine the number of real roots,>without actually determining all the roots ?> This is old stuff, and used to be routinely taught in> courses called Theory of Equations. One can use this> for more than polynomials.> Suppose one has continuous functions f_0, ..., f_k on an> interval (a, b) on the extended real line, with the property> that f_0 has only simple roots (for polynomials, the result> modifies if there are multiple roots), at roots of f_0, the> derivative exists and f_1 has the same sign as the> derivative of f_0, and for all j > 1,> f_j(x) = f_{j-1}(x)*g_j(x) - f_{j-1}(x).> Furthermore, let f_k have constant sign on (a,b). For> a <= x <=b, let s(x) be the number of sign changes in> the sequence of f's, possibly taking limits at endpoints.> The number of roots in the interval is then s(a) - s(b).> For polynomials, one can use f_1(x) to be the derivative> of f_0, and keep the degrees decreasing until a constant> is reached. This is the original Sturm sequence.> --> This address is for information only. I do not claim that these views> are those of the Statistics Department or of Purdue University.> Herman Rubin, Department of Statistics, Purdue University> hrubin@stat.purdue.edu === of roots for polynomial ?> yes, I tried calculating the Sturm sequences. Problem is that I have> to do it for the general case (a7*x^7+a6*x^6....), and the whole thing> quickly becomes unmanagable. I gave up when I got to the third long division> step.> Let me explain more, maybe there is another way to do it:> The parameters ai of the equation must be selec so that the equation is> always negative and touches 0 at at least one point, which is known. It may> touch at more points (unknown), but it must never cross 0 and becomes> possitive (0 to +oo).It appears then to be an easier problem than as originally sta. Tell me if I have this correct: You have a polynomial with realcoefficients. You want the polynomial to be nonpositive just on(0,infinity) (open interval). You don't care whether it becomespositive for negative x. The polynomial has a known root at some a >0, and may have some other roots. Then your polynomial can beexpressed as f(x) = (x-a)^2*g(x), and it remains to check whether g(x)is nonpositive on (0,infinity). Certainly the leading coefficient ofg(x) (same as f(x)) must be negative. If all coefficients of g(x) arenegative, then you're done. You have the freedom to choose exactlyone of the coefficients *independently* of the compu values of theothers, correct? Which coefficient does that tend to be?I'll try to fill in some more details later if someone else doesn'tget to === yes, I tried calculating the Sturm sequences. Problem is that Ihave> to do it for the general case (a7*x^7+a6*x^6....), and the whole thing> quickly becomes unmanagable. I gave up when I got to the third longdivision> step.Let me explain more, maybe there is another way to do it:The parameters ai of the equation must be selec so that the equationis> always negative and touches 0 at at least one point, which is known. Itmay> touch at more points (unknown), but it must never cross 0 and becomes> possitive (0 to +oo).> It appears then to be an easier problem than as originally sta.> Tell me if I have this correct: You have a polynomial with real> coefficients. You want the polynomial to be nonpositive just on> (0,infinity) (open interval). You don't care whether it becomes> positive for negative x. The polynomial has a known root at some a 0, and may have some other roots. Then your polynomial can be> expressed as f(x) = (x-a)^2*g(x), and it remains to check whether g(x)> is nonpositive on (0,infinity).Yes to all of above.> Certainly the leading coefficient of> g(x) (same as f(x)) must be negative. If all coefficients of g(x) are> negative, then you're done.No, all are not negative.> You have the freedom to choose exactly> one of the coefficients *independently* of the compu values of the> others, correct? Which coefficient does that tend to be?there are many coefficients (usually 8), I do not have 8 variables. I onlyhave 4 variables, which, in various combinations, make up the 8coefficients. The variable I vary to get the curve negative at all points,occur in about half the === :-([Q,R,E] = qr(A,0) for full matrix A, produces an economy-sizedecomposition in which E is a permutation vector, so that A(:,E) =Q*R. The column permutation E is chosen so that abs(diag(R)) isdecreasing.This command is perfect ... but i need that abs(diag(R)) isINCREASING!!!There is a specific command? If there in not a specific command, where can i find === Re: [MatLab] QR factorization ... :-(Hi Luca,[Q,R,E] = qr(A,0) for full matrix A, produces an economy-size> decomposition in which E is a permutation vector, so that A(:,E) => Q*R. The column permutation E is chosen so that abs(diag(R)) is> decreasing.> This command is perfect ... but i need that abs(diag(R)) is> INCREASING!!!This is not always possible. For example, there is no permuQR decomposition of [1 1; 0 0] so === Iterative solution of Ax=b, what's wrong?> I'm sure this was === of no help to you.Actually it _is_ of help. SorinSubject: Re: Does anyone use functional programming languagesX-Enigmail-Version: 0.76.1.0X-Enigmail-Supports: languages like Objective Caml or SISAL> are still completely unknown to scientists.> They are extremely good (if not much better) alternatives to Fortran and> C++, in terms of expressive power, easiness, and safety, and their> performance is known to be comparable with C.> Yet Fortran is still taught to physics students like it was the alpha and> the omega of numerical programming. Why is it so ?Think about natural languages: they already fill niches very well, soa new improved language like Esperanto has a slim chance ofreaching a critical mass.To get any chance of succeeding new programming languages shouldfrom the beginning provide a huge advantage compensating the loss of decadesof expertises contained in the already available libraries, in thetrained people, as well as in the compiler technology. Now to makethe situation worse, the many functional languages compete with each others.Changing of programming language is like a revolution for many largescientific computing programs because the typical development timeand usage of such programs amounts often to decades.So continuity and forward compatibility are important features.By adopting a slowly adapting evolution, Fortran provides the continuitythat is required, explaining its very long life wrt to other dead languages(Algol, Basic, APL, PL1, Pascal, Modula, ADA, etc.).Despite what you seem to believe, Fortran 95, soon 2000+, remains avery good and easy to learn language for scientific computations. =?ISO-8859-1?Q?om=22?= are still completely unknown to scientists.> They are extremely good (if not much better) alternatives to Fortran and> C++, in terms of expressive power, easiness, and safety, and their> performance is known to be comparable with C.> Yet Fortran is still taught to physics students like it was the alpha and> the omega of numerical programming. Why is it so ?Think about natural languages: they already fill niches very well, soa new improved language like Esperanto has a slim chance ofreaching a critical mass.To get any chance of succeeding new programming languages shouldfrom the beginning provide a huge advantage compensating the loss of decadesof expertises contained in the already available libraries, in thetrained people, as well as in the compiler technology. Now to makethe situation worse, the many functional languages compete with each others.Changing of programming language is like a revolution for many largescientific computing programs because the typical development timeand usage of such programs amounts often to decades.So continuity and forward compatibility are important features.By adopting a slowly adapting evolution, Fortran provides the continuitythat is required, explaining its very long life wrt to other dead languages(Algol, Basic, APL, PL1, Pascal, Modula, ADA, etc.).Despite what you seem to believe, Fortran 95, soon 2000+, remains avery good and easy === matlab programming questionWhat I want is the following:[y] = function(N) % N integerH = H(a_j1,a_j2, ..., a_jN)> % jk runs from 1 to Ny(a_j1,a_j2, ..., a_jN) = sqrt(a_j1^2 + a_j2^2 + ... + a_jN^2)Wilbert, I can't say that I completely understand what you're trying to do,> but it seems like you are asking how to index into H with N indices, when N> is not known until run time. You can use a cell array for that. Here's an> example where I create a 3x3x3 H array, but then I index into it without> explicitly writing an expression having 3 indices:> for i=1:3, for j=1:3, for k=1:3, h(i,j,k) = 100*i + 10*j + k; end; end;> end>> ind = {2,1,3};>> h(ind{:})> ans => 213If that's not the gist of your question, please let me know.> The problem is that 3 is not arbitrary. What I want is to make the> following function (or one which produces the same output):> [h] = function(N)> for j1=1:3, ..., for jN=1:3,> h(j1, ..., jN) = (j1 + ... + jN)*sqrt(j1^2 + ... + jN^2);> end; ... end;You can use a loop while 1, .... if counter(1)>3, break; end; end;in which you create a counter containing the index list, which is updalexicographically. This simulates N nes loops. I have been using this attimes, and it works well (though nes loops are slow in Matlab and shouldbe avoided if === software> Does somebody recommend journal artical or books on the design> of test cases for numerical software?http://www.mat.univie.ac.at/~neum/glopt/ presentation.html[~ is a tilde] contains some references, mostly rela to numerical software in optimization.Arnold === NeumaierSubject: Statistics on sex = HOW IS THIS men before death. Is thisstatistically possible? My girlfriend claims it is so, but her examplesfall apart with a bit of scrutiny. The average for both groups changes bythe same factor every time two people have sex. She then says it is due toan uneven population, but I argue that the population would have to be awhole hell of a lot more disproportionate for a ratio like 4:1. I think itis more likely that both groups lied on the survey, women repor lowernumbers and men repor higher === POSSIBLE????quas litteras Master of all Photosohp in continuosci.math.num-analysis:66516 scripsit eis responsurus sum> before death, the average woman has sex with 3 men before death. Is this> statistically possible? If the male heterosexual population is four times smaller than the female one,it is. I do not deem it is so, though.A few hints:i) men's arrogance and women's decence have twis the figures [In that case, half of the woman have an average of more than 21 partners. B....! ;-)]> My girlfriend claims it is soDoes she === Statistics on sex = HOW IS THIS POSSIBLE????>before death, the average woman has sex with 3 men before death. Is this>statistically possible? My girlfriend claims it is so, but her examples>fall apart with a bit of scrutiny. The average for both groups changes by>the same factor every time two people have sex. She then says it is due to>an uneven population, but I argue that the population would have to be a>whole hell of a lot more disproportionate for a ratio like 4:1. I think it>is more likely that both groups lied on the survey, women repor lower>numbers and men repor higher numbers.I guess that the average man does not have the average number of sexpartners. The average man would be the median man: half of the men had sex with12 women of more, half had sex with 11 or less.That's often what is meant by average man. Otherwise, it is notpossible, if you just count the number of first times 2 personshad sex together: F, then if there are M men and W women, F = 12M,and F = 3W, so M = 4W, and we know === Fourier MethodHi! I'm trying to solve the beam propagation equation in a nonlineardispersive media without using the slowly varying approximation. Forthis, i'm using the method illustraded by Dr. G.P.Agarwal in his bookNon-linear Fiber Optics, Chap 2, Eqn 2.4.8 with my nonlinearoperator given as per eqn 2.4.3 i.e.N= i*(nonlinearity coefft)*[ |A|^2*A +2i/w0 d/dt(|A|^2*A)-TrAd/dT(|A|^2)]If i neglect the last 2 terms, i get the standard NL SchrodingerEquation for which i have written the Matlab code. How do iincorporate the next two terms and more specifically, how do ievaluate exp(h*N(z)) with step size h? Do i take the operator toFourier space, as Dr. Agarwal has done for the Dispersive operator? Any references === whether the point is within triangle(given with 3 points)> What is the best algorithm of detecting whether the point is within> triangle(given with 3 points)?> I use the vector product: there is a triangle (x1, y1, x2, y2, x3, y3)> and a point (x, y); then> a = (x2 - x1) * (y - y1) - (y2 - y1) * (x - x1);> b = (x3 - x2) * (y - y2) - (y3 - y2) * (x - x2);> c = (x1 - x3) * (y - y3) - (y1 - y3) * (x - x3);> if a, b and c have the same sign (+ or -), then the point is> within the triangle, if they have different signs, then it is untrue.> But I think there is the better (more precise) algorithm; what would> you say?You can do this with complex arithmetic rather easily. Makes a nicesimple algorithm in Fortran or any other language that supportscomplex arithmetic.-- ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ Science knows only one commandment: contribute to science. -- === point is within triangle(given with 3 points)> What is the best algorithm of detecting whether the point is within> triangle(given with 3 points)?> I use the vector product: there is a triangle (x1, y1, x2, y2, x3, y3)> and a point (x, y); then> a = (x2 - x1) * (y - y1) - (y2 - y1) * (x - x1);> b = (x3 - x2) * (y - y2) - (y3 - y2) * (x - x2);> c = (x1 - x3) * (y - y3) - (y1 - y3) * (x - x3);> if a, b and c have the same sign (+ or -), then the point is> within the triangle, if they have different signs, then it is untrue.> But I think there is the better (more precise) algorithm; what would> you say?At my web site - see signature below, there is the Fortran program locpt.f90for determining if a point is inside a polygon.It is under the code from the Naval Surface Warfare Millerhttp://users.bigpond.net.au/amillerRetired === within triangle(given with 3 points)> What is the best algorithm of detecting whether the point is within> triangle(given with 3 points)?> I use the vector product: there is a triangle (x1, y1, x2, y2, x3, y3)> and a point (x, y); then> a = (x2 - x1) * (y - y1) - (y2 - y1) * (x - x1);> b = (x3 - x2) * (y - y2) - (y3 - y2) * (x - x2);> c = (x1 - x3) * (y - y3) - (y1 - y3) * (x - x3);> if a, b and c have the same sign (+ or -), then the point is> within the triangle, if they have different signs, then it is untrue.> But I think there is the better (more precise) algorithm; what would> you say?Let the triangle having vertices A=(x_1,y_1) , B=(x_2,y_2) , C=(x_3,y_3)and in 2D consider a point M=(x,y)Consider the determinants | x-x_3 x_2-x_3|D_1=| | | y-y_3 y_2-y_3| | x_1-x_3 x-x_3|D_2=| | | y_1-y_3 y-y_3| | x_1-x_3 x_2-x_3|D= | | . | y_1-y_3 y_2-y_3| Because ABC is a trinangle we have D=/=0 . D=0 means that AC is parallelto BC. Define p=D_1/D and q= D_2/D .Test if (*) 0 =< p =< 1 and 0=< q =< 1 .If both inequalities are verified, then M is within triangle ABC (or on sides). If all inequalities are strict,then M Sis in interior of ABC. A motivation is that when (*) hold, then there are three positive numbers p,q,r such that p+q+r=1 and x=p*x_1+q*x_2+r*x_3 y=p*y_1+q*y_2+r*y_3 that is M belongs to === occured to anyone except myself that statements like4-3=7-6 can't be true ?(The difference between 3 and 4 can't be the same as the difference between6 and 7).Or am I just having a bad day ?1) If we take the set of all possible numbers, by the definition of a numbereach one must be unique from every other one. Because if it wasn't uniquelydifferent it'd have to be the same number (Have I got this right ?)2) Because of 1) The difference between every number and every other numbermust be unique.So doesn't that mean that 4-3 can't be equal to 7-6, because if it was thedifference between the two sides of the equality is the same, hence therelationship between the two numbers on both sides of the equation must bethe same, hence both pairs of numbers on both sides of the equation must bethe same numbers. Since 7 and 4 aren't the same number, how can the === Has it occured to anyone except myself that statements like> 4-3=7-6 can't be true ?> (The difference between 3 and 4 can't be the same as the difference between> 6 and 7)....> 1) If we take the set of all possible numbers, by the definition of a number> each one must be unique from every other one. Because if it wasn't uniquely> different it'd have to be the same number (Have I got this right ?)> 2) Because of 1) The difference between every number and every other number> must be unique.> So doesn't that mean that 4-3 can't be equal to 7-6, because if it was the> difference between the two sides of the equality is the same, hence the> relationship between the two numbers on both sides of the equation must be> the same, hence both pairs of numbers on both sides of the equation must be> the same numbers. Since 7 and 4 aren't the same number, how can the equality> hold ?(This should really be in sci.math, as it doesn't haveanything to do with numerical analysis. Never mind.)Difference doesn't mean quite the same as relationship.Nor does difference, *in this context*, mean the same asit does in everyday English.When we say the difference between 4 and 3 is 1 we don'tmean that the number 1 somehow encapsulates *everything*that's not the same about 3 and 4: only that it describes*how much bigger* 4 is than 3. You're free to dislike thisuse of the word difference, but that's how it is.So: - Indeed, the way in which 4 is different from 3 isn't the same as the way in which 7 is different from 6. (For instance: 4 is a square and 3 isn't is true but 7 is a square but 6 isn't is false.) - But it doesn't follow that the difference between 4 and 3 and the difference between 7 and 6 can't be the same, because the difference between is a technical term here and means much less than it does in ordinary English.There are lots of other examples of technical termsin mathematics that have meanings quite different fromtheir meanings in ordinary language. A group is anentirely different thing from a set; a ring, acycle and a circle have nothing to do with eachother; something can be compact even if it islarger than the sun; a number can be real eventhough there is nothing in the physical world towhich it exactly matches up.To understand mathematics, you need to be able tolet go of the everyday meanings of some words, andget used to their technical meanings.-- Gareth McCaughan.sig under === POSSIBLE????X-NewsReader: GRn 3.2n February 9, 1999> before death, the average woman has sex with 3 men before death. Is this> statistically possible? My girlfriend claims it is so, but her examples> fall apart with a bit of scrutiny. The average for both groups changes by> the same factor every time two people have sex. She then says it is due to> an uneven population, but I argue that the population would have to be a> whole hell of a lot more disproportionate for a ratio like 4:1. I think it> is more likely that both groups lied on the survey, women repor lower> numbers and men repor higher === algorithm.scheduling jobs algorithm.I need to weekly schedule jobs.I have 9 people that need scheduled over a period of 6 days.The number of people that need to be working depends from day to day,and even during one day.SO i have decided to split a day in 15min periodes.ex. monday 7:30am need 3 people7:45am need 3 people8:00am need 5 people8:00am need 5 people1:00pm need 6 people...6:00pm need 3 peopleThe number of hours that people work is also different.6 people work 24 hours/working week2 people work 32 hours/working week1 person works 38 hours/working weekWhat kind of algorithm should i look for? Linaer programming? ...?I am not math wonder so any basic help === scheduling jobs algorithm.> I need to weekly schedule jobs.> I have 9 people that need scheduled over a period of 6 days.> The number of people that need to be working depends from day to day,> and even during one day.> SO i have decided to split a day in 15min periodes.> ex. monday> 7:30am need 3 people> 7:45am need 3 people> 8:00am need 5 people> 8:00am need 5 people> 1:00pm need 6 people> ...> 6:00pm need 3 people> The number of hours that people work is also different.> 6 people work 24 hours/working week> 2 people work 32 hours/working week> 1 person works 38 hours/working week> What kind of algorithm should i look for? Linaer programming? ...?> I am not math wonder so any basic help is a start.Most employees are not going to be happy with a schedulethat asks them to come in for 15 minutes or half an hour at atime, then leave and come back later.You probably want to figure out shifts for each day, suchthat the overlaid shifts give the desired depth of staff. Ie.have three people come in at 7:30am, another two shifts startat 8, etc.Don't be surprised if you have to compromise your depthchart a bit in order to make the schedule work for youremployees. You may, for example, have a need for 6 staffrather than 5 starting at 1pm. Depending on the resourcesalloca to your schedule, you might overstaff for someperiods.Anyway, the suggestion to create the daily shifts (usingpaper cutouts or side-by-side columns on graph paper),then combine the daily shifts into weekly allocations === Convergence test in the Bisection algorithm> I downloaded from Netlib a routine that computes the eps (unit round off)> for the C double type and it gives me 10^-19.> But whne I try to to arrange a convergence check in the bisection algorithm> that stop when abs(b-a) I have taken a look at fortran implementation and they usually use> abs(b-a)<100*eps. What is the purpose of that 100?> I have also arranged another test that seem to work better: I compare> abs(b-a) with abs(b-a) of the previous iteration and if they are equal it> has converged. The strange thing is that the final abs(b-a) is different for> every couple of starting b,a ! It usually attains an abs(b-a) in [3e-16 ,> 10e-16].> Note that I want the maximum precision since I want idea to use eps in a convergence test.The only sure way to achieve the maximum accuracy is to usea stopping criterion based on monotonicity.In this case, you can stop as soon as - the lower bound stops to increase strictly - the upper boundstops to decrease strictly - the interval's width stops to decrease strictlyBernard DanloyDepartment of Mathematical === Convergence test in the Bisection algorithm> I downloaded from Netlib a routine that computes the eps (unit round off)> for the C double type and it gives me 10^-19.I bet those doubles are being kept in 80-bit floating-point registers,with 64-bit mantissas ...> But whne I try to to arrange a convergence check in the bisection algorithm> that stop when abs(b-a) I have taken a look at fortran implementation and they usually use> abs(b-a)<100*eps. What is the purpose of that 100?To avoid the problem you've just encountered, which can occur evenwhen the epsilon is only slightly too small (as opposed to === compiling Scalapack codeGood morning! I've used Scalapack for a while and have never reallyhad too many problems (OK, I've pos occasionally :).Now, I'm trying to compile a code that calls Scalapack.Here are the errors messages during the link phase:[laytojb@grendel PROJ2A]$ make -f makefile.hard.grendel/home1/laytojb/bin/mpif90 -ifc -ifc -O3 -tpp7 -prefetch -unroll -align -axW -mp -fp_port -pc64 -xW -Vaxlib -o xdscaex pz_laread_mod.o pz_lawrite_mod.o pz_solver.o pz_prob_info.o pz_laread.o pz_lawrite.o myseconds.o xerbla.f -L/home1/laytojb/lib -lscalapack -lpblas_LINUX -lblacsF77init_MPI-LINUX-0 -lblacs_MPI-LINUX-0 -L/home1/laytojb/lib -lblacsF77init_MPI-LINUX-0 -lblacsCinit_MPI-LINUX-0 -llapack -lgoto_p4_512-r0.6 external subroutine XERBLA ^Comment 15 at (9:xerbla.f) : This feature is obsolescent in Fortran 9543 Lines Compiled/home1/laytojb/lib/libscalapack.a(pzgetrf.o): In function `pzgetrf_':pzgetrf.o(.text+0x17b): undefined reference to `pb_topget_'pzgetrf.o(.text+0x19a): undefined reference to `pb_topget_'pzgetrf.o(.text+0x1bc): undefined reference to `pb_topget_'pzgetrf.o(.text+0x1de): undefined reference to `pb_topset_'pzgetrf.o(.text+0x200): undefined reference to `pb_topset_'pzgetrf.o(.text+0x21f): undefined reference to `pb_topset_'pzgetrf.o(.text+0x9c3): undefined reference to `pb_topset_'pzgetrf.o(.text+0x9e2): undefined reference to `pb_topset_'/home1/laytojb/lib/libscalapack.a(pzgetrf.o)(.text +0xa04): more undefined references to `pb_topset_' follow/home1/laytojb/lib/libscalapack.a(pzgetf2.o): In function `pzgetf2_':pzgetf2.o(.text+0x177): undefined reference to `pb_topget_'make: *** [example1] Error 1The machine is a dual 2.4 GHz/P4 running Redhat Linux 7.1.The compilers are Intel Compilers 7.1-16. I built LAPACK,BLACS, and SCALAPACK as I have before (I can provide includefiles and makefiles if needed). However, for some reason, onthis machine, it always gives me link errors. I have builtSCALAPACK and successfully compiled the code on a home builtcluster and a large cluster of PIII/1.2 CPUs (also runningRedhat 7.3) with the same compilers. For same reason it justdoesn't like this === compiling Scalapack code> Good morning!> I've used Scalapack for a while and have never really> had too many problems (OK, I've pos occasionally :).> Now, I'm trying to compile a code that calls Scalapack.> Here are the errors messages during the link phase:> [laytojb@grendel PROJ2A]$ make -f makefile.hard.grendel> /home1/laytojb/bin/mpif90 -ifc -ifc -O3 -tpp7 -prefetch -unroll -align > -axW -mp -fp_port -pc64 -xW -Vaxlib -o xdscaex pz_laread_mod.o > pz_lawrite_mod.o pz_solver.o pz_prob_info.o pz_laread.o pz_lawrite.o > myseconds.o xerbla.f > -L/home1/laytojb/lib -lscalapack -lpblas_LINUX > -lblacsF77init_MPI-LINUX-0 -lblacs_MPI-LINUX-0 -L/home1/laytojb/lib > -lblacsF77init_MPI-LINUX-0 -lblacsCinit_MPI-LINUX-0 -llapack > -lgoto_p4_512-r0.6> external subroutine XERBLA> ^> Comment 15 at (9:xerbla.f) : This feature is obsolescent in Fortran 95> 43 Lines Compiled> /home1/laytojb/lib/libscalapack.a(pzgetrf.o): In function `pzgetrf_':> pzgetrf.o(.text+0x17b): undefined reference to `pb_topget_'> pzgetrf.o(.text+0x19a): undefined reference to `pb_topget_'> pzgetrf.o(.text+0x1bc): undefined reference to `pb_topget_'> pzgetrf.o(.text+0x1de): undefined reference to `pb_topset_'> pzgetrf.o(.text+0x200): undefined reference to `pb_topset_'> pzgetrf.o(.text+0x21f): undefined reference to `pb_topset_'> pzgetrf.o(.text+0x9c3): undefined reference to `pb_topset_'> pzgetrf.o(.text+0x9e2): undefined reference to `pb_topset_'> /home1/laytojb/lib/libscalapack.a(pzgetrf.o)(.text+0xa04): more > undefined references to `pb_topset_' follow> /home1/laytojb/lib/libscalapack.a(pzgetf2.o): In function `pzgetf2_':> pzgetf2.o(.text+0x177): undefined reference to `pb_topget_'> make: *** [example1] Error 1> The machine is a dual 2.4 GHz/P4 running Redhat Linux 7.1.> The compilers are Intel Compilers 7.1-16. I built LAPACK,> BLACS, and SCALAPACK as I have before (I can provide include> files and makefiles if needed). However, for some reason, on> this machine, it always gives me link errors. I have built> SCALAPACK and successfully compiled the code on a home built> cluster and a large cluster of PIII/1.2 CPUs (also running> Redhat 7.3) with the same compilers. For same reason it just> doesn't like this machine.> Any ideas? Any one?> TIA!> JeffLook at the file SLmake.inc, it should be an option -DAdd_ somewhere.I think === compiling Scalapack codegilles,> Look at the file SLmake.inc, it should be an option -DAdd_ somewhere.> I think you should === gilles,>> Look at the file SLmake.inc, it should be an option -DAdd_ somewhere.>> I think you should tak a look there.> It's define it in the CDEFS variableHere is an extract of my SLmake.inc :...# C preprocessor defs for compilation# (-DNoChange, -DAdd_, -DUpCase, or -Df77IsF2C)#CDEFS = -DAdd_ === compiling Scalapack code> Good morning!> I've used Scalapack for a while and have never really> had too many problems (OK, I've pos occasionally :).> Now, I'm trying to compile a code that calls Scalapack.> Here are the errors messages during the link phase:> [laytojb@grendel PROJ2A]$ make -f makefile.hard.grendel> /home1/laytojb/bin/mpif90 -ifc -ifc -O3 -tpp7 -prefetch -unroll -align -axW -mp -fp_port -pc64 -xW -Vaxlib -o xdscaex > pz_laread_mod.o pz_lawrite_mod.o pz_solver.o pz_prob_info.o pz_laread.o pz_lawrite.o myseconds.o xerbla.f > -L/home1/laytojb/lib -lscalapack -lpblas_LINUX -lblacsF77init_MPI-LINUX-0 -lblacs_MPI-LINUX-0 -L/home1/laytojb/lib > -lblacsF77init_MPI-LINUX-0 -lblacsCinit_MPI-LINUX-0 -llapack -lgoto_p4_512-r0.6> external subroutine XERBLA> ^> Comment 15 at (9:xerbla.f) : This feature is obsolescent in Fortran 95> 43 Lines Compiled> /home1/laytojb/lib/libscalapack.a(pzgetrf.o): In function `pzgetrf_':> pzgetrf.o(.text+0x17b): undefined reference to `pb_topget_'> pzgetrf.o(.text+0x19a): undefined reference to `pb_topget_'> pzgetrf.o(.text+0x1bc): undefined reference to `pb_topget_'> pzgetrf.o(.text+0x1de): undefined reference to `pb_topset_'> pzgetrf.o(.text+0x200): undefined reference to `pb_topset_'> pzgetrf.o(.text+0x21f): undefined reference to `pb_topset_'> pzgetrf.o(.text+0x9c3): undefined reference to `pb_topset_'> pzgetrf.o(.text+0x9e2): undefined reference to `pb_topset_'> /home1/laytojb/lib/libscalapack.a(pzgetrf.o)(.text+0xa04): more undefined references to `pb_topset_' follow> /home1/laytojb/lib/libscalapack.a(pzgetf2.o): In function `pzgetf2_':> pzgetf2.o(.text+0x177): undefined reference to `pb_topget_'> make: *** [example1] Error 1> The machine is a dual 2.4 GHz/P4 running Redhat Linux 7.1.> The compilers are Intel Compilers 7.1-16. I built LAPACK,> BLACS, and SCALAPACK as I have before (I can provide include> files and makefiles if needed). However, for some reason, on> this machine, it always gives me link errors. I have built> SCALAPACK and successfully compiled the code on a home built> cluster and a large cluster of PIII/1.2 CPUs (also running> Redhat 7.3) with the same compilers. For same reason it just> doesn't like this machine.> Any ideas? Any one?well, firstly, look where is 'pb_topset_' supposed to be definedin your succesful Redhat 7.3 build? it may also be one of those underscore problems, one library may becompiled with leading underscores and other callers not.I think Intel Fortran has options to control this. === morning! I've used Scalapack for a while and have never really> had too many problems (OK, I've pos occasionally :).> Now, I'm trying to compile a code that calls Scalapack.> Here are the errors messages during the link phase:[laytojb@grendel PROJ2A]$ make -f makefile.hard.grendel> /home1/laytojb/bin/mpif90 -ifc -ifc -O3 -tpp7 -prefetch -unroll -align -axW -mp -fp_port -pc64 -xW -Vaxlib -o xdscaex > pz_laread_mod.o pz_lawrite_mod.o pz_solver.o pz_prob_info.o pz_laread.o pz_lawrite.o myseconds.o xerbla.f > -L/home1/laytojb/lib -lscalapack -lpblas_LINUX -lblacsF77init_MPI-LINUX-0 -lblacs_MPI-LINUX-0 -L/home1/laytojb/lib > -lblacsF77init_MPI-LINUX-0 -lblacsCinit_MPI-LINUX-0 -llapack -lgoto_p4_512-r0.6> external subroutine XERBLA ^> Comment 15 at (9:xerbla.f) : This feature is obsolescent in Fortran 9543 Lines Compiled pb_topget is part of the pblas_LINUX library. Chances are the variableCDEFSin your SLmake.inc contains -Df77IsF2C , where you actually need-DAdd_ to work correctly with === /home1/laytojb/bin/mpif90Why don't you use an f77?V.-- mail me === Scalapack codeDistribution: inet> /home1/laytojb/bin/mpif90> Why don't you use an f77?Suppose I'll throw in the stock reply to that kind of question.Every f90 compiler *IS* an f77 compiler.Now it wouldn't shock me if Victor's comment might suggest somethingdeeper...namely the possibility that you might have ended up mixingcode compiled by different compilers. That can cause problems.Though perhaps I'm reading more into an f77 than is there.Afraid I don't actually have much clue as the the OP's problem.I just have a knee-jerk reaction when anyone suggests that anf90 compiler doesn't count as an f77 one.-- Richard Maine | Good judgment comes from experience;email: my first.last at org.domain | experience comes from bad judgment.org: nasa, domain: gov | -- === codeDistribution: inet/home1/laytojb/bin/mpif90Why don't you use an f77?> Suppose I'll throw in the stock reply to that kind of question.> Every f90 compiler *IS* an f77 compiler.In that they accept the same syntax. However, they may use differentsystem library calls to implement those.I can not suggest a specific scenario that would let an f77 compilersolve problem here, but I wouldn't be at all surprised === Problems compiling Scalapack codeDistribution: inet> /home1/laytojb/bin/mpif90Why don't you use an f77?Suppose I'll throw in the stock reply to that kind of question.> Every f90 compiler *IS* an f77 compiler.> In that they accept the same syntax. However, they may use different> system library calls to implement those.And so might two different f77 compilers...or two different f90compilers. The whole question of f77 vs f90 just seems completelyirrelevant to the question at hand.That a library might be compiled with a particular compilercould plausibly be an issue. Had you asked Why don't youuse the xxx f77 compiler?, where xxx was some specific choice,I wouldn't have thought it at all strange. I'd have evenconsidered it a possibly good match for the symptoms.But to me, the an in Why don't you use an f77? implies that anyf77 compiler would do, apparently just so long as it wasn't also anf90 compiler. I suppose I'm reading a lot into an vs the, but thedistinction can completely change the meaning of a sentence. Perhapsthe meaning that I saw there wasn't intended.-- Richard Maine | Good judgment comes from experience;email: my first.last at org.domain | experience comes from bad judgment.org: nasa, domain: gov | -- Mark === for x>1e10 accuratelyI have read the paper by K.C. Ng but it looks very complica for a poor programmer's mind.Is there any implementation (in C or in any other language) of the === calculate sin/cos x for x>1e10 accuratelyWhy? It is unusual === sin(1e22)X-Enigmail-Version: 0.76.1.0X-Enigmail-Supports: pgp-inline, pgp-mime> I would like to calculate sin/cos x for x>1e10 accurately> I have read the paper by K.C. Ng but it looks very complica for a poor programmer's mind.> Is there any implementation (in C or in any other language) of the for your help...I don't know what range-reduction algorithm it uses, but for this sort of thingNetlib is your friend. http://cm.bell-labs.com/netlib/fdlibm/index.htmlThere is some gorgeously yucky range-reduction code in the *rem_pio2.c === Chebyshev approximation problemI have a function f(x,y,z) which I need to approximate through numericalmethods, starting from a discrete lattice f[i,j,k].f(x, constant y, constant z) would best me modelled by a rational Chebyshev approximationf(constant x, y, constant z) can be done by simple Chebyshev polynomialsf(constant x, constant y, z) is even angular, and so is optimally described througha series of Cosine Fourier coefficients (found by FFTs or something)Therefore, the best approximation choice is Sum_{i,j,k} a_{i,j,k} T_i (x) T_j (y) Cos(kz)f(x,y,z) ~ ------------------------------------------------- Sum_{l} b_{l} T_l (x)(Where T_i () is the i-th Chebyshev polynomial)Now, I have programs to calculate the all the one-dymensional coefficients(all taken from numerical recipees).How would I put them together to get the 3-dimensional coefficientsa_{i,j,k}, b_{l} from my f[i,j,k] matrix?Had there been no RATIONAL component, I guess it would have been trivialsince everything here is orthogonal, but if I understand my numericalrecipees correctly, this is not the case with the RC algorithm.(And I absolutely NEED the RC for the first function)Is there any algorithm floating around for this problem, or has someonedealt with it before and can share the === Multi-dymensional rational Chebyshev approximation problem >I have a function f(x,y,z) which I need to approximate through numerical >methods, starting from a discrete lattice f[i,j,k]. >f(x, constant y, constant z) would best me modelled by a rational Chebyshev >approximation >f(constant x, y, constant z) can be done by simple Chebyshev polynomials >f(constant x, constant y, z) is even angular, >and so is optimally described through >a series of Cosine Fourier coefficients (found by FFTs or something) >Therefore, the best approximation choice is > Sum_{i,j,k} a_{i,j,k} T_i (x) T_j (y) Cos(kz) >f(x,y,z) ~ ------------------------------------------------- > Sum_{l} b_{l} T_l (x) >(Where T_i () is the i-th Chebyshev polynomial) >Now, I have programs to calculate the all the one-dymensional coefficients >(all taken from numerical recipees). >How would I put them together to get the 3-dimensional coefficients >a_{i,j,k}, b_{l} from my f[i,j,k] matrix? >Had there been no RATIONAL component, I guess it would have been trivial >since everything here is orthogonal, but if I understand my numerical >recipees correctly, this is not the case with the RC algorithm. >(And I absolutely NEED the RC for the first function) >Is there any algorithm floating around for this problem, or has someone >dealt with for any help >G. Torrieri >torrieri[at]physics.arizona.edunonlinear l-infinity fit with three independent variables: I know of no such thing . but if you multiply by the denominator denom= Sum_{l} b_{l} T_l (x) -eps*denom <= >f(x,y,z)- Sum_{i,j,k} a_{i,j,k} T_i (x) T_j (y) Cos(kz) <= eps*denom (*)and minimize now the error function F(eps,a,b) def=eps unknown vector of the problem is (eps,a,b)subject to the constraints you obtain by inserting your (x,y,z)-data into (*) then you obtain a nonlinear optimization problem with nonlinear inequalityconstraints. since this is feasible and will have what is called a setof finely grained constraints, a job for andre tits code cfsqp or one ofthe large scale (if there are many coeffs) nlo codes.seehttp://plato.la.asu.edu/topics/problems/ === Multi-dymensional rational Chebyshev approximation problem-> ->I have a function f(x,y,z) which I need to approximate through numerical->methods, starting from a discrete lattice f[i,j,k].->->f(x, constant y, constant z) would best me modelled by a rational Chebyshev ->approximation->->f(constant x, y, constant z) can be done by simple Chebyshev polynomials->->f(constant x, constant y, z) is even angular, ->and so is optimally described through->a series of Cosine Fourier coefficients (found by FFTs or something)->->Therefore, the best approximation choice is->-> Sum_{i,j,k} a_{i,j,k} T_i (x) T_j (y) Cos(kz)->f(x,y,z) ~ --------------------------------------------------> Sum_{l} b_{l} T_l (x)->->torrieri[at]physics.arizona.edu(skip)-> ->nonlinear l-infinity fit with three independent variables: I know of no such ->thing . but if you multiply by the denominator -> denom= Sum_{l} b_{l} T_l (x) -> -> -eps*denom <= >f(x,y,z)- Sum_{i,j,k} a_{i,j,k} T_i (x) T_j (y) Cos(kz) -> <= eps*denom (*)->and minimize now the error function -> -> F(eps,a,b) def=eps unknown vector of the problem is (eps,a,b)->subject to the constraints you obtain by inserting your (x,y,z)-data into (*) (skip) ->peterIf you really need a true L_inf norm approximant, with formal error bounds, then Peter's approach may be the way.If you are lucky, then multiplying through by the denominator, and invoking the various orthogonalityrelationships among the T_n and trig functions might get you somewhere analytically.For a rude, crude estimate, I'm fond of the following -- used it many times in 1 D, but no reasonwhy it won't work in 3D.Let the vector y be your lhs, evalua on the lattice. Let U be the matrix of numerator basis functions, evaluaon the lattice, over the set of (i,j,k). Let V be the equivalent matrix of denominator basis functions, exceptpull out the l=0 (constant term) and set it = 1.So, y = U a /(V b + 1), where the divison is element by element.Or, (y*V) b + y = U a, where the * means multiply the j'th row of V by y_j. Call that matrix W (= y*V)Now, W b + y = U a, or y= U a which can be solved via SVD in the usual manner, for the solution vector (a:b). -W bHere's scilab code for the 1-D case, using polynomials, for coefd,yfit]=ratfit(x,y,nn,nd)// fit polynomial/polynomial to x,y data.// x is column vector// leading coeff of denominator always 1.// use peval(coefn,x)./peval(coefd,x) to evaluatenx=length(x);xnarr=zeros(nx,nn+1)xdarr=zeros(nx,nd)// Numerator matrixfor i=0:nn; xnarr(1:$,i+1)=x^i;end;// Denominatorfor i=1:nd; xdarr(1:$,i) =x^i .* y ;end;// solvecoef=[-xdarr,xnarr]y// evaluate fityfit=[-xdarr,xnarr]*coef//Numerator and denominator === libraryDoes anyone know of a good multiple precision library abailable out therefor Java?Also, for AR, ARMA, MA, ARIMA & GARCH computations, should we use a multipleprecision === after a long investigation on my part that theforthnet server, had set a certain option (mod_gzip) on JPG's & GIF's,the result of which was all my web-pages on Mathematics (and other pageson my site) sometimes failed to load their images. and, my apologies,but occasionally server providers take into account only certain classesof people (notably users of M$ products, which load the imagescorrectly).Users with Netscape Communicator 4.5-4.7x would not have been able toload GIF's &JPG's. Older versions of Navigator seem to have been workingok.The problem has now been fixed, so please if you viewed these usingNetscape Communicator 4.5/4.7x, please reload them.-- Ioannishttp://users.forthnet.gr/ath/jgal/_____________________ ______________________Eventually, _everything_ is links.Just about to build a new system.This system will be used for numerical float-point computation.OS GNU/Linux or OS X(G5).I am willing to pay under 500$ for a processor, however, if performance significantly increases with price, then I might reconcider.It would be interesting to read other tips rela to am looking for advices or links.> Just about to build a new system.> This system will be used for numerical float-point computation.> OS GNU/Linux or OS X(G5).> I am willing to pay under 500$ for a processor, however, if performance> significantly increases with price, then I might reconcider.> It would be interesting to read other tips rela to Taras.Depends in part on how much precision you need. Standard mpu's have80 bit internal representations.Also, as someone no, most systems are cost-optimized so the time ofthe numerical operations is about the same as that for memory fetch/store.in effect, that if you made the math unit infinitely fast you wouldonly halve the execution time. (Roughly.)can do multiplications and adds (with one operand in memory) in 1-3 clockcycles. But this means that the prefetch queue must be full always, andthe branches taken predic correctly (by whatever branch-predictionalgorithm is being used). Having multiple pipelines helps a lot. But asC. Bond notes, this favorable state of affairs will be upset by a computationwhose statistics differ from those used in the branch-prediction algorithm.Thus you really have to test these things yourself.Note that the cpu of choice for a new parallel system being construcat one of the supercomputer centers is the Gameboy II. Have you consideredit at all?-- Julian V. ^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ Science knows only one commandment: contribute to science. -- Bertolt links.> Just about to build a new system.> This system will be used for numerical float-point computation.> OS GNU/Linux or OS X(G5).> I am willing to pay under 500$ for a processor, however, if performance> significantly increases with price, then I might reconcider.> It would be interesting to read lot.> Taras.If performance is critical, you may just have to benchmark it yourself.Overall math performance is affec by cache size, bus architecture,background threads, etc. and not just FPU performance. In some cases youcan rewrite FPU code to avoid access stalls and cache dumps, but you're notlikely to find a compiler which will do that consistently or reliably onall platforms.--There are two things you must never attempt to prove: the unprovable -- andthe === functions for conforming plate bending elementWe are developing a finite element code for plate bending elements. Thenodal degrees of freedom are u, v, w, dw/dx, dw/dy and d2w/dx*dy.The interpolation functions are Hermite functions which impose that w,dw/dx, dw/dy or d2w/dx*dy should be one in one particular node, and thatall other degrees of freedom should be zero. The formulation of thestiffness matrix is based on the traditional variational principles.This works fine for rectangular plate elements, but once the initialelement geometry is distor (for example to a trapezium), the resultsare rather inaccurate. There is about 1 % error on the displacements,but the error on the bending moments is much larger (10 % to 50 %).I do not find any documentation on how this problem is solved incommercial codes, because these distor meshes yield good results inAbaqus.Does anybody have any suggestions ?Wim Van PaepegemE-mail : Wim.VanPaepegem@UGent.be