Yeah, thanks, i knew about that one. Only .....
a float divide or float multiply takes ~ 400 cycles on a PIC chip. At 20 MHz, doing it that way takes LOTS AND LOTS of time.... but then again, i only need to check it three times a second.... i''ll run sims on it , but i have a feeling it will take WAY too long.
Thanks very much for ur time. I know it is not game related... sorry :|
ARRRRRRRRRRGH - horrible (and possibly unsolvable ) problem
quote:
pow = exp(b * x0);
fx = i + (-b * a * x0 * pow + a * (1 - pow));
fpx = -2 * b * a * pow - b * b * a * x0 * pow;
x0 -= fx / fpx;
One can still shave those formulas a lot. Instead of solving x0, solve b*x0 and divide with b once you get the result. Above has 12 multiples and 1 divide per loop, but it can be reduced to (2 divs and 1 multi) or (1 div and 3 multi). For the exp() part, I suggest using a small table and some interpolation that produces good results in tests (try different ones, e.g. cubic interpolation). It's a pretty well behaving function so I'd expect to get good results that way, in ridiculously low time compared to McLaurin series.
[edited by - civguy on March 5, 2004 3:58:01 AM]
I just want to add this note...
If you need to calculate it 3 times per second, you will use the last value as the new iteration startpoint, eg. your value is likely to be accurate and just need some adjustment and thus wouldn''t take up that much cpu time. (if the changes of the parameters aren''t large).
But you are probably doing that already...
If you need to calculate it 3 times per second, you will use the last value as the new iteration startpoint, eg. your value is likely to be accurate and just need some adjustment and thus wouldn''t take up that much cpu time. (if the changes of the parameters aren''t large).
But you are probably doing that already...
behold the power of Maple:
and all that in just 3.7 seconds.
but it seems to me by the way this lambertw function is defined, ie as a function of itsself, its just shoving the problem away. look it up id say.
[edited by - eelco on March 5, 2004 12:22:31 PM]
solve( I+(-b*a*v*exp(b*v)+a*(1-exp(b*v))) = 0,v); (I + a) exp(1) LambertW(--------------) - 1 a ---------------------------- b and all that in just 3.7 seconds.
but it seems to me by the way this lambertw function is defined, ie as a function of itsself, its just shoving the problem away. look it up id say.
[edited by - eelco on March 5, 2004 12:22:31 PM]
Author
(I + a) exp(1) LambertW(--------------) - 1 a ---------------------------- bSo if i solve this, it will give me an exact answer? Rooting around on the internet gives me
/* ANSI C code for W(x). K M Briggs 98 Feb 12 http://www-epidem.plantsci.cam.ac.uk/~kbriggs/LambertW.c Based on Halley iteration. Converges rapidly for all valid x. */double W(const double x) { int i; double p,e,t,w,eps=4.0e-16; /* eps=desired precision */ if (x<-0.36787944117144232159552377016146086) { fprintf(stderr,"x=%g is < -1/e, exiting.\n",x); exit(1); } if (0.0==x) return 0.0; /* get initial approximation for iteration... */ if (x<1.0) { /* series near 0 */ p=sqrt(2.0*(2.7182818284590452353602874713526625*x+1.0)); w=-1.0+p-p*p/3.0+11.0/72.0*p*p*p; } else w=log(x); /* asymptotic */ if (x>3.0) w-=log(w); for (i=0; i<20; i++) { /* Halley loop */ e=exp(w); t=w*e-x; t/=e*(w+1.0)-0.5*(w+2.0)*t/(w+1.0); w-=t; if (fabs(t)<eps*(1.0+fabs(w))) return w; /* rel-abs error */ } /* never gets here */ fprintf(stderr,"No convergence at x=%g\n",x); exit(1);}.... i''ll try and work it out on computer... but it looks nice... (only same problem with log and exp... )
Lookup table is looking good at the moment...
Thanks for all your help
Author
SOLVED!!!!!!!!!!!!
Thankyou guys (and gals) VERY VERY VERY much for all ur help....
Cheerz
JB
//#include "math.h"int i;long l;#define M_LN2 0.69314718055994530942double scratch;#define M_EXP1 2.718281828459045/* (I + a) exp(1) LambertW(--------------) - 1 a ---------------------------- b */// ret = b^x, x (= {J};double ipow(double b, int x) { //Move into temp / "scratch" register scratch = b; //Iterate for (i = 1; i < x; i++) //Multiply scratch*=b; //Return calculated value (b^x, x (= {J}; ) return scratch;}double factorial(int x) { scratch = x--; while (x>0) scratch *= (x--); return scratch;}double exp(double x) { long ipart; double exponent,answer; //e^x = 2^(x/ln(2)); //Find the integer component ipart = (long)x;/* //Work out 2 to the ipart... ipart = 1<<ipart;*/ //Ok, work out x / ln(2). exponent = x / M_LN2; //Work out floating point part exponent = x - (long)x; //Set to 1 (e factorial shite) answer = 1; //Decrement i for (i =1; i < 6; i++) { answer += ( ipow(exponent,i) / factorial(i) ); }; return ipow(M_EXP1,ipart) * answer;}double fabs(const double x) { //Is x smaller than zero (ie. negative) if (x < 0.0f) //Multiply by negative one to make it positive return x*-1; //Just plain ol'' x return x;}/* ANSI C code for W(x). K M Briggs 98 Feb 12 http://www-epidem.plantsci.cam.ac.uk/~kbriggs/LambertW.c Based on Halley iteration. Converges rapidly for all valid x. */double LambertW(const double x) { int i; double p,e,t,w,eps=1e-3; /* eps=desired precision */ if (x<-0.36787944117144232159552377016146086) {return(-1); } //Let''s not root around... get it ??? w = 1; for (i=0; i<20; i++) { /* Halley loop */ e=exp(w); t=w*e-x; t/=e*(w+1.0)-0.5*(w+2.0)*t/(w+1.0); w-=t; if (fabs(t)<eps*(1.0+fabs(w))) return w; /* rel-abs error */ } }int main( ) { /* (I + a) exp(1) LambertW(--------------) - 1 a ---------------------------- b */#define I 0.65f#define a 0.000001f#define b 0.65f double v; //Solve equation v= (LambertW( (I+a)*M_EXP1/a) - 1)/b; return v;}Thankyou guys (and gals) VERY VERY VERY much for all ur help....
Cheerz
JB
GOOD! I''m glad you found a solution. Now, since this was confirmed not to be game development related, and since the poster has his answer, I''m closing the thread!
Graham Rhodes
Senior Scientist
Applied Research Associates, Inc.
Graham Rhodes
Senior Scientist
Applied Research Associates, Inc.
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
