Jerry Avins wrote: (snip)> g[n] = (N^2/g[n-1] + g[n-1])/2> g[0] = 1 > g[1] = 50.00500 Far from convergence, convergence is slow. > g[2] = 26.24010 The error is approximately halved. > g[3] = 15.02553 > g[4] = 10.84043 Good to 2 digits. When close, convergence > g[5] = 10.03258 Good to 3 digits. is fast. > g[6] = 10.00005 Good to 6 digits. > 6[7] = 10.0000000001 11 digitsMore specifically, digits[n]=1-log10(abs(g[n]-10)) digits[0]=.046 digits[1]=-0.6 digits[2]=-0.2 digits[3]=0.3 digits[4]=1.08 digits[5]=2.49 digits[6]=5.3 digits[7]=11 -- glen
Fast Cube Root Using C33
Started by ●February 26, 2009
Reply by ●March 1, 20092009-03-01
Reply by ●March 2, 20092009-03-02
On Fri, 27 Feb 2009 06:54:02 -0800 (PST), pnachtwey <pnachtwey@gmail.com> wrote:>On Feb 26, 10:35�pm, Tim Wescott <t...@seemywebsite.com> wrote: >> On Thu, 26 Feb 2009 19:04:01 -0800, pnachtwey wrote: >> > I need a fast cube root routine. �I started with the obvious >> >http://metamerist.com/cbrt/cbrt.htm >> > What is key is getting the initial estimate using the Kahan bit hack but >> > that only works with IEEE floating point and not with TI floating point. >> > �I can't believe that I am the first to do this since the C33 is >> > ancient. �The obvious solution is to convert the TI float to the IEEE >> > float, apply the hack, and then convert back to TI float but there must >> > be a better way. � Actually doing another iteration is probably faster. >> > The problem is that sign of the mantissa is at bit 23. >> >> > Does anyboy know of a link? �Thanks. >> >> > Peter Nachtwey >> >> Dunno the details (I started doing DSP long after the C33 went obsolete, >> and I'm too cool* to use a floating point DSP anyway). �So I can't help >> you directly, but I _can_ jump in and kibbitz. >> >> Have you considered extracting the exponent, doing a quick table look-up >> on the resulting mantissa bits, and using whatever the table coughs up as >> a seed for N-R? >> >> * or perhaps my employers are too cheap. �Hmm. >> >> --http://www.wescottdesign.com >The exponent look up table is a good alternative because is quick but >a memory hog. I was looking for a variation of the Kahan "bit >hack". I might try that and use the mantissa to interpolate between >the values. At least I will know. I think this can be done without a >divide. > >Peter NachtweySorry if this is naive, but I would have imagined that just multiplying the exponent by 0.3333 would produce a decent initial value?
Reply by ●March 2, 20092009-03-02
Scott Hemphill wrote:> Martin Eisenberg <martin.eisenberg@udo.edu> writes:>> Additionally, you might gain from computing the *inverse* cube >> root because that has a polynomial Newton function, reducing >> the number of float divisions to one inversion at the end. > > No divisions are required at the end. To amplify, the function > for which we are trying to find a root is: > > f(x) = x^-3 - n> When we are done, we can calculate: > > n^(1/3) = n*x^2It actually depends on accuracy requirements. Say the relative error in x after the loop is at most R; then 1/x stays essentially as accurate, whereas the typical relative error in n*x^2 about doubles to 2*R. This is clearly most bothersome if one uses only a few iterations but it holds down to full accuracy (R = 1/2 ulp). Martin -- I used to think the brain was the most important organ in the body, until I realized who was telling me that. --Emo Phillips
Reply by ●March 2, 20092009-03-02
Tony wrote:> Sorry if this is naive, but I would have imagined that just > multiplying the exponent by 0.3333 would produce a decent initial value?Well, it is a floating point value, so there is an integer exponent (usually of two) and a significand (or fraction). You divide the integer exponent, then use the remainder to improve the initial value. There is a tradeoff between a better initial value and more cycles of Newton-Raphson. The optimal point is somewhat system dependent. -- glen
Reply by ●March 3, 20092009-03-03
On Mar 2, 1:21�am, Tony <t...@nowhere.com> wrote:> On Fri, 27 Feb 2009 06:54:02 -0800 (PST), pnachtwey <pnacht...@gmail.com> wrote: > >On Feb 26, 10:35�pm, Tim Wescott <t...@seemywebsite.com> wrote: > >> On Thu, 26 Feb 2009 19:04:01 -0800, pnachtwey wrote: > >> > I need a fast cube root routine. �I started with the obvious > >> >http://metamerist.com/cbrt/cbrt.htm > >> > What is key is getting the initial estimate using the Kahan bit hack but > >> > that only works with IEEE floating point and not with TI floating point. > >> > �I can't believe that I am the first to do this since the C33 is > >> > ancient. �The obvious solution is to convert the TI float to the IEEE > >> > float, apply the hack, and then convert back to TI float but there must > >> > be a better way. � Actually doing another iteration is probably faster. > >> > The problem is that sign of the mantissa is at bit 23. > > >> > Does anyboy know of a link? �Thanks. > > >> > Peter Nachtwey > > >> Dunno the details (I started doing DSP long after the C33 went obsolete, > >> and I'm too cool* to use a floating point DSP anyway). �So I can't help > >> you directly, but I _can_ jump in and kibbitz. > > >> Have you considered extracting the exponent, doing a quick table look-up > >> on the resulting mantissa bits, and using whatever the table coughs up as > >> a seed for N-R? > > >> * or perhaps my employers are too cheap. �Hmm. > > >> --http://www.wescottdesign.com > >The exponent look up table is a good alternative because is quick but > >a memory hog. �I was looking for a variation of the Kahan "bit > >hack". � I might try that and use the mantissa to interpolate between > >the values. �At least I will know. �I think this can be done without a > >divide. > > >Peter Nachtwey > > Sorry if this is naive, but I would have imagined that just multiplying the exponent by > 0.3333 would produce a decent initial value?Basically that is what I did. The TI C33 has a fast_itof intrinsic instruction to convert the exponent back to an integer. The problem is that the mantissa gets messed up so I cleared those bits. I also added an offset to the exponent before multiplying by 0.3333333. This reduced the max error. There wasn't much math behind what I did. Tried different methods and tried to minimize the worst case. I just looking for quick and dirty because the iterations are cheap. Peter Nachtwey
Reply by ●March 3, 20092009-03-03
pnachtwey wrote: ...> Basically that is what I did. The TI C33 has a fast_itof intrinsic > instruction to convert the exponent back to an integer. The problem > is that the mantissa gets messed up so I cleared those bits. I also > added an offset to the exponent before multiplying by 0.3333333. This > reduced the max error. There wasn't much math behind what I did. > Tried different methods and tried to minimize the worst case. I just > looking for quick and dirty because the iterations are cheap.If you post the code, I'll copy it and thank you. I have no use for it now, but I do have 'C32 and 'C33 DSKits. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
