14575 http://www.wolfram.com/products/applications/excel_link/ Jens > > > is there any way of linking Mathematica with Excel? > > Lu.92s ==== I am getting the same problem (though Mathematica 4.1). Has anyone any ideas? Aron. > >I am facing the problem in starting the link to math kernel from within >Excel. > >Specifically, when I get the message link failed to open when I click >Launch button in the 'Start Mathematica Link' dialog box. I have tried using >Multilink too so that I am able to access the kernel from Mathematica and >Excel simultaneously but I still get the same message. > >I am currently using Mathematica 4.0 and Excel 2002 (Excel XP). The programs >Mathematica and Excel otherwise appear to work fine. I am using Windows XP >home edition on Dell 8100 laptop. I have used the Mathematica Link for MS >Excel for Excel 2000 in my installation as there were no specific files for >Excel 2002. Given that I was able to successfully add the menus within excel >I think the addin should work fine but it does not? > >All help would be sincerely appreciated. > >Sincerely, > >Tahir Sheikh. ==== Download the 2002 file from this page. I use it for Excel 2002 Service Pack 2. http://support.wolfram.com/applicationpacks/excel_link/excelxp.html > > I am getting the same problem (though Mathematica 4.1). > Has anyone any ideas? > > Aron. > > > >I am facing the problem in starting the link to math kernel from within >Excel. > >Specifically, when I get the message link failed to open when I click >Launch button in the 'Start Mathematica Link' dialog box. I have tried > using >Multilink too so that I am able to access the kernel from Mathematica > and >Excel simultaneously but I still get the same message. > >I am currently using Mathematica 4.0 and Excel 2002 (Excel XP). The > programs >Mathematica and Excel otherwise appear to work fine. I am using Windows > XP >home edition on Dell 8100 laptop. I have used the Mathematica Link for > MS >Excel for Excel 2000 in my installation as there were no specific files > for >Excel 2002. Given that I was able to successfully add the menus within > excel >I think the addin should work fine but it does not? > >All help would be sincerely appreciated. > >Sincerely, > >Tahir Sheikh. > > > ==== I did not request any accuracy for f. I set the accuracy of the numerical components of the expression f. You cannot request the accuracy of the result of your computation in Mathematica, you can only set the accuracy of the input and later check what accuracy of the output results form it. In my last message on this topic I tried to explain this in the plainest and simplest way I could think of. There is nothing more left for me to say. I feel like Sisyphus but unlike him I can at least give up! Andrzej Kozlowski >> >> [...] >> >> I would say this is correct and show that SetPrecision is very useful >> indeed. It tells you (what of course you ought to already know in this >> case anyway) that machine precision will not give you a realiable >> answer in this case. If you know your numbers with a great deal of >> accuracy you can get an accurate answer: >> >> In[24]:= >> f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - >> 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 100]; >> a=SetPrecision[77617.,100]; b = SetPrecision[33096.,100]; >> >> >> In[26]:= >> {f, Precision[f]} >> >> Out[26]= >> {-0.82739605994682136814116509547981629199903311578438481991 >> 781484167246798617832`61.2597, 61} >> > > Congratulations! You just requested accuracy of 100 for f and got 61 ( > to convince yourself add Accuracy[f] to In[26]). If In[24] one > replaces SetAccuracy by SetPrecision the result is similar. > > PK > >> Again you can be pretty sure that you got an accurate answer, provided >> of course your original setting of precision was valid. >> >> Honestly, to say that SetPrecision and SetAccuaracy are useless is one >> of the silliest thing I have read on this list in years. >> >> > >> Andrzej Kozlowski >> Yokohama, Japan >> http://www.mimuw.edu.pl/~akoz/ >> http://platon.c.u-tokyo.ac.jp/andrzej/ > > > Andrzej Kozlowski Yokohama, Japan http://www.mimuw.edu.pl/~akoz/ http://platon.c.u-tokyo.ac.jp/andrzej/ Reply-To: ==== Here's a more intuitive method, perhaps: a + b/c + d*(Sqrt[e]/c) == 0 f = %[[1, 3]] %% /. f -> g First@Solve[%, g] f^2 == (g^2 /. %) However, it occurs to me you might want a more general method to collect a radical on one side and then square both sides. If so, here's a clumsy first attempt: expr = a + b/c + d*(Sqrt[e]/c) == 0 f = First@Cases[expr, Power[a_, Rational[b_, c_]], Infinity] power = First@Cases[f, Rational[b_, c_] -> Rational[c, b], Infinity] coefficient = First@Cases[expr, Times[a_, f] -> a, Infinity] Solve[expr /. coefficient f -> g, g][[1, 1]] g^2 == (g^2 /. %) % /. g -> coefficient f DrBob -----Original Message----- DrBob -----Original Message----- want the equation to be written out with terms grouped by powers of x, but I think I can do that part :) I'll be very grateful to anyone who can give me some pointers. Or, at least point me to some tutorial in the Mathematica documentation. I've been looking over the documentation and I found Appendix A.5 in The Mathematica Book, but that doesn't help me. I _need_ some examples. I did find a couple of well-written posts in this newsgroup, but not quite close enough to what I want. Troy. =-=-=-=-=-=-=-=-=-= FYI, here's the expression I'm working with. denom = Sqrt[(B^2 - r^2)^2 + 4*(r^2)*(b^2)] cnu = (2*b^2 - B^2 + r^2)/denom snu = -2*b*Sqrt[B^2 - b^2]/denom sif = 2*r*b/denom cif = (r^2 - B^2)/denom pdr = -Cos[ds]*Sin[q]*(snu*cif + cnu*sif) - Sin[ds]*(cnu*cif - snu*sif) 0 == -(B^2 - b^2)*V^2/(r^2) + (((B*V)^2)/( r^2) - 2*w*b*V*Cos[q]*Cos[ds] + (w* r)^2 - (w*r*pdr)^2)*(Cos[qr])^2 Although I said it's a polynomial in x, it's really a polynomial in b that I'm after. ==== > > On Friday, October 4, 2002, at 06:01 PM, DrBob > >[...] > > I would say this is correct and show that > SetPrecision is very useful > indeed. It tells you (what of course you ought > to already know in this > case anyway) that machine precision will not > give you a realiable > answer in this case. If you know your numbers > with a great deal of > accuracy you can get an accurate answer: > > In[24]:= > f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - > b^6 - > 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 100]; > a=SetPrecision[77617.,100]; b = > SetPrecision[33096.,100]; > > > In[26]:= > {f, Precision[f]} > > Out[26]= > > {-0.82739605994682136814116509547981629199903311578438481991 > 781484167246798617832`61.2597, 61} > > > Congratulations! You just requested accuracy of > 100 for f and got 61 ( > to convince yourself add Accuracy[f] to In[26]). > If In[24] one > replaces SetAccuracy by SetPrecision the result is > similar. > > PK > [...] > > One has (initially) an accuracy of 100 for an > expression that contains > variables. > > In[25]:= Clear[a,b,f] > > In[26]:= f = SetAccuracy[333.75*b^6 + > a^2*(11*a^2*b^2 - b^6 - > 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 100]; > > In[27]:= Accuracy[f] > Out[27]= 100. > > Now we assign values to some indeterminants in f. > > In[28]:= a = SetPrecision[77617.,100]; b = > SetPrecision[33096.,100]; > > In[29]:= {f, Precision[f], Accuracy[f]} > Out[29]= > {-0.8273960599468213681411650954798162919990331157843848199178148, > > 61.2599, 61.3422} > > The precision and accuracy has dropped. This is all > according to > standard numerical analysis regarding cancellation > error. You'll find it > in any textbook on the topic. > Assume that I want accuracy and precision of 100 for f. You advice me to make experiments to find out, what should be the initial precision and accuracy of a and b to reach the requested accuracy and precision for f. Notice, that you cannot just repeat I[26], we saw already what happens. I have to re-type I[24], I[25], I[26], I[27], I[28], and I[29] as many times as needed to get f with accuracy and precision 100. Dan, you simply advocate to do MANUAL WORK that should be done by machine. Let's suppose that in the above example I just want 60 digits not 61. Precisely, I want 60 digits and nothing or zeros afterwards. Let's see if I could use SetAccuracy. In[30]:= SetAccuracy[%, 60] Out[30]= -0.82739605994682136814116509547981629199903311578438481991781 In[31]:= % // FullForm Out[30]//FullForm= -0.827396059946821368141165095479816291999033115784384819917814841672467988` 59.9177 Oops, it did not work (as expected). Let's highlight with mouse the expression in Out[30] and copy to a new cell. Oops, we got -0.827396059946821368141165095479816291999033115784384819917814841672467988` 59.9177 again. Let's change Out[30] to a text cell and then copy. In[31]:= -0.82739605994682136814116509547981629199903311578438481991781 Out[31]= -0.82739605994682136814116509547981629199903311578438481991781 Success? Not so fast. In[32]:= % // FullForm Out[32]//FullForm= -0.8273960599468213681411650954798162919990331157843848199178099999999999986 35 08`59.2041 Dan, is there any simple way to get what I want? As I repeated already number of times, at this stage of the development of computer technology, software should do it for me (!). We both know that this is doable. Some of the textbooks that you just advised me to read describe it. As a developer of Mathematica, tell us why do you consider this to be a bad idea? Peter Kosta > As for what happens when you artificially raise > precision (or accuracy) > of machine numbers far beyond that guaranteed by > their internal > representation, that falls into to category of > garbage in, garbage out. > It is, howoever, valid to use SetPrecision to raise > precision in > (typically iterative) algorithms where significance > arithmetic might be > unduly pessimistic due to incorrect assumptions > about uncorollatedness > of numerical error. Examples of such usage have > appeared in this news > group. > > > Daniel Lichtblau > Wolfram Research __________________________________________________ Do you Yahoo!? http://faith.yahoo.com ==== > Are there any known issues with simpy treating the JLink.jar as a Java > extension as follows? > cp JLink.jar $JAVA_HOME/jre/lib/ext? > According to my understanding of the discussion in the Java Tutorial on > extensions, that should work: It does; starting with M4.2, J/Link 2.0 gets preinstalled and comes with a 1.4 murphee ==== Are there any known issues with simpy treating the JLink.jar as a Java extension as follows? cp JLink.jar $JAVA_HOME/jre/lib/ext? According to my understanding of the discussion in the Java Tutorial on extensions, that should work: http://java.sun.com/docs/books/tutorial/ext/basics/install.html Commants? STH . ==== >Are there any known issues with simpy treating the JLink.jar as a Java >extension as follows? >cp JLink.jar $JAVA_HOME/jre/lib/ext? > >According to my understanding of the discussion in the Java Tutorial on >extensions, that should work: >http://java.sun.com/docs/books/tutorial/ext/basics/install.html > >Commants? You should not do this. Code from the jre/lib/ext directory is trusted, so this poses a security risk from malicious applets. Leave JLink.jar where it lives in the JLink directory. If you want it to be available to all Java programs on your system, add its location to your CLASSPATH environment variable (this is not a security risk, as remote applets cannot load classes from CLASSPATH). Todd Gayley Wolfram Research ==== The key is in using the command Factor with the option Extension: In[1]:= Factor[x^4 + x^3 + x^2 + x + 1, Extension -> {GoldenRatio}] Out[1]= -((-1 - x + GoldenRatio*x - x^2)*(1 + GoldenRatio*x + x^2)) For manual verification you should keep in mind that: 1/GoldenRatio = GoldenRatio - 1 Germ.87n Buitrago ----- Original Message ----- > (x^2 + GoldenRatio x + 1) ( x^2 - 1/GoldenRatio x + 1) > > What instructions do I need to execute to achieve this output? > > -Steve Earth > Harker School > http://www.harker.org/ > ==== Actually, including 1/GoldenRatio in the extension leads to an unnecessarily complicated formula. In this case there is no real need to so, since by definition In[30]:= Unevaluated[1/GoldenRatio==GoldenRatio-1]//FullSimplify Out[30]= True If one really insists on having the answer in the form proposed in Steve's original posting one can simply do: (Collect[#1, x] & ) /@ Factor[x^4 + x^3 + x^2 + x + 1, Extension -> {GoldenRatio}] /. -1 + GoldenRatio -> 1/GoldenRatio (-(-1 + x/GoldenRatio - x^2))*(1 + GoldenRatio*x + x^2) > > In[]:=Factor[x^4 + x^3 + x^2 + x + 1, Extension -> {GoldenRatio, > 1/GoldenRatio}] > Out[]=-((-3 - 2*x + Sqrt[5]*x + GoldenRatio*x - 3*x^2)* > (3 + x + Sqrt[5]*x + GoldenRatio*x + 3*x^2))/9 > > Jens > >> >> Greetings MathGroup, >> >> My name is Steve Earth, and I am a new subscriber to this list and >> also a >> new user of Mathematica; so please forgive this rather simple >> question... >> >> I would like to enter the quartic x^4 + x^3 + x^2 + x + 1 into >> Mathematica >> and have it be able to tell me that it factors into >> >> (x^2 + GoldenRatio x + 1) ( x^2 - 1/GoldenRatio x + 1) >> >> What instructions do I need to execute to achieve this output? >> >> -Steve Earth >> Harker School >> http://www.harker.org/ > > > Andrzej Kozlowski Yokohama, Japan http://www.mimuw.edu.pl/~akoz/ http://platon.c.u-tokyo.ac.jp/andrzej/ ==== I want to apply a function to every k-th element of a long list and add the result to the k+1 element. [Actually k = 3 and I just want to multiply myList[[k]] by a constant (independent of k) and add the result to myList[[k+1]] for every value of k that's divisible by 3.] Is there a way to do this -- or in general to get at every k-th element of a list -- that's faster and more elegant than writing a brute force Do[] loop or using Mod[] operators, and that will take advantage of native List operators, but still not be too recondite? I've been thinking about multiplying a copy of myList by a mask list {0,0,1,0,0,1,..} to generate a masked copy and approaches like that. Better ways??? ==== Take[list, am, n, sa] gives elements m through n in steps of s. > I want to apply a function to every k-th element of a long list and > add the result to the k+1 element. > > [Actually k = 3 and I just want to multiply myList[[k]] by a > constant (independent of k) and add the result to myList[[k+1]] for > every value of k that's divisible by 3.] > > Is there a way to do this -- or in general to get at every k-th > element of a list -- that's faster and more elegant than writing a brute > force Do[] loop or using Mod[] operators, and that will take > advantage of native List operators, but still not be too recondite? > > I've been thinking about multiplying a copy of myList by a mask list > {0,0,1,0,0,1,..} to generate a masked copy and approaches like that. > Better ways??? ==== lst= Range[14] {1,2,3,4,5,6,7,8,9,10,11,12,13,14} A list of positions in lst ( for your purpose Range[1, Length[lst], 3] will do) ps= {3,5,10}; The following applies h to each ps element in lst and adds the result to the following element (lst[[ps+1]]=h/@lst[[ps]]+lst[[ps+1]];lst) {1,2,3,4+h[3],5,6+h[5],7,8,9,10,11+h[10],12,13,14} -- Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 > I want to apply a function to every k-th element of a long list and > add the result to the k+1 element. > > [Actually k = 3 and I just want to multiply myList[[k]] by a > constant (independent of k) and add the result to myList[[k+1]] for > every value of k that's divisible by 3.] > > Is there a way to do this -- or in general to get at every k-th > element of a list -- that's faster and more elegant than writing a brute > force Do[] loop or using Mod[] operators, and that will take > advantage of native List operators, but still not be too recondite? > > I've been thinking about multiplying a copy of myList by a mask list > {0,0,1,0,0,1,..} to generate a masked copy and approaches like that. > Better ways??? > ==== Consider the following approach, whish uses the MapAt command, that is Map with 'mapping-position' control. dummyFun={#,trueFun@#}& list={a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p} spec=Partition[Range[3,Length[list],3],1] MapAt[dummyFun,list,spec] %//Flatten Hope that is what you want, Borut | I want to apply a function to every k-th element of a long list and | add the result to the k+1 element. | | [Actually k = 3 and I just want to multiply myList[[k]] by a | constant (independent of k) and add the result to myList[[k+1]] for | every value of k that's divisible by 3.] | | Is there a way to do this -- or in general to get at every k-th | element of a list -- that's faster and more elegant than writing a brute | force Do[] loop or using Mod[] operators, and that will take | advantage of native List operators, but still not be too recondite? | | I've been thinking about multiplying a copy of myList by a mask list | {0,0,1,0,0,1,..} to generate a masked copy and approaches like that. | Better ways??? | Reply-To: kuska@informatik.uni-leipzig.de ==== something like: With[{k=3}, Flatten[ {#[[2]] + c*#[[1]], #[[3]]} & /@ Partition[lst, k, k, {1, 1}], 1] ] ?? Jens > > I want to apply a function to every k-th element of a long list and > add the result to the k+1 element. > > [Actually k = 3 and I just want to multiply myList[[k]] by a > constant (independent of k) and add the result to myList[[k+1]] for > every value of k that's divisible by 3.] > > Is there a way to do this -- or in general to get at every k-th > element of a list -- that's faster and more elegant than writing a brute > force Do[] loop or using Mod[] operators, and that will take > advantage of native List operators, but still not be too recondite? > > I've been thinking about multiplying a copy of myList by a mask list > {0,0,1,0,0,1,..} to generate a masked copy and approaches like that. > Better ways??? Reply-To: ==== Daniel, >>The precision/accuracy tracking mechanism will generally let you know, in some fashion, that you have no trustworthy digits. But it is up to the user to check that sort of thing. In this case Mathematica did NOT let us know, in any fashion, that we had no trustworthy digits. Precision and Accuracy outputs were completely misleading. (16 and -5 respectively.) Even Andrzej Kozlowski, who's adept in Mathematica, thought that would be meaningful, and never came up with a better way to check (other than using infinite precision for numbers that probably aren't known that exactly). Peter Kosta demonstrated that he could get a completely erroneous answer with Infinite precision. I blame the problem primarily, and I don't think there's any way to make the answer meaningful. That's not Mathematica's fault at all, and users need to be aware of that old maxim: garbage in, garbage out. comes up with a 22-digit result, it doesn't take much sophistication to realize the answer can't have 16-digit precision. Here's an even more extreme result: f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 50]; a = 77617.; b = 33096.; f Precision[f] -1.180591620717411303424`71.0721*^21 71 71.0721 digits of precision? I don't think so!! We can do the following instead: x = Interval[333.75]; y = Interval[5.5]; a = Interval[77617.]; b = Interval[33096.]; x*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + y*b^8 + a/(2*b) Interval[{-4.486248158726164*^22, 4.2501298345826815*^22}] and that looks like the right answer, finally!! I like that method. However, that doesn't change the fact that Accuracy, Precision, and SetAccuracy appear to be completely useless. I haven't seen an example in which they did what anyone (but you) thought they should do. Bobby -----Original Message----- > you're not aware there's a problem, it lets you go on your merry way, > working with noise. > > Bobby Mathematica is not a mind reader. But the evaluation sequence, while complicated, is reasonably well documented. If you perform machine arithmetic, or for that matter significance arithmetic, and there is massive cancellation error, no use of SetAccuracy after the fact will fix it. The precision/accuracy tracking mechanism will generally let you know, in some fashion, that you have no trustworthy digits. But it is up to the user to check that sort of thing. It is not obvious to me what sort of error the software might notice to report. If you have a concise example of input, and expected output, I can look further. I've not seen anything in this thread that struck me as a failure of the software to warn the user, but maybe I missed something. Daniel ==== GentleBeings I have a straightforward implementation of successive approximations but I cannot seem to froce the code to find the correct solution when I have trig or exponentials involved. Can the assembled wisdom point to straghtforward fixes I know FindRoot works the object is to teach programming and successive approx, tho. kenf Below is the code Clear[f, g, gi, lim, r, rr, fr, gir, a, b, c, d, conv]; Plot[{x * ((x + 3)), 10*Sin[x]}, {x, 0.01, 2.4}, PlotStyle -> {{RGBColor[1, 0, 0], Thickness[ .006]}, {RGBColor[0, 0, 1], Thickness[ .006]}} ]; rr = FindRoot[x * ((x + 3)) == 10*Sin[x], {x, 2, 0.01, 2.4}]; f[a_] := a * ((a + 3)) /; a > 0; g[b_] := 10. * Sin[b] /; b > 0; gi[c_] := ArcSin[0.1*c] /; c > 0; Print[Actual root is , rr]; lim = 10; r = 2.0; conv = 10^-4; For[i = 1, i < lim, i++, { fr = f[r]; gir = gi[fr]; d = Abs[N[gir] - r]; i If[d < conv, Break[]]; r = gir; Print[The value of x = , r, found after , i, iterations,, with a tolerence , d, n] } ] Print[The value of x = , r, found after , i, iterations,, with a tolerence , d, n] Every man, woman and responsible child has an unalienable individual, civil, Constitutional and human right to obtain, own, and carry, openly or concealed, any weapon -- rifle, shotgun, handgun, machine gun, anything -- any time, any place, without asking anyone's permission. L. Neil Smith ==== > I want to apply a function to every k-th element of a long list and > add the result to the k+1 element. > > [Actually k = 3 and I just want to multiply myList[[k]] by a > constant (independent of k) and add the result to myList[[k+1]] for > every value of k that's divisible by 3.] > > Is there a way to do this -- or in general to get at every k-th > element of a list -- that's faster and more elegant than writing a > brute > force Do[] loop or using Mod[] operators, and that will take > advantage of native List operators, but still not be too recondite? > > I've been thinking about multiplying a copy of myList by a mask > list > {0,0,1,0,0,1,..} to generate a masked copy and approaches like that. > Better ways??? > > Here is a generalization of what you've asked: f[l_List,c_,k_Integer,p_Integer]:=Flatten[Block[{r=#[[k]] c},Join[Take[#,k-1],{r,#[[k+1]]+ r},Drop[#,k+1]]]&/@Partition[l,p]]/;(Mod[Length[l],p]==0&&k (Expand[#1] == Expand[#2] & )] {(1 + (1/2)*(1 - Sqrt[5])*x + x^2)* (1 + (1/2)*(1 + Sqrt[5])*x + x^2)} Andrzej Kozlowski Toyama International University JAPAN http://sigma.tuins.ac.jp/~andrzej/ > Steve > The notebook given after NOTEBOOK below contains functions for > factoring and > partial fractioning. > Here is an application to your problem: the first stage avoids our > needing > to know anything about the answer. > > fc=FactorR[x^4+x^3+x^2+x+1,x] > > (1 - (1/2)*(-1 - Sqrt[5])*x + x^2)* > (1 - (1/2)*(-1 + Sqrt[5])*x + x^2) > > Now we need to get rid of Sqrt[5] in terms of GoldenRatio. > This is rather messy: > > Simplify/@(fc/. Sqrt[5][Rule]2 GoldenRatio-1) > > (1 + x - GoldenRatio*x + x^2)*(1 + GoldenRatio*x + x^2) > > Simplify/@(%/.-GoldenRatio[Rule] 1/GoldenRatio -1) > > (1 + x/GoldenRatio + x^2)*(1 + GoldenRatio*x + x^2) > > > Another example > > PartialFractionsR[(1 + x)x/(1 - 3*x + x^2), x] > > 1 - (2*(-1 + 4*x))/((3 + Sqrt[5] - 2*x)*(-3 + 2*x)) + > (2*(-1 + 4*x))/((-3 + 2*x)*(-3 + Sqrt[5] + 2*x)) > > Simplify[%] > > (x*(1 + x))/(1 - 3*x + x^2) > > NOTEBOOK: to make a notebook from the following, copy from the next > line to > the line preceding XXX and paste into a new Mathematica notebook. > > Notebook[{ > > Cell[CellGroupData[{ > Cell[Factors and PartialFractions, Subtitle], > > Cell[Allan Hayes, 16 August 2001, Text], > > Cell[< > Here are some functions for factoring and expressing in partial > fractions over the reals and over the complex numbers. > >, Text], > > Cell[BoxData[ > (Quit)], Input], > > Cell[BoxData[{ > (Off[General::spell1, General::spell]), n, > ((FactorC::usage = polynomial in x with complex coefficients, gives its factorization > over the complex numbers.n > The output may include Root objects which may be evaluated with > ToRadicals or N.>;)n), n, > ((FactorR::usage = polynomial in x with real coefficients, gives its factorization over > the reals.n > The output may include Root objects which may be evaluated with > ToRadicals or N.>;)n), n, > ((PartialFractionsC::usage = where ratl is a rational in x with complex coefficients, gives its > factorization over the complex numbers.n > The output may include Root objects which may be evaluated with > ToRadicals or N.>;)n), n, > ((PartialFractionsR::usage = where ratl is a rational in x with real coefficients, gives its > factorization over the real numbers.n > The output may include Root objects which may be evaluated with > ToRadicals or N.>;)), n, > (On[General::spell1, General::spell])}], Input, > InitializationCell->True], > > Cell[TextData[{ > FactorC[p_, x_] := , > StyleBox[(*over complex numbers*), > FontFamily->Arial, > FontWeight->Plain], > nTimes @@ Cases[Roots[p == 0, x, n Cubics -> False], u_ == > v_ -> x - v]n nFactorR[p_, x_] := , > StyleBox[(*over reals, coefficients must be real*), > FontFamily->Arial, > FontWeight->Plain], > n (Times @@ Join[Cases[#1, u_ == v_ /; Im[v] == 0 :> n x > - v], Cases[#1, u_ == v_ /; Im[v] > 0 :> n x^2 - x*2*Re[v] + > Abs[v]^2]] & )[n Roots[p == 0, x, Cubics -> False]] > }], Input, > InitializationCell->True], > > Cell[TextData[{ > PartialFractionsC[p_, x_] := , > StyleBox[(*over complex numbers*), > FontFamily->Arial, > FontWeight->Plain], > n(#+Apart[#2/FactorC[#3,x]])&@@Flatten[{PolynomialReduce[#,#2], > #2}]&[Numerator[#],Denominator[#]]&[Together[p]]n n > PartialFractionsR[p_, x_] := , > StyleBox[(*over reals, coefficients must be real*), > FontFamily->Arial, > FontWeight->Plain], > n(#+Apart[#2/FactorR[#3,x]])&@@Flatten[{PolynomialReduce[#,#2], > #2}]&[Numerator[#],Denominator[#]]&[Together[p]] > }], Input, > InitializationCell->True], > > Cell[CellGroupData[{ > > Cell[PROGRAMMING NOTES, Subsubsection], > > Cell[TextData[{ > The option , > StyleBox[Cubics->False, > FontFamily->Courier], > is used to keep the roots of cubics in , > StyleBox[Root[....], > FontFamily->Courier], > form. This is better for computation.n, > StyleBox[Re[v], > FontFamily->Courier], > and , > StyleBox[Abs[v]^2, > FontFamily->Courier], > are used rather than , > StyleBox[v+Conjugate[v] , > FontFamily->Courier], > and , > StyleBox[v*Conjugate[v], > FontFamily->Courier], > to prevent , > StyleBox[Apart, > FontFamily->Courier], > from factorising , > StyleBox[x^2 - x*2*Re[v] + Abs[v]^2], > FontFamily->Courier], > back to complex form. > }], Text] > }, Closed]], > > Cell[CellGroupData[{ > > Cell[EXAMPLES, Subsubsection], > > Cell[pol = Expand[(x - 1)*(x + 1)^2*(x^2 + x + 1)^2*(x^2 + 4)]; , > Input], > > Cell[CellGroupData[{ > > Cell[f1 = FactorC[pol, x], Input], > > Cell[BoxData[ > ((((-1) + x)) (((-2) [ImaginaryI] + > x)) ((2 [ImaginaryI] + > x)) ((1 + x))^2 (((((-1)))^(1/3) + x))^2 > (((-(((-1)))^(2/3)) + x))^2)], Output] > }, Open ]], > > Cell[CellGroupData[{ > > Cell[f2 = FactorR[pol, x], Input], > > Cell[BoxData[ > ((((-1) + x)) ((1 + x))^2 ((4 + > x^2)) ((1 + x + x^2))^2)], Output] > }, Open ]], > > Cell[CellGroupData[{ > > Cell[f3 = FactorR[x^3 + x + 1, x], Input], > > Cell[BoxData[ > (((x - Root[1 + #1 + #1^3 &, 1])) ((x^2 - > 2 x Root[(-1) + 2 #1 + 8 #1^3 &, 1] + > Root[(-1) - #1^4 + #1^6 &, 2]^2)))], Output] > }, Open ]], > > Cell[< > Root objects appear because of the option Cubics->False in Roots. > We can sometimes get radical forms, but notice the complication. > >, Text], > > Cell[CellGroupData[{ > > Cell[ToRadicals[f3], Input], > > Cell[BoxData[ > (((((2/(3 (((-9) + @93)))))^(1/3) - ((1/2 > (((-9) + @93))))^(1/3)/3^(2/3) + x)) ((1/3 + > 1/3 ((29/2 - (3 @93)/2))^(1/3) + > 1/3 ((1/2 ((29 + 3 @93))))^(1/3) - > 2 ((((1/2 ((9 + @93))))^(1/3)/(2 > 3^(2/3)) - > 1/(2^(2/3) ((3 ((9 + > @93))))^(1/3)))) x + x^2)))], Output] > }, Open ]], > > Cell[Inexact forms can be found, from f3 :, Text], > > Cell[CellGroupData[{ > > Cell[N[f3], Input], > > Cell[BoxData[ > (((((0.6823278038280193`)([InvisibleSpace])) + > x)) ((((1.4655712318767682`)([InvisibleSpace])) - > 0.6823278038280193` x + x^2)))], Output] > }, Open ]], > > Cell[or directly, Text], > > Cell[CellGroupData[{ > > Cell[f3 = FactorR[x^3 + x + 1//N, x], Input], > > Cell[BoxData[ > (((((0.6823278038280193`)([InvisibleSpace])) + > x)) ((((1.4655712318767682`)([InvisibleSpace])) - > 0.6823278038280193` x + x^2)))], Output] > }, Open ]], > > Cell[Partial fractions, Text], > > Cell[CellGroupData[{ > > Cell[pf1 = PartialFractionsR[(2 + x)/pol, x], Input], > > Cell[BoxData[ > (1/(60 (((-1) + x))) - 1/(10 ((1 + x))^2) - > 39/(100 ((1 + x))) + ((-54) - 31 x)/(4225 ((4 + > x^2))) + ((-1) + 3 x)/(13 ((1 + x + x^2))^2) + (44 + > 193 x)/(507 ((1 + x + x^2))))], Output] > }, Open ]], > > Cell[CellGroupData[{ > > Cell[pf2 = PartialFractionsR[(1 + x)x/(1 - 3*x + x^2), x], > Input], > > Cell[< > 1 - (2*(-1 + 4*x))/((3 + Sqrt[5] - 2*x)*(-3 + 2*x)) + > (2*(-1 + 4*x))/((-3 + 2*x)*(-3 + Sqrt[5] + 2*x)) > >, Output] > }, Open ]], > > Cell[CellGroupData[{ > > Cell[BoxData[ > (Simplify[%])], Input], > > Cell[(x*(1 + x))/(1 - 3*x + x^2), Output] > }, Open ]], > > Cell[Partial fractions will often involve Root objects , Text], > > Cell[CellGroupData[{ > > Cell[pf3 = PartialFractionsR[(1 + x)/(x^3 - x + 1), x], Input], > > Cell[BoxData[ > (((1 + > Root[1 - #1 + #1^3 &, > 1]))/((((x - > Root[1 - #1 + #1^3 &, > 1])) ((Root[1 - #1 + #1^3 &, 1]^2 - > 2 Root[1 - #1 + #1^3 &, > 1] Root[(-1) - 2 #1 + 8 #1^3 &, 1] + > Root[(-1) + #1^4 + #1^6 &, 2]^2)))) + ((x + > Root[1 - #1 + #1^3 &, 1] + > x Root[1 - #1 + #1^3 &, 1] - > 2 Root[(-1) - 2 #1 + 8 #1^3 &, 1] - > Root[(-1) + #1^4 + #1^6 &, 2]^2))/(((((-x^2) + > 2 x Root[(-1) - 2 #1 + 8 #1^3 &, 1] - > Root[(-1) + #1^4 + #1^6 &, 2]^2)) ((Root[1 - > #1 + #1^3 &, 1]^2 - > 2 Root[1 - #1 + #1^3 &, > 1] Root[(-1) - 2 #1 + 8 #1^3 &, 1] + > Root[(-1) + #1^4 + #1^6 &, 2]^2)))))], > Output] > }, Open ]], > > Cell[This can in fact be put in radical form:, Text], > > Cell[CellGroupData[{ > > Cell[ToRadicals[pf3], Input], > > Cell[BoxData[ > (((1 - ((2/(3 ((9 - @69)))))^(1/3) - ((1/2 ((9 > - @69))))^(1/3)/3^(2/3)))/(((((-(1/3)) + > 1/3 ((25/2 - (3 @69)/2))^(1/3) + > 1/3 ((1/2 ((25 + 3 @69))))^(1/3) + > (((-((2/(3 ((9 - @69)))))^(1/3)) - ((1/2 ((9 - > @69))))^(1/3)/3^(2/3)))^2 - > 2 (((-((2/(3 ((9 - @69)))))^(1/3)) - > ((1/2 ((9 - @69))))^(1/3)/3^(2/3))) ((1/24 ((864 > - 96 @69))^(1/3) + ((1/2 ((9 + @69))))^(1/3)/(2 > 3^(2/3)))))) ((((2/(3 ((9 - @69)))))^(1/3) + ((1 > /2 ((9 - @69))))^(1/3)/3^(2/3) + x)))) + ((1/3 - > 1/3 ((25/2 - (3 @69)/2))^(1/3) - ((2/(3 > ((9 - @69)))))^(1/3) - ((1/2 ((9 - @69))))^(1/3)/3 > ^(2/3) - 1/3 ((1/2 ((25 + 3 @69))))^(1/3) - > 2 ((1/24 ((864 - 96 @69))^(1/3) + ((1/2 > ((9 + @69))))^(1/3)/(2 3^(2/3)))) + > x + (((-((2/(3 ((9 - @69)))))^(1/3)) - > ((1/2 ((9 - @69))))^(1/3)/3^(2/3))) x))/(((((-(1 > /3)) + 1/3 ((25/2 - (3 @69)/2))^(1/3) + > 1/3 ((1/2 ((25 + 3 @69))))^(1/3) + > (((-((2/(3 ((9 - @69)))))^(1/3)) - ((1/2 ((9 - > @69))))^(1/3)/3^(2/3)))^2 - > 2 (((-((2/(3 ((9 - @69)))))^(1/3)) - > ((1/2 ((9 - @69))))^(1/3)/3^(2/3))) ((1/24 ((864 > - 96 @69))^(1/3) + ((1/2 ((9 + @69))))^(1/3)/(2 > 3^(2/3)))))) ((1/3 - > 1/3 ((25/2 - (3 @69)/2))^(1/3) - > 1/3 ((1/2 ((25 + 3 @69))))^(1/3) + > 2 ((1/24 ((864 - 96 @69))^(1/3) + ((1/2 > ((9 + @69))))^(1/3)/(2 3^(2/3)))) x - > x^2)))))], Output] > }, Closed]], > > Cell[We could have found the inexact form directly., Text], > > Cell[CellGroupData[{ > > Cell[BoxData[ > (PartialFractionsR[((1 + x))/((x^3 - x + 1)) // N, > x])], Input], > > Cell[BoxData[ > ((-(0.07614206365252976`/(((1.324717957244746`)( > [InvisibleSpace])) + > 1.` x))) + (((0.7982664819556426`)( > [InvisibleSpace])) + 0.07614206365252976` > x)/(((0.754877666246693`)([InvisibleSpace])) - > 1.324717957244746` x + 1.` x^2))], Output] > }, Open ]] > }, Closed]] > }, Open ]] > }, > ScreenRectangle->{{0, 1024}, {0, 709}}, > AutoGeneratedPackage->None, > WindowSize->{534, 628}, > WindowMargins->{{199, Automatic}, {0, Automatic}}, > ShowCellLabel->False, > StyleDefinitions -> Default.nb > ] > > XXX > -- > Allan > > --------------------- > Allan Hayes > Mathematica Training and Consulting > Leicester UK > www.haystack.demon.co.uk > hay@haystack.demon.co.uk > Voice: +44 (0)116 271 4198 > > >> Greetings MathGroup, >> >> My name is Steve Earth, and I am a new subscriber to this list and >> also a >> new user of Mathematica; so please forgive this rather simple >> question... >> >> I would like to enter the quartic x^4 + x^3 + x^2 + x + 1 into >> Mathematica > >> and have it be able to tell me that it factors into >> >> (x^2 + GoldenRatio x + 1) ( x^2 - 1/GoldenRatio x + 1) >> >> What instructions do I need to execute to achieve this output? >> >> -Steve Earth >> Harker School >> http://www.harker.org/ >> > > > > > > > > ==== In my earlier posting I used Union and SameTest to replace two equivalent answers (arising form the symmetry of the equation) by a single one. However, the way I did, while givign the right answer, it makes little logical sense since in general applying Expand would make any factorizations the same leaving us always with just a single one. WIthout using Union at all we get: (a + b*x + x^2)*(c + d*x + x^2) /. Select[SolveAlways[x^4 + x^3 + x^2 + x + 1 == (a + b*x + x^2)*(c + d*x + x^2), x], FreeQ[#1, I] & ] {(1 + (1/2 - Sqrt[5]/2)*x + x^2)* (1 + (1/2)*(1 + Sqrt[5])*x + x^2), (1 + (1/2)*(1 - Sqrt[5])*x + x^2)* (1 + (1/2)*(1 + Sqrt[5])*x + x^2)} Since having two identical answers differing only in the way they are written out is a bit of a nuisance, a way to get rid of one of them which does not suffer from the illogicality of my first approach is: Union[(a + b*x + x^2)*(c + d*x + x^2) /. Select[SolveAlways[x^4 + x^3 + x^2 + x + 1 == (a + b*x + x^2)*(c + d*x + x^2), x], FreeQ[#1, I] & ], SameTest -> (Simplify[First[#1] == First[#2]] & )] {(1 + (1/2)*(1 - Sqrt[5])*x + x^2)* (1 + (1/2)*(1 + Sqrt[5])*x + x^2)} > There is an equivalent approach that will give the answer without > knowing it in advance. All we need to know is that any quartic can be > factored over the reals as a product of two quadratics, so: > > > Union[(a + b*x + x^2)*(c + d*x + x^2) /. > Select[SolveAlways[x^4 + x^3 + x^2 + x + 1 == > (a + b*x + x^2)*(c + d*x + x^2), x], FreeQ[#1, I] & ], > SameTest -> (Expand[#1] == Expand[#2] & )] > > > {(1 + (1/2)*(1 - Sqrt[5])*x + x^2)* > (1 + (1/2)*(1 + Sqrt[5])*x + x^2)} > > > Andrzej Kozlowski > Toyama International University > JAPAN > http://sigma.tuins.ac.jp/~andrzej/ > > > > > >> Steve >> The notebook given after NOTEBOOK below contains functions for >> factoring and >> partial fractioning. >> Here is an application to your problem: the first stage avoids our >> needing >> to know anything about the answer. >> >> fc=FactorR[x^4+x^3+x^2+x+1,x] >> >> (1 - (1/2)*(-1 - Sqrt[5])*x + x^2)* >> (1 - (1/2)*(-1 + Sqrt[5])*x + x^2) >> >> Now we need to get rid of Sqrt[5] in terms of GoldenRatio. >> This is rather messy: >> >> Simplify/@(fc/. Sqrt[5][Rule]2 GoldenRatio-1) >> >> (1 + x - GoldenRatio*x + x^2)*(1 + GoldenRatio*x + x^2) >> >> Simplify/@(%/.-GoldenRatio[Rule] 1/GoldenRatio -1) >> >> (1 + x/GoldenRatio + x^2)*(1 + GoldenRatio*x + x^2) >> >> >> Another example >> >> PartialFractionsR[(1 + x)x/(1 - 3*x + x^2), x] >> >> 1 - (2*(-1 + 4*x))/((3 + Sqrt[5] - 2*x)*(-3 + 2*x)) + >> (2*(-1 + 4*x))/((-3 + 2*x)*(-3 + Sqrt[5] + 2*x)) >> >> Simplify[%] >> >> (x*(1 + x))/(1 - 3*x + x^2) >> >> NOTEBOOK: to make a notebook from the following, copy from the next >> line to >> the line preceding XXX and paste into a new Mathematica notebook. >> >> Notebook[{ >> >> Cell[CellGroupData[{ >> Cell[Factors and PartialFractions, Subtitle], >> >> Cell[Allan Hayes, 16 August 2001, Text], >> >> Cell[< >> Here are some functions for factoring and expressing in partial >> fractions over the reals and over the complex numbers. >> >, Text], >> >> Cell[BoxData[ >> (Quit)], Input], >> >> Cell[BoxData[{ >> (Off[General::spell1, General::spell]), n, >> ((FactorC::usage = > polynomial in x with complex coefficients, gives its factorization >> over the complex numbers.n >> The output may include Root objects which may be evaluated with >> ToRadicals or N.>;)n), n, >> ((FactorR::usage = > polynomial in x with real coefficients, gives its factorization over >> the reals.n >> The output may include Root objects which may be evaluated with >> ToRadicals or N.>;)n), n, >> ((PartialFractionsC::usage = > where ratl is a rational in x with complex coefficients, gives its >> factorization over the complex numbers.n >> The output may include Root objects which may be evaluated with >> ToRadicals or N.>;)n), n, >> ((PartialFractionsR::usage = > where ratl is a rational in x with real coefficients, gives its >> factorization over the real numbers.n >> The output may include Root objects which may be evaluated with >> ToRadicals or N.>;)), n, >> (On[General::spell1, General::spell])}], Input, >> InitializationCell->True], >> >> Cell[TextData[{ >> FactorC[p_, x_] := , >> StyleBox[(*over complex numbers*), >> FontFamily->Arial, >> FontWeight->Plain], >> nTimes @@ Cases[Roots[p == 0, x, n Cubics -> False], u_ == >> v_ -> x - v]n nFactorR[p_, x_] := , >> StyleBox[(*over reals, coefficients must be real*), >> FontFamily->Arial, >> FontWeight->Plain], >> n (Times @@ Join[Cases[#1, u_ == v_ /; Im[v] == 0 :> n x >> >> - v], Cases[#1, u_ == v_ /; Im[v] > 0 :> n x^2 - x*2*Re[v] + >> Abs[v]^2]] & )[n Roots[p == 0, x, Cubics -> False]] >> }], Input, >> InitializationCell->True], >> >> Cell[TextData[{ >> PartialFractionsC[p_, x_] := , >> StyleBox[(*over complex numbers*), >> FontFamily->Arial, >> FontWeight->Plain], >> n(#+Apart[#2/FactorC[#3,x]])&@@Flatten[{PolynomialReduce[#,#2], >> #2}]&[Numerator[#],Denominator[#]]&[Together[p]]n n >> PartialFractionsR[p_, x_] := , >> StyleBox[(*over reals, coefficients must be real*), >> FontFamily->Arial, >> FontWeight->Plain], >> n(#+Apart[#2/FactorR[#3,x]])&@@Flatten[{PolynomialReduce[#,#2], >> #2}]&[Numerator[#],Denominator[#]]&[Together[p]] >> }], Input, >> InitializationCell->True], >> >> Cell[CellGroupData[{ >> >> Cell[PROGRAMMING NOTES, Subsubsection], >> >> Cell[TextData[{ >> The option , >> StyleBox[Cubics->False, >> FontFamily->Courier], >> is used to keep the roots of cubics in , >> StyleBox[Root[....], >> FontFamily->Courier], >> form. This is better for computation.n, >> StyleBox[Re[v], >> FontFamily->Courier], >> and , >> StyleBox[Abs[v]^2, >> FontFamily->Courier], >> are used rather than , >> StyleBox[v+Conjugate[v] , >> FontFamily->Courier], >> and , >> StyleBox[v*Conjugate[v], >> FontFamily->Courier], >> to prevent , >> StyleBox[Apart, >> FontFamily->Courier], >> from factorising , >> StyleBox[x^2 - x*2*Re[v] + Abs[v]^2], >> FontFamily->Courier], >> back to complex form. >> }], Text] >> }, Closed]], >> >> Cell[CellGroupData[{ >> >> Cell[EXAMPLES, Subsubsection], >> >> Cell[pol = Expand[(x - 1)*(x + 1)^2*(x^2 + x + 1)^2*(x^2 + 4)]; , >> Input], >> >> Cell[CellGroupData[{ >> >> Cell[f1 = FactorC[pol, x], Input], >> >> Cell[BoxData[ >> ((((-1) + x)) (((-2) [ImaginaryI] + >> x)) ((2 [ImaginaryI] + >> x)) ((1 + x))^2 (((((-1)))^(1/3) + x))^2 >> (((-(((-1)))^(2/3)) + x))^2)], Output] >> }, Open ]], >> >> Cell[CellGroupData[{ >> >> Cell[f2 = FactorR[pol, x], Input], >> >> Cell[BoxData[ >> ((((-1) + x)) ((1 + x))^2 ((4 + >> x^2)) ((1 + x + x^2))^2)], Output] >> }, Open ]], >> >> Cell[CellGroupData[{ >> >> Cell[f3 = FactorR[x^3 + x + 1, x], Input], >> >> Cell[BoxData[ >> (((x - Root[1 + #1 + #1^3 &, 1])) ((x^2 - >> 2 x Root[(-1) + 2 #1 + 8 #1^3 &, 1] + >> Root[(-1) - #1^4 + #1^6 &, 2]^2)))], Output] >> }, Open ]], >> >> Cell[< >> Root objects appear because of the option Cubics->False in Roots. >> We can sometimes get radical forms, but notice the complication. >> >, Text], >> >> Cell[CellGroupData[{ >> >> Cell[ToRadicals[f3], Input], >> >> Cell[BoxData[ >> (((((2/(3 (((-9) + @93)))))^(1/3) - ((1/2 >> (((-9) + @93))))^(1/3)/3^(2/3) + x)) ((1/3 + >> 1/3 ((29/2 - (3 @93)/2))^(1/3) + >> 1/3 ((1/2 ((29 + 3 @93))))^(1/3) - >> 2 ((((1/2 ((9 + @93))))^(1/3)/(2 >> 3^(2/3)) - >> 1/(2^(2/3) ((3 ((9 + >> @93))))^(1/3)))) x + x^2)))], Output] >> }, Open ]], >> >> Cell[Inexact forms can be found, from f3 :, Text], >> >> Cell[CellGroupData[{ >> >> Cell[N[f3], Input], >> >> Cell[BoxData[ >> (((((0.6823278038280193`)([InvisibleSpace])) + >> x)) ((((1.4655712318767682`)([InvisibleSpace])) - >> 0.6823278038280193` x + x^2)))], Output] >> }, Open ]], >> >> Cell[or directly, Text], >> >> Cell[CellGroupData[{ >> >> Cell[f3 = FactorR[x^3 + x + 1//N, x], Input], >> >> Cell[BoxData[ >> (((((0.6823278038280193`)([InvisibleSpace])) + >> x)) ((((1.4655712318767682`)([InvisibleSpace])) - >> 0.6823278038280193` x + x^2)))], Output] >> }, Open ]], >> >> Cell[Partial fractions, Text], >> >> Cell[CellGroupData[{ >> >> Cell[pf1 = PartialFractionsR[(2 + x)/pol, x], Input], >> >> Cell[BoxData[ >> (1/(60 (((-1) + x))) - 1/(10 ((1 + x))^2) - >> 39/(100 ((1 + x))) + ((-54) - 31 x)/(4225 ((4 + >> x^2))) + ((-1) + 3 x)/(13 ((1 + x + x^2))^2) + (44 + >> >> 193 x)/(507 ((1 + x + x^2))))], Output] >> }, Open ]], >> >> Cell[CellGroupData[{ >> >> Cell[pf2 = PartialFractionsR[(1 + x)x/(1 - 3*x + x^2), x], >> Input], >> >> Cell[< >> 1 - (2*(-1 + 4*x))/((3 + Sqrt[5] - 2*x)*(-3 + 2*x)) + >> (2*(-1 + 4*x))/((-3 + 2*x)*(-3 + Sqrt[5] + 2*x)) >> >, Output] >> }, Open ]], >> >> Cell[CellGroupData[{ >> >> Cell[BoxData[ >> (Simplify[%])], Input], >> >> Cell[(x*(1 + x))/(1 - 3*x + x^2), Output] >> }, Open ]], >> >> Cell[Partial fractions will often involve Root objects , Text], >> >> Cell[CellGroupData[{ >> >> Cell[pf3 = PartialFractionsR[(1 + x)/(x^3 - x + 1), x], Input], >> >> Cell[BoxData[ >> (((1 + >> Root[1 - #1 + #1^3 &, >> 1]))/((((x - >> Root[1 - #1 + #1^3 &, >> 1])) ((Root[1 - #1 + #1^3 &, 1]^2 - >> 2 Root[1 - #1 + #1^3 &, >> 1] Root[(-1) - 2 #1 + 8 #1^3 &, 1] + >> Root[(-1) + #1^4 + #1^6 &, 2]^2)))) + ((x + >> Root[1 - #1 + #1^3 &, 1] + >> x Root[1 - #1 + #1^3 &, 1] - >> 2 Root[(-1) - 2 #1 + 8 #1^3 &, 1] - >> Root[(-1) + #1^4 + #1^6 &, 2]^2))/(((((-x^2) + >> 2 x Root[(-1) - 2 #1 + 8 #1^3 &, 1] - >> Root[(-1) + #1^4 + #1^6 &, 2]^2)) ((Root[1 - >> #1 + #1^3 &, 1]^2 - >> 2 Root[1 - #1 + #1^3 &, >> 1] Root[(-1) - 2 #1 + 8 #1^3 &, 1] + >> Root[(-1) + #1^4 + #1^6 &, 2]^2)))))], >> Output] >> }, Open ]], >> >> Cell[This can in fact be put in radical form:, Text], >> >> Cell[CellGroupData[{ >> >> Cell[ToRadicals[pf3], Input], >> >> Cell[BoxData[ >> (((1 - ((2/(3 ((9 - @69)))))^(1/3) - ((1/2 ((9 >> - @69))))^(1/3)/3^(2/3)))/(((((-(1/3)) + >> 1/3 ((25/2 - (3 @69)/2))^(1/3) + >> 1/3 ((1/2 ((25 + 3 @69))))^(1/3) + >> (((-((2/(3 ((9 - @69)))))^(1/3)) - ((1/2 ((9 - >> @69))))^(1/3)/3^(2/3)))^2 - >> 2 (((-((2/(3 ((9 - @69)))))^(1/3)) - >> ((1/2 ((9 - @69))))^(1/3)/3^(2/3))) ((1/24 ((864 >> - 96 @69))^(1/3) + ((1/2 ((9 + @69))))^(1/3)/(2 >> 3^(2/3)))))) ((((2/(3 ((9 - @69)))))^(1/3) + ((1 >> /2 ((9 - @69))))^(1/3)/3^(2/3) + x)))) + ((1/3 - >> 1/3 ((25/2 - (3 @69)/2))^(1/3) - ((2/(3 >> ((9 - @69)))))^(1/3) - ((1/2 ((9 - @69))))^(1/3)/3 >> ^(2/3) - 1/3 ((1/2 ((25 + 3 @69))))^(1/3) - >> 2 ((1/24 ((864 - 96 @69))^(1/3) + ((1/2 >> ((9 + @69))))^(1/3)/(2 3^(2/3)))) + >> x + (((-((2/(3 ((9 - @69)))))^(1/3)) - >> ((1/2 ((9 - @69))))^(1/3)/3^(2/3))) x))/(((((-(1 >> /3)) + 1/3 ((25/2 - (3 @69)/2))^(1/3) + >> 1/3 ((1/2 ((25 + 3 @69))))^(1/3) + >> (((-((2/(3 ((9 - @69)))))^(1/3)) - ((1/2 ((9 - >> @69))))^(1/3)/3^(2/3)))^2 - >> 2 (((-((2/(3 ((9 - @69)))))^(1/3)) - >> ((1/2 ((9 - @69))))^(1/3)/3^(2/3))) ((1/24 ((864 >> - 96 @69))^(1/3) + ((1/2 ((9 + @69))))^(1/3)/(2 >> 3^(2/3)))))) ((1/3 - >> 1/3 ((25/2 - (3 @69)/2))^(1/3) - >> 1/3 ((1/2 ((25 + 3 @69))))^(1/3) + >> 2 ((1/24 ((864 - 96 @69))^(1/3) + ((1/2 >> >> ((9 + @69))))^(1/3)/(2 3^(2/3)))) x - >> x^2)))))], Output] >> }, Closed]], >> >> Cell[We could have found the inexact form directly., Text], >> >> Cell[CellGroupData[{ >> >> Cell[BoxData[ >> (PartialFractionsR[((1 + x))/((x^3 - x + 1)) // N, >> x])], Input], >> >> Cell[BoxData[ >> ((-(0.07614206365252976`/(((1.324717957244746`)( >> [InvisibleSpace])) + >> 1.` x))) + (((0.7982664819556426`)( >> [InvisibleSpace])) + 0.07614206365252976` >> x)/(((0.754877666246693`)([InvisibleSpace])) - >> 1.324717957244746` x + 1.` x^2))], Output] >> }, Open ]] >> }, Closed]] >> }, Open ]] >> }, >> ScreenRectangle->{{0, 1024}, {0, 709}}, >> AutoGeneratedPackage->None, >> WindowSize->{534, 628}, >> WindowMargins->{{199, Automatic}, {0, Automatic}}, >> ShowCellLabel->False, >> StyleDefinitions -> Default.nb >> ] >> >> XXX >> -- >> Allan >> >> --------------------- >> Allan Hayes >> Mathematica Training and Consulting >> Leicester UK >> www.haystack.demon.co.uk >> hay@haystack.demon.co.uk >> Voice: +44 (0)116 271 4198 >> >> > Greetings MathGroup, > > My name is Steve Earth, and I am a new subscriber to this list and > also a > new user of Mathematica; so please forgive this rather simple > question... > > I would like to enter the quartic x^4 + x^3 + x^2 + x + 1 into > Mathematica >> > and have it be able to tell me that it factors into > > (x^2 + GoldenRatio x + 1) ( x^2 - 1/GoldenRatio x + 1) > > What instructions do I need to execute to achieve this output? > > -Steve Earth > Harker School > http://www.harker.org/ > >> >> >> >> >> >> >> >> > > Andrzej Kozlowski Yokohama, Japan http://www.mimuw.edu.pl/~akoz/ http://platon.c.u-tokyo.ac.jp/andrzej/ Reply-To: ==== Let's be realistic. If you want 60 digits of precision, too bad! -- in the real world. There's nothing we can measure that closely. Drug concentrations in clinical trials are generally measured within 15%, for instance. Even machine precision is more than can be realistically expected in any application I can think of. Even getting a satellite to Jupiter probably involves more error in the final result than machine precision. (If not, it's because we rely on ongoing corrections and natural factors that put the satellite where it should be, such as gravity drawing it toward each rendezvous -- not on that kind of precision in propulsion or guidance.) So... unless all numerics in a problem have a theoretical origin, and could be represented in Mathematica as Infinite precision expressions... all this talk of higher-precision computation seems futile. The realistic question is this: given that I have confidence, say, in 6 digits of precision for the inputs of an expression, how many digits can I trust in the end result? Giving inputs MORE precision than they deserve isn't the best way to answer that question. Here are two methods of answering it in Mathematica. One uses bignums and the other uses Intervals. Repetitious trial and error are NOT required either way. BIGNUMS: ClearAll[a, b, f, x, y] f = x*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + y*b^8 + a/(2*b); a = SetPrecision[77617, 6]; b = SetPrecision[33096, 6]; x = SetPrecision[33375/100, 6]; y = SetPrecision[55/10, 6]; InputForm[f] -4.1396`-12.5121*^19 Several previous solutions have set the precision or accuracy of f before giving a and b (and possibly x and y) values. That results in making the exponents imprecise along with all coefficients (not just x and y), which may or may not be what we want. INTERVALS: I'll first digress to figure out what Interval is equivalent to 6-digit precision. You might not actually do this if you like the Interval method, but you have to decide SOMEHOW what Interval width to use. nums = {77617, 33096, 33375/100, 55/10}; (Interval[SetPrecision[#1, 6]] & ) /@ nums /. Interval[{a_, b_}] :> 2*((b - a)/(b + a)); InputForm[%] {3.2209438653903`0.207*^-6, 3.7768914672467`0.2761*^-6, 2.9260299625467`0.1652*^-6, 2.7743252840909`0.1421*^-6} For these numbers, # + Interval[{-1,1}]*#/630000& gets us very close, so I'll use that. The second method therefore is: g = #1 + Interval[{-1, 1}]*(#1/630000) & ; a = g[77617]; b = g[33096]; x = g[333.75]; y = g[5.5]; f = x*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + y*b^8 + a/(2*b) (Min[f] + Max[f])/2 Interval[{-2.136361928005054*^32, 2.1363651195928296*^32}] 1.5957938878075177*^26 Either method shows the answer has no trustworthy digits, but I think the second result is far easier to interpret. Here's another example, using the Sin function, whose derivative is Cos, whose magnitude is bounded by one. The precision of the Sin of a number should be GREATER than the precision of the number itself, especially when Cos is small. a = SetPrecision[Pi/2, 6]; InputForm[Sin[a]] 0.9999999999990905052982256654`11.6078 a = g[N[Pi/2]]; Sin[a] - 1 Interval[{-3.1084024243455137*^-12, 0}] I hope this was worthwhile to someone. DrBob -----Original Message----- >[...] > > I would say this is correct and show that > SetPrecision is very useful > indeed. It tells you (what of course you ought > to already know in this > case anyway) that machine precision will not > give you a realiable > answer in this case. If you know your numbers > with a great deal of > accuracy you can get an accurate answer: > > In[24]:= > f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - > b^6 - > 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 100]; > a=SetPrecision[77617.,100]; b = > SetPrecision[33096.,100]; > > > In[26]:= > {f, Precision[f]} > > Out[26]= > > {-0.82739605994682136814116509547981629199903311578438481991 > 781484167246798617832`61.2597, 61} > > > Congratulations! You just requested accuracy of > 100 for f and got 61 ( > to convince yourself add Accuracy[f] to In[26]). > If In[24] one > replaces SetAccuracy by SetPrecision the result is > similar. > > PK > [...] > > One has (initially) an accuracy of 100 for an > expression that contains > variables. > > In[25]:= Clear[a,b,f] > > In[26]:= f = SetAccuracy[333.75*b^6 + > a^2*(11*a^2*b^2 - b^6 - > 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 100]; > > In[27]:= Accuracy[f] > Out[27]= 100. > > Now we assign values to some indeterminants in f. > > In[28]:= a = SetPrecision[77617.,100]; b = > SetPrecision[33096.,100]; > > In[29]:= {f, Precision[f], Accuracy[f]} > Out[29]= > {-0.8273960599468213681411650954798162919990331157843848199178148, > > 61.2599, 61.3422} > > The precision and accuracy has dropped. This is all > according to > standard numerical analysis regarding cancellation > error. You'll find it > in any textbook on the topic. > Assume that I want accuracy and precision of 100 for f. You advice me to make experiments to find out, what should be the initial precision and accuracy of a and b to reach the requested accuracy and precision for f. Notice, that you cannot just repeat I[26], we saw already what happens. I have to re-type I[24], I[25], I[26], I[27], I[28], and I[29] as many times as needed to get f with accuracy and precision 100. Dan, you simply advocate to do MANUAL WORK that should be done by machine. Let's suppose that in the above example I just want 60 digits not 61. Precisely, I want 60 digits and nothing or zeros afterwards. Let's see if I could use SetAccuracy. In[30]:= SetAccuracy[%, 60] Out[30]= -0.82739605994682136814116509547981629199903311578438481991781 In[31]:= % // FullForm Out[30]//FullForm= -0.827396059946821368141165095479816291999033115784384819917814841672467 988` 59.9177 Oops, it did not work (as expected). Let's highlight with mouse the expression in Out[30] and copy to a new cell. Oops, we got -0.827396059946821368141165095479816291999033115784384819917814841672467 988` 59.9177 again. Let's change Out[30] to a text cell and then copy. In[31]:= -0.82739605994682136814116509547981629199903311578438481991781 Out[31]= -0.82739605994682136814116509547981629199903311578438481991781 Success? Not so fast. In[32]:= % // FullForm Out[32]//FullForm= -0.827396059946821368141165095479816291999033115784384819917809999999999 998635 08`59.2041 Dan, is there any simple way to get what I want? As I repeated already number of times, at this stage of the development of computer technology, software should do it for me (!). We both know that this is doable. Some of the textbooks that you just advised me to read describe it. As a developer of Mathematica, tell us why do you consider this to be a bad idea? Peter Kosta > As for what happens when you artificially raise > precision (or accuracy) > of machine numbers far beyond that guaranteed by > their internal > representation, that falls into to category of > garbage in, garbage out. > It is, howoever, valid to use SetPrecision to raise > precision in > (typically iterative) algorithms where significance > arithmetic might be > unduly pessimistic due to incorrect assumptions > about uncorollatedness > of numerical error. Examples of such usage have > appeared in this news > group. > > > Daniel Lichtblau > Wolfram Research __________________________________________________ Do you Yahoo!? http://faith.yahoo.com ==== Bobby, The example that I gave in my previous posting does not make my point, or it makes it in rather a hidden way, but it does show something interesting about computation with bigfloats. I'll explain what I mean by this and then give an example that does make the point directly. First, the previous example: Sin[#1]*10^#1*Log[1+10^(-#1)]&[15.9] -0.336629 Sin[#1]*10^#1*Log[1+10^(-#1)]&[SetPrecision[15.9,20]] Precision[%] -0.190858581374189370 17.7558 Sin[#1]*10^#1*Log[1+10^(-#1)]&[SetPrecision[15.9,7]] Precision[%] -0.19086 4.70309 (**) It looks as if the internal computations must be to a higher precision than 4 and that they start at SetPrecision[15.9,7]//FullForm 15.9000000000000003553`7 With MaxError = 10^-Accuracy[sp]//FullForm 1.590000000000001`*^-6 Roughly speaking, not more than the first seven digits are asserted to be correct. Now, the new example (taken from Stan Wagon, Programming Tips, Mathematica in Education and Research Volume 7, Number 2, 1988 p50) Clear[x] ser= Normal[Series[Cos[x],{x,0,200}]]; x= 75.0; ser//FullForm -2.7019882604300525`*^15 Probably not reliable. Set precision to 20: x= SetPrecision[75.0, 20]; (a=ser)//FullForm -16928.799183047`1.4688 MaxError= 10^-Accuracy[a] 575.263 Not good enough. Raise the precision to x= SetAccuracy[75.0,40]; (a=ser)//FullForm 1.0807905977573169155`7.7627 MaxError = 10^-Accuracy[a] !(1.3999657487996298`*^-6) That is 1.3999657487996298 10^-6 Good enough Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 ==== Bobby, You note some important points but there are situations in which raising accuracy is a practical necessity. 1) It may be that a calculation (perhaps describing a real world phenonenon) necessitates our raising the precision of inuts so as to work to high precision internally even when the output is not very sensitive to changes in input. Compare the following graphs pts = Table[({#1, Sin[#1]*10^#1*Log[1 + 10^(-#1)]} & )[x], {x, 15., 20., 0.1}]; ListPlot[pts, PlotJoined -> True] pts = Table[({#1, Sin[#1]*10^#1*Log[1 + 10^(-#1)]} & )[ SetAccuracy[x, 20]], {x, 15., 20., 0.1}]; ListPlot[pts, PlotJoined -> True] 2) Another situation in which high precision numbers are needed is when a computation with exact numbers would be slow or might use up all the memory available. We then need to replace the exact numbers with high precision bigfloats which causes the number of digits used internally to be restricted and the maximum error in the output to be reported. The N function will raise the precision of the exact numbers to try and reach the requested precision. There are also concerns with, for example, how the replacing bigfloats are related to the original numbers. How important this depends on the use to which the calculation is being put and how sensitive the output is to changes in inputs - in the plotting above it is not important. -- Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 > Let's be realistic. If you want 60 digits of precision, too bad! -- in > the real world. There's nothing we can measure that closely. Drug > concentrations in clinical trials are generally measured within 15%, for > instance. Even machine precision is more than can be realistically > expected in any application I can think of. Even getting a satellite to > Jupiter probably involves more error in the final result than machine > precision. (If not, it's because we rely on ongoing corrections and > natural factors that put the satellite where it should be, such as > gravity drawing it toward each rendezvous -- not on that kind of > precision in propulsion or guidance.) > > So... unless all numerics in a problem have a theoretical origin, and > could be represented in Mathematica as Infinite precision expressions... > all this talk of higher-precision computation seems futile. > > The realistic question is this: given that I have confidence, say, in 6 > digits of precision for the inputs of an expression, how many digits can > I trust in the end result? Giving inputs MORE precision than they > deserve isn't the best way to answer that question. Here are two > methods of answering it in Mathematica. One uses bignums and the > other uses Intervals. > > Repetitious trial and error are NOT required either way. > > BIGNUMS: > > ClearAll[a, b, f, x, y] > f = x*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + y*b^8 + a/(2*b); > a = SetPrecision[77617, 6]; > b = SetPrecision[33096, 6]; > x = SetPrecision[33375/100, 6]; > y = SetPrecision[55/10, 6]; > InputForm[f] > > -4.1396`-12.5121*^19 > > Several previous solutions have set the precision or accuracy of f > before giving a and b (and possibly x and y) values. That results in > making the exponents imprecise along with all coefficients (not just x > and y), which may or may not be what we want. > > INTERVALS: > > I'll first digress to figure out what Interval is equivalent to 6-digit > precision. You might not actually do this if you like the Interval > method, but you have to decide SOMEHOW what Interval width to use. > > nums = {77617, 33096, 33375/100, 55/10}; > (Interval[SetPrecision[#1, 6]] & ) /@ nums /. Interval[{a_, b_}] :> > 2*((b - a)/(b + a)); > InputForm[%] > > {3.2209438653903`0.207*^-6, > 3.7768914672467`0.2761*^-6, > 2.9260299625467`0.1652*^-6, > 2.7743252840909`0.1421*^-6} > > For these numbers, # + Interval[{-1,1}]*#/630000& gets us very close, so > I'll use that. The second method therefore is: > > g = #1 + Interval[{-1, 1}]*(#1/630000) & ; > a = g[77617]; > b = g[33096]; > x = g[333.75]; > y = g[5.5]; > f = x*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + y*b^8 + a/(2*b) > (Min[f] + Max[f])/2 > > Interval[{-2.136361928005054*^32, 2.1363651195928296*^32}] > > 1.5957938878075177*^26 > > Either method shows the answer has no trustworthy digits, but I think > the second result is far easier to interpret. > > Here's another example, using the Sin function, whose derivative is Cos, > whose magnitude is bounded by one. The precision of the Sin of a number > should be GREATER than the precision of the number itself, especially > when Cos is small. > > a = SetPrecision[Pi/2, 6]; > InputForm[Sin[a]] > > 0.9999999999990905052982256654`11.6078 > > a = g[N[Pi/2]]; > Sin[a] - 1 > > Interval[{-3.1084024243455137*^-12, 0}] > > I hope this was worthwhile to someone. > > DrBob > > -----Original Message----- > > > > > On Friday, October 4, 2002, at 06:01 PM, DrBob > >[...] > > I would say this is correct and show that > SetPrecision is very useful > indeed. It tells you (what of course you ought > to already know in this > case anyway) that machine precision will not > give you a realiable > answer in this case. If you know your numbers > with a great deal of > accuracy you can get an accurate answer: > > In[24]:= > f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - > b^6 - > 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 100]; > a=SetPrecision[77617.,100]; b = > SetPrecision[33096.,100]; > > > In[26]:= > {f, Precision[f]} > > Out[26]= > > > {-0.82739605994682136814116509547981629199903311578438481991 > 781484167246798617832`61.2597, 61} > > > Congratulations! You just requested accuracy of > 100 for f and got 61 ( > to convince yourself add Accuracy[f] to In[26]). > If In[24] one > replaces SetAccuracy by SetPrecision the result is > similar. > > PK > [...] > > One has (initially) an accuracy of 100 for an > expression that contains > variables. > > In[25]:= Clear[a,b,f] > > In[26]:= f = SetAccuracy[333.75*b^6 + > a^2*(11*a^2*b^2 - b^6 - > 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 100]; > > In[27]:= Accuracy[f] > Out[27]= 100. > > Now we assign values to some indeterminants in f. > > In[28]:= a = SetPrecision[77617.,100]; b = > SetPrecision[33096.,100]; > > In[29]:= {f, Precision[f], Accuracy[f]} > Out[29]= > > {-0.8273960599468213681411650954798162919990331157843848199178148, > > 61.2599, 61.3422} > > The precision and accuracy has dropped. This is all > according to > standard numerical analysis regarding cancellation > error. You'll find it > in any textbook on the topic. > > > Assume that I want accuracy and precision of 100 for > f. You advice me to make experiments to find out, what > should be the initial precision and accuracy of a and > b to reach the requested accuracy and precision for f. > Notice, that you cannot just repeat I[26], we saw > already what happens. I have to re-type I[24], I[25], > I[26], I[27], I[28], and I[29] as many times as needed > to get f with accuracy and precision 100. > > Dan, you simply advocate to do MANUAL WORK that should > be done by machine. > > Let's suppose that in the above example I just want 60 > digits not 61. Precisely, I want 60 digits and nothing > or zeros afterwards. Let's see if I could use > SetAccuracy. > > In[30]:= > SetAccuracy[%, 60] > > Out[30]= > -0.82739605994682136814116509547981629199903311578438481991781 > > In[31]:= > % // FullForm > > Out[30]//FullForm= > -0.827396059946821368141165095479816291999033115784384819917814841672467 > 988` > 59.9177 > > Oops, it did not work (as expected). Let's highlight > with mouse the expression in Out[30] and copy to a new > cell. Oops, we got > -0.827396059946821368141165095479816291999033115784384819917814841672467 > 988` > 59.9177 > again. Let's change Out[30] to a text cell and then > copy. > > In[31]:= > -0.82739605994682136814116509547981629199903311578438481991781 > > Out[31]= > -0.82739605994682136814116509547981629199903311578438481991781 > > Success? Not so fast. > > In[32]:= > % // FullForm > > Out[32]//FullForm= > -0.827396059946821368141165095479816291999033115784384819917809999999999 > 998635 > 08`59.2041 > > Dan, is there any simple way to get what I want? > > As I repeated already number of times, at this stage > of the development of computer technology, software > should do it for me (!). We both know that this is > doable. Some of the textbooks that you just advised me > to read describe it. As a developer of Mathematica, > tell us why do you consider this to be a bad idea? > > Peter Kosta > > As for what happens when you artificially raise > precision (or accuracy) > of machine numbers far beyond that guaranteed by > their internal > representation, that falls into to category of > garbage in, garbage out. > It is, howoever, valid to use SetPrecision to raise > precision in > (typically iterative) algorithms where significance > arithmetic might be > unduly pessimistic due to incorrect assumptions > about uncorollatedness > of numerical error. Examples of such usage have > appeared in this news > group. > > > Daniel Lichtblau > Wolfram Research > > > __________________________________________________ > Do you Yahoo!? > http://faith.yahoo.com > > > > Reply-To: ==== Here's a start, with a more complicated function: lst = Range[50]; f = #1^2 & ; k = 7; r = Range[k + 1, Length[lst], k]; lst[[r]] = lst[[r]] + f /@ lst[[r - 1]]; lst {1, 2, 3, 4, 5, 6, 7, 57, 9, 10, 11, 12, 13, 14, 211, 16, 17, 18, 19, 20, 21, 463, 23, 24, 25, 26, 27, 28, 813, 30, 31, 32, 33, 34, 35, 1261, 37, 38, 39, 40, 41, 42, 1807, 44, 45, 46, 47, 48, 49, 2451} or: g[lst_List, k_Integer?Positive, f_Function] := Module[ {result = lst, r = Range[k + 1, Length[lst], k]}, result[[r]] = result[[r]] + f /@ lst[[r - 1]]; result ] lst = Range[50]; lst = g[lst, 7, #1^2 & ] or: lst = Range[50]; lst + Drop[Prepend[MapIndexed[If[Mod[First@#2, k] == 0, f@#1, 0] &, lst], 0], -1] DrBob -----Original Message----- advantage of native List operators, but still not be too recondite? I've been thinking about multiplying a copy of myList by a mask list {0,0,1,0,0,1,..} to generate a masked copy and approaches like that. Better ways??? ==== I'm playing around with Mathematica just trying to see what happens if... One thing I came up with was lst1={a,b,c}; lst2={d,e,f}; Map[lst1,lst2} which resulted in the following rather unusual looking expression(?): {{a, b, c}[d], {a, b, c}[e], {a, b, c}[f]} I'm wondering if such a List represents something 'meaningful'. Any opinions? STH . ==== > I'm playing around with Mathematica just trying to see what happens if... One thing > I came up with was lst1={a,b,c}; lst2={d,e,f}; Map[lst1,lst2} which > resulted in the following rather unusual looking expression(?): > > {{a, b, c}[d], {a, b, c}[e], {a, b, c}[f]} > > I'm wondering if such a List represents something 'meaningful'. Any > opinions? Through /@ {{a, b, c}[d], {a, b, c}[e], {a, b, c}[f]} {{a[d], b[d], c[d]}, {a[e], b[e], c[e]}, {a[f], b[f], c[f]}} which might be useful. -- Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 ==== >I want to apply a function to every k-th element of a long list and >add the result to the k+1 element. > >[Actually k = 3 and I just want to multiply myList[[k]] by a >constant (independent of k) and add the result to myList[[k+1]] for >every value of k that's divisible by 3.] lst = Table[Random[], {20}] fact = 2 Here's a matrix method. f represents an indexed element of a matrix. f[x_,x_] := If[Mod[x,3]==0, fact, 1] f[x_, y_] := If[Mod[x,3]==0&&y[Equal]x+1, 2, 0] Create a matrix from the elements. arr = Array[f,{Length[lst], Length[lst]}]; Then matrix multiply. newlst1 = lst.arr Here's another way with the highly underrated MapIndexed. Create pairs of the nth and n-1th values. pairs = Transpose[{lst, Prepend[Drop[lst,-1],0]}] Create a function that takes the nth value (val), the n-1th value (prevval), and the index (num). multlst[{val_, prevval_}, {num_}] := Switch[Mod[num, 3], 0, fact*val, 1,val+ fact*prevval, 2, val ] Then, MapIndexed across the pairs. newlst2=MapIndexed[multlst, pairs] -------------------------------------------------------------- Omega Consulting The final answer to your Mathematica needs Spend less time searching and more time finding. http://www.wz.com/internet/Mathematica.html ==== I have two equations that I have solved for: x[n_] := 2331 + 8 n y[n_] := -3108 - 11n I want to include only solutions which are non-negative, that is x >= 0 and y >= 0. In this example we can do 2331 + 8n > = 0 and solve for n, n >= -291.375 and -3108 - 11 n >= 0 and solve for n, n <= -282.545 So we have -291.375 <= n <= -282.545. The integer solution set here is for n = {-290, -289, -288, -287, -286, -285, -284, -283}. So in this case we have 8 non-negative solutions. Given that I can supply x[n] and y[n], how do I go about finding the set n? ==== Here is a way of doing what you want (output cells are indented): <=0,y[n]>=0},n] -(2331/8) <= n <= -(3108/11) FullForm[soln] Inequality[Rational[-2331,8],LessEqual,n,LessEqual,Rational[-3108,11]] Range[Apply[Sequence,{Ceiling[soln[[1]]],Floor[soln[[-1]]]}]] {-291, -290, -289, -288, -287, -286, -285, -284, -283} -- Steve Luttrell West Malvern, UK > > I have two equations that I have solved for: > > x[n_] := 2331 + 8 n > y[n_] := -3108 - 11n > > I want to include only solutions which are non-negative, that is x >= 0 and > y >= 0. > > In this example we can do 2331 + 8n > = 0 and solve for n, n >= -291.375 > and -3108 - 11 n >= 0 and solve for n, n <= -282.545 > > So we have -291.375 <= n <= -282.545. > > The integer solution set here is for n = > {-290, -289, -288, -287, -286, -285, -284, -283}. > > So in this case we have 8 non-negative solutions. > > Given that I can supply x[n] and y[n], how do I go about finding the set n? > > > > ==== Let's first consider your original problem and take a small list as an example: mylist = {a, b, c, d, e, f, g}; As you pointed out there is a rather obvious and natural way to do it using the Do loop. In[5]:= Do[mylist[[3*i + 1]] = k*mylist[[3*i]] + mylist[[3*i + 1]], {i, 1, Length[mylist]/3}]; mylist Out[5]= {a, b, c, d + c*k, e, f, g + f*k} One can also do it using something like your second approach. Notice the need to re-set mylist which got changed by the Do loop: In[6]:= mylist={a,b,c,d,e,f,g}; In[7]:= RotateRight[Table[If[Mod[i,3]==0,k, 0],{i,1,Length[mylist]}]*mylist]+mylist Out[7]= {a,b,c,d+c k,e,f,g+f k} There is an ambiguity that appears if the length of the list is exactly divisible by 3. In that case k times the last element should be added to the next element. In this case the first approach (using Do) will produce an error message while the second will interpret the next element to mean the first element of the list. One can fix that but I shan't bother to do so and assume that the length of the list is not divisible by 3. Now let's take a large list and compare the performance of the two approaches: In[8]:= k = 2; In[9]:= mylist = Table[Random[Integer, {1, 9}], {10000}]; In[11]:= mylist1 = mylist; In[12]:= Timing[list1 = (Do[mylist[[3*i + 1]] = k*mylist[[3*i]] + mylist[[3*i + 1]], {i, 1, Length[mylist]/3}]; mylist); ] Out[12]= {0.21999999999999997*Second, Null} In[13]:= mylist = mylist1; In[14]:= Timing[list2 = RotateRight[Table[If[Mod[i, 3] == 0, k, 0], {i, 1, Length[mylist]}]*mylist] + mylist; ] Out[14]= {0.020000000000000018*Second, Null} In[15]:= list1 == list2 Out[15]= True So you were right, at least in this implementation the second approach turns out to be much faster. Note also however, that when you change the problem to the one that you originally stated the two approaches will no longer give the same answer. The reason is that your statement is ambiguous: apply a function to every k-th element of a long list and add the result to the k+1 element can either mean that you want to take the original k-th element of the original list, multiply by a constant and add to the next one, or do the same with the already altered k-th element (by the previous step of the procedure). The Do loop approach will do the latter and the other the former. Andrzej Kozlowski Yokohama, Japan http://www.mimuw.edu.pl/~akoz/ http://platon.c.u-tokyo.ac.jp/andrzej/ > I want to apply a function to every k-th element of a long list and > add the result to the k+1 element. > > [Actually k = 3 and I just want to multiply myList[[k]] by a > constant (independent of k) and add the result to myList[[k+1]] for > every value of k that's divisible by 3.] > > Is there a way to do this -- or in general to get at every k-th > element of a list -- that's faster and more elegant than writing a > brute > force Do[] loop or using Mod[] operators, and that will take > advantage of native List operators, but still not be too recondite? > > I've been thinking about multiplying a copy of myList by a mask > list > {0,0,1,0,0,1,..} to generate a masked copy and approaches like that. > Better ways??? > > > ==== Can Mathematica factor the polynomial p1=x^6+9/14*x^5+9/28*x^4+3/35*x^3+9/700*x^2+9/8750*x+3/87500; without a priori knowledge of the Extension field? ==== Carlos, Futher to my previous posting (which gave the code for the function FactorR used below), here is a complete factorisation by radicals. I also test that the product of the factors gives the original polynomial. We want to factor the polynomial p1 = x^6 + (9/14)*x^5 + (9/28)*x^4 + (3/35)*x^3 + (9/700)*x^2 + (9/8750)*x + 3/87500; in radicals. We can't expect this to be easy or even possible in terms of radicals (the general quintic is not solvable interms of radicals). But, using the function FactorR given in my posting, Re:factoring quartic over radicals, sent a few days ago (08/012/02) , we get p2 = FactorR[p1, x] (x^2 - 2*x*Root[3 + 45*#1 + 225*#1^2 + 700*#1^3 & , 1] + Root[-3 + 225*#1^2 - 5625*#1^4 + 87500*#1^6 & , 2]^2)* (x^2 - 2*x*Root[1827 + 65340*#1 + 974700*#1^2 + 7824000*#1^3 + 36360000*#1^4 + 100800000*#1^5 + 156800000*#1^6 & , 2] + Root[9 - 1350*#1^2 + 84375*#1^4 - 3056250*#1^6 - 11250000*#1^8 - 984375000*#1^10 + 7656250000*#1^12 & , 3]^2)* (x^2 - 2*x*Root[1827 + 65340*#1 + 974700*#1^2 + 7824000*#1^3 + 36360000*#1^4 + 100800000*#1^5 + 156800000*#1^6 & , 1] + Root[9 - 1350*#1^2 + 84375*#1^4 - 3056250*#1^6 - 11250000*#1^8 - 984375000*#1^10 + 7656250000*#1^12 & , 4]^2) Try to change the root objects to radical form: p3 = p2 /. r_Root :> ToRadicals[r] (3/140 + (1/140)*(13/5)^(2/3)*3^(1/3) - (1/140)*(13/5)^(1/3)*3^(2/3) - 2*(-(3/28) - (1/28)*(13/5)^(2/3)*3^(1/3) + (1/28)*(13/5)^(1/3)*3^(2/3))*x + x^2)*(x^2 - 2*x*Root[1827 + 65340*#1 + 974700*#1^2 + 7824000*#1^3 + 36360000*#1^4 + 100800000*#1^5 + 156800000*#1^6 & , 2] + Root[9 - 1350*#1 + 84375*#1^2 - 3056250*#1^3 - 11250000*#1^4 - 984375000*#1^5 + 7656250000*#1^6 & , 1])* (x^2 - 2*x*Root[1827 + 65340*#1 + 974700*#1^2 + 7824000*#1^3 + 36360000*#1^4 + 100800000*#1^5 + 156800000*#1^6 & , 1] + Root[9 - 1350*#1 + 84375*#1^2 - 3056250*#1^3 - 11250000*#1^4 - 984375000*#1^5 + 7656250000*#1^6 & , 2]) We succeeded with the first factor: f1 = p3[[1]] 3/140 + (1/140)*(13/5)^(2/3)*3^(1/3) - (1/140)*(13/5)^(1/3)*3^(2/3) - 2*(-(3/28) - (1/28)*(13/5)^(2/3)*3^(1/3) + (1/28)*(13/5)^(1/3)*3^(2/3))*x + x^2 The product of the other two factors, in a form avoiding root objects, is easily found by division: q = PolynomialQuotient[p1, f1, x] 3/3500 + ((13/5)^(1/3)*3^(2/3))/3500 + (3/175 + (1/175)*(13/5)^(1/3)*3^(2/3))*x + (9/70 - (1/140)*(13/5)^(2/3)*3^(1/3) + (1/28)*(13/5)^(1/3)*3^(2/3))*x^2 + (3/7 - (1/14)*(13/5)^(2/3)*3^(1/3) + (1/14)*(13/5)^(1/3)*3^(2/3))*x^3 + x^4 Try FactorR on this f23 = FactorR[q, x] (x^2 + ((-(1/2))*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)*3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920] + Im[(1/2)*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)*3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920 + (-2250 + 25*13^(2/3)*15^(1/3) - 125*13^(1/3)*15^(2/3))/8750 + (3*(7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2)/ 1225000000 - (-((2*(300 + 20*13^(1/3)*15^(2/3)))/4375) + (1/17500)*((7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))* ((2250 - 25*13^(2/3)*15^(1/3) + 125*13^(1/3)*15^(2/3))/4375 - (7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2/ 306250000)))/(4*Sqrt[-(117/1960) - (9/784)*(13/5)^(2/3)*3^(1/3) - (39*(13/5)^(1/3)*3^(2/3))/3920])]])^2 - 2*x*((-7500 + 250*13^(2/3)*15^(1/3) - 250*13^(1/3)*15^(2/3))/70000 + Re[(1/2)*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)*3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920 + (-2250 + 25*13^(2/3)*15^(1/3) - 125*13^(1/3)*15^(2/3))/8750 + (3*(7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2)/ 1225000000 - (-((2*(300 + 20*13^(1/3)*15^(2/3)))/4375) + (1/17500)*((7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))* ((2250 - 25*13^(2/3)*15^(1/3) + 125*13^(1/3)*15^(2/3))/4375 - (7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2/ 306250000)))/(4*Sqrt[-(117/1960) - (9/784)*(13/5)^(2/3)*3^(1/3) - (39*(13/5)^(1/3)*3^(2/3))/3920])]]) + ((-7500 + 250*13^(2/3)*15^(1/3) - 250*13^(1/3)*15^(2/3))/70000 + Re[(1/2)*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)*3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920 + (-2250 + 25*13^(2/3)*15^(1/3) - 125*13^(1/3)*15^(2/3))/8750 + (3*(7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2)/ 1225000000 - (-((2*(300 + 20*13^(1/3)*15^(2/3)))/4375) + (1/17500)*((7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))* ((2250 - 25*13^(2/3)*15^(1/3) + 125*13^(1/3)*15^(2/3))/4375 - (7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2/ 306250000)))/(4*Sqrt[-(117/1960) - (9/784)*(13/5)^(2/3)*3^(1/3) - (39*(13/5)^(1/3)*3^(2/3))/3920])]])^2)* (x^2 + ((1/2)*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)*3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920] + Im[(-(1/2))*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)*3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920 + (-2250 + 25*13^(2/3)*15^(1/3) - 125*13^(1/3)*15^(2/3))/8750 + (3*(7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2)/ 1225000000 + (-((2*(300 + 20*13^(1/3)*15^(2/3)))/4375) + (1/17500)*((7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))* ((2250 - 25*13^(2/3)*15^(1/3) + 125*13^(1/3)*15^(2/3))/4375 - (7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2/ 306250000)))/(4*Sqrt[-(117/1960) - (9/784)*(13/5)^(2/3)*3^(1/3) - (39*(13/5)^(1/3)*3^(2/3))/3920])]])^2 - 2*x*((-7500 + 250*13^(2/3)*15^(1/3) - 250*13^(1/3)*15^(2/3))/70000 + Re[(-(1/2))*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)*3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920 + (-2250 + 25*13^(2/3)*15^(1/3) - 125*13^(1/3)*15^(2/3))/8750 + (3*(7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2)/ 1225000000 + (-((2*(300 + 20*13^(1/3)*15^(2/3)))/4375) + (1/17500)*((7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))* ((2250 - 25*13^(2/3)*15^(1/3) + 125*13^(1/3)*15^(2/3))/4375 - (7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2/ 306250000)))/(4*Sqrt[-(117/1960) - (9/784)*(13/5)^(2/3)*3^(1/3) - (39*(13/5)^(1/3)*3^(2/3))/3920])]]) + ((-7500 + 250*13^(2/3)*15^(1/3) - 250*13^(1/3)*15^(2/3))/70000 + Re[(-(1/2))*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)*3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920 + (-2250 + 25*13^(2/3)*15^(1/3) - 125*13^(1/3)*15^(2/3))/8750 + (3*(7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2)/ 1225000000 + (-((2*(300 + 20*13^(1/3)*15^(2/3)))/4375) + (1/17500)*((7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))* ((2250 - 25*13^(2/3)*15^(1/3) + 125*13^(1/3)*15^(2/3))/4375 - (7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2/ 306250000)))/(4*Sqrt[-(117/1960) - (9/784)*(13/5)^(2/3)*3^(1/3) - (39*(13/5)^(1/3)*3^(2/3))/3920])]])^2) We try to get rid of the parts Re[.] and Im[.]:, f231 = f23 /. z:(_Re | _Im) :> ToRadicals[FullSimplify[z]] (((-7500 + 250*13^(2/3)*15^(1/3) - 250*13^(1/3)*15^(2/3))/70000 - (1/280)*Sqrt[3*(-390 + 13*13^(2/3)*15^(1/3) + 15*13^(1/3)*15^(2/3))])^2 + ((1/2)*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)*3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920] + (1/280)*Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)])^2 - 2*((-7500 + 250*13^(2/3)*15^(1/3) - 250*13^(1/3)*15^(2/3))/70000 - (1/280)*Sqrt[3*(-390 + 13*13^(2/3)*15^(1/3) + 15*13^(1/3)*15^(2/3))])*x + x^2)*(((-7500 + 250*13^(2/3)*15^(1/3) - 250*13^(1/3)*15^(2/3))/70000 + (1/280)*Sqrt[3*(-390 + 13*13^(2/3)*15^(1/3) + 15*13^(1/3)*15^(2/3))])^2 + ((-(1/2))*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)*3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920] + (1/280)*Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)])^2 - 2*((-7500 + 250*13^(2/3)*15^(1/3) - 250*13^(1/3)*15^(2/3))/70000 + (1/280)*Sqrt[3*(-390 + 13*13^(2/3)*15^(1/3) + 15*13^(1/3)*15^(2/3))])*x + x^2) We now have the ramaining two factors in radical form, but a little simplification helps: f232 = f231 /. (n_)?NumericQ :> Simplify[n] ((1/78400)*(30 - 13^(2/3)*15^(1/3) + 13^(1/3)*15^(2/3) + Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)])^2 + (1/78400)*(Sqrt[1170 + 45*13^(2/3)*15^(1/3) + 39*13^(1/3)*15^(2/3)] + Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)])^2 - (1/140)*(-30 + 13^(2/3)*15^(1/3) - 13^(1/3)*15^(2/3) - Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)])*x + x^2)* ((1/78400)*(-30 + 13^(2/3)*15^(1/3) - 13^(1/3)*15^(2/3) + Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)])^2 + (1/78400)*(Sqrt[1170 + 45*13^(2/3)*15^(1/3) + 39*13^(1/3)*15^(2/3)] - Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)])^2 - (1/140)*(-30 + 13^(2/3)*15^(1/3) - 13^(1/3)*15^(2/3) + Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)])*x + x^2) TEST Test if the product of the factors is equal to p1: prd1 = Collect[Expand[f232*f1], x] 172077/3841600000 - (4959*(13/5)^(2/3)*3^(1/3))/768320000 + (117*(13/5)^(1/3)*3^(2/3))/27440000 - (1/1920800000)* (3*Sqrt[1170 + 45*13^(2/3)*15^(1/3) + 39*13^(1/3)*15^(2/3)]* Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)]* Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)]) + (1/3841600000)*((13/5)^(1/3)*3^(2/3)*Sqrt[1170 + 45*13^(2/3)*15^(1/3) + 39*13^(1/3)*15^(2/3)]*Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)]*Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)]) + (491193/384160000 - (7731*(13/5)^(2/3)*3^(1/3))/76832000 + (117*(13/5)^(1/3)*3^(2/3))/2744000 - (1/384160000)* (9*Sqrt[1170 + 45*13^(2/3)*15^(1/3) + 39*13^(1/3)*15^(2/3)]* Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)]* Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)]) - (1/384160000)*((13/5)^(2/3)*3^(1/3)*Sqrt[1170 + 45*13^(2/3)*15^(1/3) + 39*13^(1/3)*15^(2/3)]*Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)]*Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)]) + (1/192080000)*((13/5)^(1/3)*3^(2/3)* Sqrt[1170 + 45*13^(2/3)*15^(1/3) + 39*13^(1/3)*15^(2/3)]* Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)]* Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)]))*x + (1087173/76832000 - (12771*(13/5)^(2/3)*3^(1/3))/30732800 + (5967*(13/5)^(1/3)*3^(2/3))/30732800 - (1/76832000)* (9*Sqrt[1170 + 45*13^(2/3)*15^(1/3) + 39*13^(1/3)*15^(2/3)]* Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)]* Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)]) - (1/153664000)*(3*(13/5)^(2/3)*3^(1/3)*Sqrt[1170 + 45*13^(2/3)*15^(1/3) + 39*13^(1/3)*15^(2/3)]*Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)]*Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)]) + (1/153664000)*(3*(13/5)^(1/3)*3^(2/3)* Sqrt[1170 + 45*13^(2/3)*15^(1/3) + 39*13^(1/3)*15^(2/3)]* Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)]* Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)]))*x^2 + (23871/274400 - (99*(13/5)^(2/3)*3^(1/3))/109760 + (663*(13/5)^(1/3)*3^(2/3))/ 548800 - (1/2744000)*(Sqrt[1170 + 45*13^(2/3)*15^(1/3) + 39*13^(1/3)*15^(2/3)]*Sqrt[-1170 + 39*13^(2/3)*15^(1/3) + 45*13^(1/3)*15^(2/3)]*Sqrt[1170 + 73*13^(2/3)*15^(1/3) + 67*13^(1/3)*15^(2/3)]))*x^3 + (9*x^4)/28 + (9*x^5)/14 + x^6 prd1 /. (n_)?NumericQ :> ToRadicals[FullSimplify[n]] 172077/3841600000 - (4959*(13/5)^(2/3)*3^(1/3))/768320000 + (117*(13/5)^(1/3)*3^(2/3))/27440000 - (9*(234 + 221*(13/5)^(1/3)*3^(2/3) - 33*13^(2/3)*15^(1/3)))/384160000 + (117*(-165 + 17*13^(2/3)*15^(1/3) + 6*13^(1/3)*15^(2/3)))/3841600000 + (9*x)/8750 + (9*x^2)/700 + (3*x^3)/35 + (9*x^4)/28 + (9*x^5)/14 + x^6 Together[%] (3 + 90*x + 1125*x^2 + 7500*x^3 + 28125*x^4 + 56250*x^5 + 87500*x^6)/87500 Apart[%] 3/87500 + (9*x)/8750 + (9*x^2)/700 + (3*x^3)/35 + (9*x^4)/28 + (9*x^5)/14 + x^6 This is p1: p1 3/87500 + (9*x)/8750 + (9*x^2)/700 + (3*x^3)/35 + (9*x^4)/28 + (9*x^5)/14 + x^6 ------------------ It is ususlly better to try to reduce a difference to zero than to reduce one form to another tst1 = Collect[Expand[f232*f1 - p1], x] tst2 = tst1 /. (n_)?NumericQ :> ToRadicals[FullSimplify[n]] 8073/768320000 - (4959*(13/5)^(2/3)*3^(1/3))/768320000 + (117*(13/5)^(1/3)*3^(2/3))/27440000 - (9*(234 + 221*(13/5)^(1/3)*3^(2/3) - 33*13^(2/3)*15^(1/3)))/384160000 + (117*(-165 + 17*13^(2/3)*15^(1/3) + 6*13^(1/3)*15^(2/3)))/3841600000 Together[%] 0 -- Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 > Can Mathematica factor the polynomial > > p1=x^6+9/14*x^5+9/28*x^4+3/35*x^3+9/700*x^2+9/8750*x+3/87500; > > without a priori knowledge of the Extension field? > ==== Carlos, p1=x^6+9/14*x^5+9/28*x^4+3/35*x^3+9/700*x^2+9/8750*x+3/87500; We can't expect this to be easy or even possible in terms of radicals (the general quintic is not solvable interms of radicals). But, using the function FactorR given in my posting, Re:factoring quartic over radicals, sent a few days (08/012/02) ago, we get p2=FactorR[p1,x] (x^2 - 2*x*Root[3 + 45*#1 + 225*#1^2 + 700*#1^3 & , 1] + Root[-3 + 225*#1^2 - 5625*#1^4 + 87500*#1^6 & , 2]^2)* (x^2 - 2*x*Root[1827 + 65340*#1 + 974700*#1^2 + 7824000*#1^3 + 36360000*#1^4 + 100800000*#1^5 + 156800000*#1^6 & , 2] + Root[9 - 1350*#1^2 + 84375*#1^4 - 3056250*#1^6 - 11250000*#1^8 - 984375000*#1^10 + 7656250000*#1^12 & , 3]^2)* (x^2 - 2*x*Root[1827 + 65340*#1 + 974700*#1^2 + 7824000*#1^3 + 36360000*#1^4 + 100800000*#1^5 + 156800000*#1^6 & , 1] + Root[9 - 1350*#1^2 + 84375*#1^4 - 3056250*#1^6 - 11250000*#1^8 - 984375000*#1^10 + 7656250000*#1^12 & , 4]^2) Try to express the factors in terms of radicals p3=ToRadicals/@p2 (3/140 + (1/140)*(13/5)^(2/3)*3^(1/3) - (1/140)*(13/5)^(1/3)*3^(2/3) - 2*(-(3/28) - (1/28)*(13/5)^(2/3)*3^(1/3) + (1/28)*(13/5)^(1/3)*3^(2/3))*x + x^2)* (x^2 - 2*x*Root[1827 + 65340*#1 + 974700*#1^2 + 7824000*#1^3 + 36360000*#1^4 + 100800000*#1^5 + 156800000*#1^6 & , 2] + Root[9 - 1350*#1 + 84375*#1^2 - 3056250*#1^3 - 11250000*#1^4 - 984375000*#1^5 + 7656250000*#1^6 & , 1])* (x^2 - 2*x*Root[1827 + 65340*#1 + 974700*#1^2 + 7824000*#1^3 + 36360000*#1^4 + 100800000*#1^5 + 156800000*#1^6 & , 1] + Root[9 - 1350*#1 + 84375*#1^2 - 3056250*#1^3 - 11250000*#1^4 - 984375000*#1^5 + 7656250000*#1^6 & , 2]) We succeeded with the first factor f1=p3[[1]] 3/140 + (1/140)*(13/5)^(2/3)*3^(1/3) - (1/140)*(13/5)^(1/3)*3^(2/3) - 2*(-(3/28) - (1/28)*(13/5)^(2/3)*3^(1/3) + (1/28)*(13/5)^(1/3)*3^(2/3))*x + x^2 The product of the other two factors, in a form avoiding root objects, is easily found: q=PolynomialQuotient[p1,f1,x] 3/3500 + ((13/5)^(1/3)*3^(2/3))/3500 + (3/175 + (1/175)*(13/5)^(1/3)*3^(2/3))*x + (9/70 - (1/140)*(13/5)^(2/3)*3^(1/3) + (1/28)*(13/5)^(1/3)*3^(2/3))*x^2 + (3/7 - (1/14)*(13/5)^(2/3)*3^(1/3) + (1/14)*(13/5)^(1/3)*3^(2/3))*x^3 + x^4 5*(3 + 10*x*(6 + 5*x*(9 + 10*x*(3 + 7*x))))) Try FactorR on this f23=FactorR[q,x] (x^2 + ((-(1/2))*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)* 3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920] + Im[(1/2)*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)* 3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920 + (-2250 + 25*13^(2/3)*15^(1/3) - 125*13^(1/3)* 15^(2/3))/8750 + (3*(7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)* 15^(2/3))^2)/1225000000 - (-((2*(300 + 20*13^(1/3)*15^(2/3)))/4375) + (1/17500)*((7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))*((2250 - 25*13^(2/3)* 15^(1/3) + 125*13^(1/3)*15^(2/3))/4375 - (7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2/306250000)))/ (4*Sqrt[-(117/1960) - (9/784)*(13/5)^(2/3)* 3^(1/3) - (39*(13/5)^(1/3)*3^(2/3))/ 3920])]])^2 - 2*x*((-7500 + 250*13^(2/3)*15^(1/3) - 250*13^(1/3)*15^(2/3))/70000 + Re[(1/2)*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)* 3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920 + (-2250 + 25*13^(2/3)*15^(1/3) - 125*13^(1/3)* 15^(2/3))/8750 + (3*(7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)* 15^(2/3))^2)/1225000000 - (-((2*(300 + 20*13^(1/3)*15^(2/3)))/4375) + (1/17500)*((7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))*((2250 - 25*13^(2/3)* 15^(1/3) + 125*13^(1/3)*15^(2/3))/4375 - (7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2/306250000)))/ (4*Sqrt[-(117/1960) - (9/784)*(13/5)^(2/3)* 3^(1/3) - (39*(13/5)^(1/3)*3^(2/3))/ 3920])]]) + ((-7500 + 250*13^(2/3)*15^(1/3) - 250*13^(1/3)* 15^(2/3))/70000 + Re[(1/2)*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)* 3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920 + (-2250 + 25*13^(2/3)*15^(1/3) - 125*13^(1/3)* 15^(2/3))/8750 + (3*(7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)* 15^(2/3))^2)/1225000000 - (-((2*(300 + 20*13^(1/3)*15^(2/3)))/4375) + (1/17500)*((7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))*((2250 - 25*13^(2/3)* 15^(1/3) + 125*13^(1/3)*15^(2/3))/4375 - (7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2/306250000)))/ (4*Sqrt[-(117/1960) - (9/784)*(13/5)^(2/3)* 3^(1/3) - (39*(13/5)^(1/3)*3^(2/3))/ 3920])]])^2)* (x^2 + ((1/2)*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)* 3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920] + Im[(-(1/2))*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)* 3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920 + (-2250 + 25*13^(2/3)*15^(1/3) - 125*13^(1/3)* 15^(2/3))/8750 + (3*(7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)* 15^(2/3))^2)/1225000000 + (-((2*(300 + 20*13^(1/3)*15^(2/3)))/4375) + (1/17500)*((7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))*((2250 - 25*13^(2/3)* 15^(1/3) + 125*13^(1/3)*15^(2/3))/4375 - (7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2/306250000)))/ (4*Sqrt[-(117/1960) - (9/784)*(13/5)^(2/3)* 3^(1/3) - (39*(13/5)^(1/3)*3^(2/3))/ 3920])]])^2 - 2*x*((-7500 + 250*13^(2/3)*15^(1/3) - 250*13^(1/3)*15^(2/3))/70000 + Re[(-(1/2))*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)* 3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920 + (-2250 + 25*13^(2/3)*15^(1/3) - 125*13^(1/3)* 15^(2/3))/8750 + (3*(7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)* 15^(2/3))^2)/1225000000 + (-((2*(300 + 20*13^(1/3)*15^(2/3)))/4375) + (1/17500)*((7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))*((2250 - 25*13^(2/3)* 15^(1/3) + 125*13^(1/3)*15^(2/3))/4375 - (7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2/306250000)))/ (4*Sqrt[-(117/1960) - (9/784)*(13/5)^(2/3)* 3^(1/3) - (39*(13/5)^(1/3)*3^(2/3))/ 3920])]]) + ((-7500 + 250*13^(2/3)*15^(1/3) - 250*13^(1/3)* 15^(2/3))/70000 + Re[(-(1/2))*Sqrt[117/1960 + (9/784)*(13/5)^(2/3)* 3^(1/3) + (39*(13/5)^(1/3)*3^(2/3))/3920 + (-2250 + 25*13^(2/3)*15^(1/3) - 125*13^(1/3)* 15^(2/3))/8750 + (3*(7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)* 15^(2/3))^2)/1225000000 + (-((2*(300 + 20*13^(1/3)*15^(2/3)))/4375) + (1/17500)*((7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))*((2250 - 25*13^(2/3)* 15^(1/3) + 125*13^(1/3)*15^(2/3))/4375 - (7500 - 250*13^(2/3)*15^(1/3) + 250*13^(1/3)*15^(2/3))^2/306250000)))/ (4*Sqrt[-(117/1960) - (9/784)*(13/5)^(2/3)* 3^(1/3) - (39*(13/5)^(1/3)*3^(2/3))/ 3920])]])^2) We still use Re and Im -- Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 > Can Mathematica factor the polynomial > > p1=x^6+9/14*x^5+9/28*x^4+3/35*x^3+9/700*x^2+9/8750*x+3/87500; > > without a priori knowledge of the Extension field? > ==== Dear MathGroup Experts, Is there a slick way to use Mathematica to plot the phase portrait of a pair of (straightforward) differential equations? I'd also like to see a plot of the associated vector/direction field. I found a good package called DynPac, but it's very cumbersome to use, and I was hoping there's another tool out there that my students would feel more comfortable using. Jason -- Jason Miller, Ph.D. Division of Mathematics and Computer Science Truman State University 100 East Normal St. Kirksville, MO 63501 http://vh216801.truman.edu 660.785.7430 ==== >I have two equations that I have solved for: > >x[n_] := 2331 + 8 n >y[n_] := -3108 - 11n > >I want to include only solutions which are non-negative, that is x >= 0 >and >y >= 0. > >In this example we can do 2331 + 8n > = 0 and solve for n, n >= -291.375 >and -3108 - 11 n >= 0 and solve for n, n <= -282.545 > >So we have -291.375 <= n <= -282.545. > >The integer solution set here is for n = >{-290, -289, -288, -287, -286, -285, -284, -283}. > >So in this case we have 8 non-negative solutions. > >Given that I can supply x[n] and y[n], how do I go about finding the set >n? > Needs[Algebra`InequalitySolve`]; x[n_] := 2331 + 8n; y[n_] := -3108 - 11n; rng = InequalitySolve[{x[n] >= 0, y[n] >= 0}, n]; soln = Range[Ceiling[rng[[1]]], rng[[-1]]] {-291, -290, -289, -288, -287, -286, -285, -284, -283} Length[soln] 9 Bob Hanlon ==== Factoring without specifying the extension does not really make sense. Of course Mathematica can easily factor yur polynomial into linear factors over the complex numbers (with the help of Solve), but I suspect you are really asking for is factoring over the reals. This is harder and needs more human input. But anyway, Mathematica can do this, or at least I have done it using Mathematica. In fact if you are satisfied with a numerical answer Mathematica can do alone and in seconds: In[1]:= Simplify[N[x^6 + (9/14)*x^5 + (9/28)*x^4 + (3/35)*x^3 + (9/700)*x^2 + (9/8750)*x + 3/87500]] Out[1]= 1.*(0.010974992601737198 + 0.20255610310498295*x + x^2)*(0.020476912388332692 + 0.2047691238833268*x + x^2)*(0.15256133957420948 + 0.23553191586883315*x + x^2) But I have in fact been foolish enough to compute the exact answer too. I do not propose to post it here for it's absolutely horrible (expressed in terms of Root objects) and quite useless. However if you really want to see it I can send it to you privately. Andrzej Kozlowski Yokohama, Japan http://www.mimuw.edu.pl/~akoz/ http://platon.c.u-tokyo.ac.jp/andrzej/ > To: mathgroup@smc.vnet.net > > Can Mathematica factor the polynomial > > p1=x^6+9/14*x^5+9/28*x^4+3/35*x^3+9/700*x^2+9/8750*x+3/87500; > > without a priori knowledge of the Extension field? > > > ==== > > Can Mathematica factor the polynomial > > p1=x^6+9/14*x^5+9/28*x^4+3/35*x^3+9/700*x^2+9/8750*x+3/87500; > > without a priori knowledge of the Extension field? The endless sci.math.symbolic thread strikes MathGroup! Not exactly possible with no prior knowledge. Factor must work with a given field, and the default is the rationals. You might direct it, say by using the discriminant of the polynomial (as pointed out by Peter Montgomery and Stephen Forrest on sci.math.symbolic. One may do this in Mathematica as: p1 = x^6+9/14*x^5+9/28*x^4+3/35*x^3+9/700*x^2+9/8750*x+3/87500; I cribbed code for Discriminant right from www.mathworld.com: Discriminant[p_?PolynomialQ,x_] := With[{n = Exponent[p,x]}, Cancel[((-1)^(n(n-1)/2)Resultant[p,D[p,x],x])/Coefficient[p,x,n]^(2n-1)]] In[3]:= InputForm[Factor[p1, Extension->Sqrt[Discriminant[p1,x]]]] Out[3]//InputForm= ((-15*I + Sqrt[195] - (225*I)*x + 15*Sqrt[195]*x - (1125*I)*x^2 + 75*Sqrt[195]*x^2 - (3500*I)*x^3)*(15*I + Sqrt[195] + (225*I)*x + 15*Sqrt[195]*x + (1125*I)*x^2 + 75*Sqrt[195]*x^2 + (3500*I)*x^3))/12250000 Daniel Lichtblau Wolfram Research ==== On second thoughts: any one who really wants to see what the exact answer is, can evaluate the following: sols = Select[{a, b} /. SolveAlways[ x^6 + (9/14)*x^5 + (9/28)*x^4 + (3/35)*x^3 + (9/700)*x^2 + (9/8750)*x + 3/87500 == (x^2 + a*x + b)*(x^4 + c*x^3 + d*x^2 + e*x + f), x], FreeQ[N[#1], _Complex] & ] Times @@ (Map[x^2 + {x, 1}.# &, sols]) N[%] (0.010974992601737203 + 0.20255610310498245*x + x^2)* (0.020476912388332696 + 0.20476912388332683*x + x^2)* (0.15256133957420942 + 0.23553164887424147*x + x^2) Andrzej Kozlowski Yokohama, Japan http://www.mimuw.edu.pl/~akoz/ http://platon.c.u-tokyo.ac.jp/andrzej/ > Factoring without specifying the extension does not really make sense. > Of course Mathematica can easily factor yur polynomial into linear > factors over the complex numbers (with the help of Solve), but I > suspect you are really asking for is factoring over the reals. This is > harder and needs more human input. But anyway, Mathematica can do > this, or at least I have done it using Mathematica. In fact if you are > satisfied with a numerical answer Mathematica can do alone and in > seconds: > > In[1]:= > Simplify[N[x^6 + (9/14)*x^5 + (9/28)*x^4 + (3/35)*x^3 + (9/700)*x^2 + > (9/8750)*x + 3/87500]] > > Out[1]= > 1.*(0.010974992601737198 + 0.20255610310498295*x + > x^2)*(0.020476912388332692 + > 0.2047691238833268*x + x^2)*(0.15256133957420948 + > 0.23553191586883315*x + x^2) > > But I have in fact been foolish enough to compute the exact answer > too. I do not propose to post it here for it's absolutely horrible > (expressed in terms of Root objects) and quite useless. However if you > really want to see it I can send it to you privately. > > Andrzej Kozlowski > Yokohama, Japan > http://www.mimuw.edu.pl/~akoz/ > http://platon.c.u-tokyo.ac.jp/andrzej/ > > > > >> To: mathgroup@smc.vnet.net >> >> Can Mathematica factor the polynomial >> >> p1=x^6+9/14*x^5+9/28*x^4+3/35*x^3+9/700*x^2+9/8750*x+3/87500; >> >> without a priori knowledge of the Extension field? >> >> >> > > Reply-To: ==== Those aren't equations; they're functions, and we don't solve functions, we solve equations. What equations do you want to solve? DrBob -----Original Message----- and -3108 - 11 n >= 0 and solve for n, n <= -282.545 So we have -291.375 <= n <= -282.545. The integer solution set here is for n = {-290, -289, -288, -287, -286, -285, -284, -283}. So in this case we have 8 non-negative solutions. Given that I can supply x[n] and y[n], how do I go about finding the set n? ==== Is there a logically fundamental difference between functional and procedural programming? What I mean to ask is, can we do exactly the same thing with purely functional approaches as we can with purely procedural approaches? Is this basically the recursive verses iterative distinction? Why would one chose one approach over the other? STH . ==== > Let's be realistic. If you want 60 digits of precision, too bad! -- in > the real world. There's nothing we can measure that closely. Drug > concentrations in clinical trials are generally measured within 15%, for > instance. Even machine precision is more than can be realistically > expected in any application I can think of. Even getting a satellite to > Jupiter probably involves more error in the final result than machine > precision. (If not, it's because we rely on ongoing corrections and > natural factors that put the satellite where it should be, such as > gravity drawing it toward each rendezvous -- not on that kind of > precision in propulsion or guidance.) > > So... unless all numerics in a problem have a theoretical origin, and > could be represented in Mathematica as Infinite precision expressions... > all this talk of higher-precision computation seems futile. ...Except in the very rare instance when one needs to do intermediate calculations with, e.g., 60 digits of precision in order to get only a few correct digits of the final answer. The length of this thread is surely proof of the need for a definitive reference on the topic. Has there been a Mathematica-centered numerical analysis book published since Skeel & Keiper? --- Selwyn Hollis ==== > In[24]:= > f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - > 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 100]; > a=SetPrecision[77617.,100]; b = SetPrecision[33096.,100]; > > In[26]:= > {f, Precision[f]} > > Out[26]= > {-0.82739605994682136814116509547981629199903311578438481991 > 781484167246798617832`61.2597, 61} > > Congratulations! You just requested accuracy of > 100 for f and got 61 ( > to convince yourself add Accuracy[f] to In[26]). There is no request in that example for an accuracy of 100 in the result. The only request is for an accuracy of 100 in the input. > In[26]:= f = SetAccuracy[333.75*b^6 + > a^2*(11*a^2*b^2 - b^6 - > 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 100]; > > In[27]:= Accuracy[f] > Out[27]= 100. > > Now we assign values to some indeterminants in f. > > In[28]:= a = SetPrecision[77617.,100]; b = SetPrecision[33096.,100]; > > In[29]:= {f, Precision[f], Accuracy[f]} > Out[29]= > {-0.8273960599468213681411650954798162919990331157843848199178148, > 61.2599, 61.3422} > > The precision and accuracy has dropped. This is all > according to > standard numerical analysis regarding cancellation > error. You'll find it in any textbook on the topic. > > > Assume that I want accuracy and precision of 100 for > f. You advice me to make experiments to find out, what > should be the initial precision and accuracy of a and > b to reach the requested accuracy and precision for f. > Notice, that you cannot just repeat I[26], we saw > already what happens. I have to re-type I[24], I[25], > I[26], I[27], I[28], and I[29] as many times as needed > to get f with accuracy and precision 100. > > Dan, you simply advocate to do MANUAL WORK that should > be done by machine. You do not have to do any of this manually. The machine (Mathematica) will do all of this, usually using built-in functions. The N function, for example, will automatically adjust the working precision to give you a precision that you request, provided that doing so doesn't involve making up arbitrary digits. The example above starts out with machine numbers (333.75, 5.5, etc.), uses SetPrecision and SetAccuracy to make up arbitrary digits to pad those numbers out to some specified number of digits, and then does some simple arithmetic. If the goal is to get some specified number of digits in the result, and it is ok to make up arbitrary digits like this to achieve that goal, then the only manual work required to achive that goal is to apply SetAccuracy or SetPrecision to the result, to tell the computer that that is what you want. > Let's suppose that in the above example I just want 60 > digits not 61. Precisely, I want 60 digits and nothing > or zeros afterwards. Let's see if I could use > SetAccuracy. > > In[30]:= SetAccuracy[%, 60] > > Out[30]= -0.82739605994682136814116509547981629199903311578438481991781 > > In[31]:= % // FullForm > > Out[30]//FullForm= > -0.827396059946821368141165095479816291999033115784384819917814841672467988` > 59.9177 > > Oops, it did not work (as expected). If you could explain what you were expecting I am sure there are many contributors to this group who could explain to you why it did not do that. > Let's highlight > with mouse the expression in Out[30] and copy to a new > cell. Oops, we got > -0.827396059946821368141165095479816291999033115784384819917814841672467988` > 59.9177 again. Let's change Out[30] to a text cell and then copy. > > In[31]:= -0.82739605994682136814116509547981629199903311578438481991781 > > Out[31]= -0.82739605994682136814116509547981629199903311578438481991781 > > Success? Not so fast. If you could describe what you were trying to achieve with all of that copying and pasting and such I am again sure that there are many contributors to this group who could describe how to do it. It is very unlikely that the process will involve any copying and pasting or detours through text cells. > In[32]:= > % // FullForm > > Out[32]//FullForm= > -0.82739605994682136814116509547981629199903311578438481991780999999999999863 5 > 08`59.2041 > > Dan, is there any simple way to get what I want? Probably the answer is yes, but you will have to describe more clearly what you want. > As I repeated already number of times, at this stage > of the development of computer technology, software > should do it for me (!). If what you want to do is get a certain number of digits in the result, and it is ok to make up arbitrary digits as in the examples above, then you can do that by simply applying SetPrecision or SetAccuracy to the result. If you want the computer to automatically adjust the working precision to give a certain precision in the result, you can do that using N. If you want something else, and you can describe what that is, then probably someone can describe how to get Mathematica to do that for you. Dave Withoff Wolfram Research ==== Putting in Print[Some header text, TableForm[ Table[----------]]] puts the header text and the Table in the same cell, but with a largish gap (3 lines?) between the text and the Table's header lines. Hardly a serious problem, pretty trivial in fact, but a bit ugly and seems as if it really shouldn't work this way. Any way to get rid of this? siegman@stanford.edu Reply-To: ==== The first problem is that ArcSin[x/10] isn't a left-inverse of 10 Sin[x] in the region of the root. That's why your code converges on another root. A better inverse for your purpose is the gi below: f = #(# + 3) &; g = 10Sin@# &; gi = Pi - ArcSin[#/10] &; Using that inverse, there's still a problem, though. Your approximation of the root with r = gi@f@r fails because the derivative of gi@f at the root is gi'[f@rr]f'[rr] -2.02342 That's greater than one in magnitude, so distance from the root is magnified rather than diminished. Instead, try r = fi@g@r Where fi[y_] := Evaluate[x /. Last@Solve[f@x == y, x]] This should work, since fi'[g@rr]g'[rr] -0.49421355442685166 Sure enough, it does work: Values = NestList[fi[g[#1]] & , 2., 12] {2., 1.8679332339369226, 1.9368280058437026, 1.9040490698559083, 1.9205018812387413, 1.9124384783620436, 1.916439313493727, 1.9144659935649395, 1.9154421880638148, 1.914959973607326, 1.915198347545482, 1.9150805538598723, 1.9151387724997768} (f[#1] - g[#1] & ) /@ values {0.9070257317431825, -0.4688124734947827, 0.2242366717647286, -0.11228304957089463, 0.05509675095190758, -0.02732121417963107, 0.0134795615740817, -0.006667318794729482, 0.003293718665990042, -0.0016281317372435211, 0.0008045637255396088, -0.00039764607942949226, 0.0001965172491082967} Absolute error is cut in half at each iteration, as expected with a derivative for fi@g near -1/2. The negative sign causes the error to alternate in sign. We can also look at the log-absolute value of the error as follows: (Log[Abs[f[#1] - g[#1]]] & ) /@ values {-0.09758445910053692, -0.7575524337892089, -1.4950532145264737, -2.18673236744676, -2.8986645309519163, -3.6000918023816326, -4.306580698204814, -5.010537479670821, -5.715738058891756, -6.420322094992857, -7.125210383308283, -7.829948195959684, -8.534760348785936} Rest[%] - Drop[%, -1] {-0.659967974688672, -0.7375007807372648, -0.6916791529202861, -0.7119321635051565, -0.7014272714297163, -0.7064888958231816, -0.7039567814660064, -0.7052005792209357, -0.7045840361011004, -0.7048882883154262, -0.704737812651401, -0.7048121528262516} The difference in logarithm of the absolute error approaches Log@Abs[fi'[g@rr]g'[rr]] -0.7047875587967565 If you want to use this for teaching, try to use as little code as possible -- as I have above -- and always try to avoid Do loops in Mathematica on principle. DrBob -----Original Message----- kenf Below is the code Clear[f, g, gi, lim, r, rr, fr, gir, a, b, c, d, conv]; Plot[{x * ((x + 3)), 10*Sin[x]}, {x, 0.01, 2.4}, PlotStyle -> {{RGBColor[1, 0, 0], Thickness[ .006]}, {RGBColor[0, 0, 1], Thickness[ .006]}} ]; rr = FindRoot[x * ((x + 3)) == 10*Sin[x], {x, 2, 0.01, 2.4}]; f[a_] := a * ((a + 3)) /; a > 0; g[b_] := 10. * Sin[b] /; b > 0; gi[c_] := ArcSin[0.1*c] /; c > 0; Print[Actual root is , rr]; lim = 10; r = 2.0; conv = 10^-4; For[i = 1, i < lim, i++, { fr = f[r]; gir = gi[fr]; d = Abs[N[gir] - r]; i If[d < conv, Break[]]; r = gir; Print[The value of x = , r, found after , i, iterations,, with a tolerence , d, n] } ] Print[The value of x = , r, found after , i, iterations,, with a tolerence , d, n] Every man, woman and responsible child has an unalienable individual, civil, Constitutional and human right to obtain, own, and carry, openly or concealed, any weapon -- rifle, shotgun, handgun, machine gun, anything -- any time, any place, without asking anyone's permission. L. Neil Smith ==== f[t_] = {t + Sqrt[3] Sin[t], 2 Cos[t], t Sqrt[3] - Sin[t]} So It's basically a vector whose coordinates are determined based on the values you pass in. Then I took the derivative by just typing f', which outputs {1 + Sqrt[3] Cos[#1], -2 Sin[#1], Sqrt[3] - Cos[#1]}& What I'd like to do is have Mathematica calculate the norm of this as it would any vector, so that I can play with the norm function. As it turns out, the norm in this case is identical to Sqrt[8], so it would be nice if Mathematica could figure that out. Is it possible to do this? ==== [My previous reply seems to have gone astray - at least it has not come back to me. Here is a slightly edited repeat] You have computed f[t_] = {t + Sqrt[3] Sin[t], 2 Cos[t], t Sqrt[3] - Sin[t]}; and then found f' {1 + Sqrt[3]*Cos[#1], -2*Sin[#1], Sqrt[3] - Cos[#1]} & To get the function for the norm of the derivative we can use norm = Evaluate/@(Simplify/@(Sqrt[#.#]&/@(f'))) 2*Sqrt[2] & We map the usual functions for calclulating and simplifying the norm inside Function[.] (which is the full form of (.)& and then map the function Evaluation to make the result evaluate -- this is needed since Function has the attribute HoldAll. Please note that the parentheses are essential. -- Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 > > f[t_] = {t + Sqrt[3] Sin[t], 2 Cos[t], t Sqrt[3] - Sin[t]} > > So It's basically a vector whose coordinates are determined based on the > values you pass in. > > Then I took the derivative by just typing f', which outputs > > {1 + Sqrt[3] Cos[#1], -2 Sin[#1], Sqrt[3] - Cos[#1]}& > > > What I'd like to do is have Mathematica calculate the norm of this as it > would any vector, so that I can play with the norm function. As it turns > out, the norm in this case is identical to Sqrt[8], so it would be nice if > Mathematica could figure that out. Is it possible to do this? > > > ==== > f[t_] = {t + Sqrt[3] Sin[t], 2 Cos[t], t Sqrt[3] - Sin[t]} > > So It's basically a vector whose coordinates are determined based on the > values you pass in. > > Then I took the derivative by just typing f', which outputs > > {1 + Sqrt[3] Cos[#1], -2 Sin[#1], Sqrt[3] - Cos[#1]}& > > > What I'd like to do is have Mathematica calculate the norm of this as it > would any vector, so that I can play with the norm function. I don't think whether the definition is prompt (=) or delayed (:=) is as important as supplied the argument [t] in In := Simplify[Sqrt[f'[t] . f'[t]]] Out = 2 Sqrt[2] Or you can construct a new pure function with In := Evaluate[Simplify[Sqrt[f'[#] . f'[#]]]]& Out = 2 Sqrt[2] & Of course, there is not much left to play with :) Tom Burton ==== >I want to apply a function to every k-th element of a long list and >add the result to the k+1 element. > >[Actually k = 3 and I just want to multiply myList[[k]] by a >constant (independent of k) and add the result to myList[[k+1]] for >every value of k that's divisible by 3.] If I understand what you are trying to do correctly, the following should work (f[#[[1]]+#[[2]])&/@Partition[Drop[list,k-1],2,k] Here f is the function you want applied to the k-th element ==== MathGroup/Mathematica Newsgroup users - in the past 24 hours or so and have not seen your post, please resend/repost it. Disk drive problems caused some of the posts in the queue to be lost. This has been fixed. Sorry for the delay. Steve Christensen Moderator ==== >Please allow me to summarize what I've learned in the recent discussion, and >retract my claim that Accuracy, Precision, and SetAccuracy are useless. >Numbers come in three varieties The technical term is flavors. >- machine precision, Infinite precision, >and bignum or bigfloat. Bignums and bigfloats (synonymous?) Actually bignums can also refer to integers too large to represent as machine integers. But I tend to use bignum when I really mean bigfloat, and I suspect this sloppy practice may be common. >aren't called that in the Help Browser, but they're the result of using >N[expr,k] or >SetAccuracy[expr,k] where k is bigger than machine precision. >If k <= machine precision, the result is a machine precision number, even >if you know the expression isn't that precise. >If, when you use N or SetAccuracy as described above, the expression >contains undefined symbols, you get an expression with all its numerics >replaced by bignums of the indicated precision. When the symbols are >defined later, if ANY of them are machine precision, the expression is >computed with machine arithmetic - with the side-effect that coefficients >that originally were Infinite-precision are now only machine precision. >That is, x^2 might have become x^2.0000000000000000000000000000000000 >but later became x^2., for instance. I think this is correct in cases where all symbolic stuff gets replaced by numeric values. In general there is a sort of coercion to lowest precision, with the caveat that machine floats pollute everything. >If all the symbols have been set to bignum or Infinite precision values, >the computation will be done taking precision into account, and the result >has a Precision or Accuracy that makes sense. In all other cases, >Precision returns Infinity for entirely Infinite-precision expressions >and 16 for everything else. I'm not sure I understand this last sentence. My interpretation: Computations that are exact will have infinite precision. Computations in machine arithmetic will claim a precision of 16. If that is what you are claiming, then yes, that's what Mathematica is doing (but see my last remarks). >When one of the experts says significance arithmetic that's what they >mean - using SetAccuracy or N to give things more than 16 digits, leaving >no machine precision numbers anywhere in the expression, and using Accuracy >or Precision, which ARE meaningful in that case, to judge the result. >(It's meaningful if all your inputs really have more than 16 digits of >precision, that is.) I'm as guilty as anyone else in this thread, perhaps more so, of being too loose with the technical jargon. Also I am not certain what version 4 makes of SetAccuracy/SetPrecision in terms of significance arithmetic. In the development kernel they will force everything in sight to have the indicated precision, whether justified or not. This may well introduce error even with exact input, e.g. in cases where intermediate computations would require higher precision in order to get an end result with the requested precision or accuracy. N[], on the other hand, will handle that and, except in pathological circumstances, will give a result with the correct precision. As another minor point, arbitrary precision numbers are simply (tautologically?) numbers that may have arbitrarily large precision (subject to software limitations). Significance arithmetic refers to a particular model of manipulating such numbers with a mechanism for tracking precision. There are other models, in particular fixed precision arithmetic; that we use the former, by default, is an occasional source of sturm und drang in this news group. I'm sure the distinction between arbitrary precision numbers and significance arithmetic has at least minor relevance to this thread, and I imagine I've helped to confuse the issue in some places by using the terms almost interchangeably. >You can't use significance arithmetic to determine how much precision a >result has if your inputs have 16 or 15 or 2 digits of precision. One can, if the numbers are really bignums (of low precision, naturally). What one cannot do at present is create such low precision numbers via N[]. >In the example we've been looking at, you can give the inputs MORE accuracy >than you really believe they have, and still get back 0 digits from >Precision at the end, so there are clearly no trustworthy digits when >you use the original inputs either. If an expression is on the razor's >edge, and has lost only a few digits of precision, that wouldn't work >so well. >Oddly enough, significance arithmetic in the Browser doesn't take you >to any of that. Instead, it takes you to Interval arithmetic, a more >sophisticated method, which may give a more accurate gauge of how much >precision you really have, and WILL deal with machine precision numbers >and numbers with even less precision. It does a very good job on the >example. However, it isn't very suitable for Complex numbers, matrices, >etc. NSolve and NIntegrate probably can't handle it, either. I have filed a suggestion in-house that the documentation on significance arithmetic take one to the section on arbitrary precision numbers (3.1.5), as that would be more appropriate. Note that that section, while primarily concerned with the significance arithmetic model, also briefly mentions fixed precision bignum arithmetic. >Daniel Lichtblau promises that all this will be clearer in the next release. I'm not sure I'd go that far. What I will claim is that the distinction between machine numbers and bignums will be more transparent to users. At present if one does, say, N[number,precision] then one will get a machine number if prec<=$MachinePrecision. We have made a change so that this will no longer be the case. I am not prepared to go into details at this time (sorry). Perhaps more important for everyday use, and certainly more pertinent for this thread, Precision[] will distinguish between bignums of 16 digits precision and machine numbers. Again, I have to defer on details. At the very least I think the pitfall of believing a claim of 16 digits precision for machine numbers will be removed. Daniel Lichtblau Wolfram Research ==== Bobby, One point: >.... bigfloats ... [are] the result of using N[expr,k] or SetAccuracy[expr,k] where k is bigger than machine precision. If k <= > machine > precision, the result is a machine precision number. We get bigfloats with k<=machine precision with SetAccuracy and SetPrecision but not with N: Example a=SetPrecision[2.3,5] 2.3000 Precision[a] 5. Precision[a^2000] 1.69897 Also, of course, when more than machine precision significant digits are given a=1.01234567891234500; Precision[a] 17.301 Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 Reply-To: ==== Please allow me to summarize what I've learned in the recent discussion, and retract my claim that Accuracy, Precision, and SetAccuracy are useless. Numbers come in three varieties - machine precision, Infinite precision, and bignum or bigfloat. Bignums and bigfloats (synonymous?) aren't called that in the Help Browser, but they're the result of using N[expr,k] or SetAccuracy[expr,k] where k is bigger than machine precision. If k <= machine precision, the result is a machine precision number, even if you know the expression isn't that precise. If, when you use N or SetAccuracy as described above, the expression contains undefined symbols, you get an expression with all its numerics replaced by bignums of the indicated precision. When the symbols are defined later, if ANY of them are machine precision, the expression is computed with machine arithmetic - with the side-effect that coefficients that originally were Infinite-precision are now only machine precision. That is, x^2 might have become x^2.0000000000000000000000000000000000 but later became x^2., for instance. If all the symbols have been set to bignum or Infinite precision values, the computation will be done taking precision into account, and the result has a Precision or Accuracy that makes sense. In all other cases, Precision returns Infinity for entirely Infinite-precision expressions and 16 for everything else. When one of the experts says significance arithmetic that's what they mean - using SetAccuracy or N to give things more than 16 digits, leaving no machine precision numbers anywhere in the expression, and using Accuracy or Precision, which ARE meaningful in that case, to judge the result. (It's meaningful if all your inputs really have more than 16 digits of precision, that is.) You can't use significance arithmetic to determine how much precision a result has if your inputs have 16 or 15 or 2 digits of precision. In the example we've been looking at, you can give the inputs MORE accuracy than you really believe they have, and still get back 0 digits from Precision at the end, so there are clearly no trustworthy digits when you use the original inputs either. If an expression is on the razor's edge, and has lost only a few digits of precision, that wouldn't work so well. Oddly enough, significance arithmetic i