/* * findroot.c * * Copyright (C) 1989, 1991, Mark Podlipec, Craig E. Kolb * All rights reserved. * * This software may be freely copied, modified, and redistributed * provided that this copyright notice is preserved on all copies. * * You may not distribute this software, in whole or in part, as part of * any commercial product without the express consent of the authors. * * There is no warranty or other guarantee of fitness of this software * for any purpose. It is provided solely "as is". * * $Id$ * * $Log$ * $Revision$ * Initial version. * */ #include "libcommon/common.h" /* * POLY_SIZE should be set equal to or larger than the highest degree * polynomial to be solved by this routine. This structure could've * been an argument to FindRoot, thus making it more flexible, but * I decided it would be easier just to have it here and not have to worry * about it in all the routines that call it. * * Some space optimization could be done, but I didn't think the savings * would be worth the effort. * */ #define POLY_SIZE 8 typedef struct STRUCT_POLY_STACK { Float roots[POLY_SIZE]; int valid[POLY_SIZE]; Float f[POLY_SIZE]; } POLY_STACK; static POLY_STACK gg[POLY_SIZE]; static int r_loop=0; #define TRUE 1 #define FALSE 0 #define VALID 1 /* #define FR_EPS 1e-8 #define FR_EPS_2 1e-9 */ #define FR_EPS 1e-7 #define FR_EPS_2 1e-8 #define FR_EPS_3 1e-24 #define FRIsZero(x) ( ((x) > -FR_EPS) && ((x) < FR_EPS) ) #define FRIsZero3(x) ( ((x) > -FR_EPS_3) && ((x) < FR_EPS_3) ) /* * This is used to determine the sign of x */ #define FR_SGN(x) ( ((x)<0.0)?-1:(((x)>0.0)?1:0) ) /* * Evaluate the degree _Deg polynomial _p at point _x and put * the value into _y(note: when calling give _y not &_y) */ #define FR_EVALPOLY(_p,_Deg,_x,_y) \ { int _i; \ _y = (_p)[(_Deg)]; \ for(_i=((_Deg)-1);_i>=0;_i--) _y=(_x)*_y + (_p)[_i]; \ } #ifndef AMIGA Float sqrt(); #endif /* * This is a specialized routine to find the smallest root in the * interval a,b. FindRoot will return a 1 if a root is found and * a 0 if no roots are found. If a root is found, it's value will * be put in *root. */ int FindRoot(ff,degree,a,b,root,flag) Float *ff; /* array of coefficients of the polynomial */ /* ex: 5x^2 - 2x + 3.0 ff[0] = 3.0 ff[1] = -2.0 ff[2] = 5.0 */ int degree; /* the degree of polynomial */ Float a,b; /* the interval to look for roots in */ Float *root; /* where to put a root if found */ int flag; /* Should always be -1 when called. This is used to initialize derivatives the 1st time called but not when FindRoot calls itself subsequently.*/ { int i,loop; Float *f1,*f0; /* * Just in case something really goes wrong bomb out before stack space * gets recursivley eaten away. (This happened during initial debug and * I can't think a situation that would cause it now but it doesn't hurt. * * r_loop should be larger than maximum degree squared. */ if (r_loop > (POLY_SIZE * (POLY_SIZE+1)) ) { fprintf(stderr,"FindRoot Error: r_loop=%ld degree=%ld flag=%ld\n", r_loop,degree,flag); exit(0); } /* not a valid polynomial as far as we're concerned. */ if (degree==0) return FALSE; /* Some useful variables. f0=ff() and f1=ff'(). */ f0 = gg[degree].f; f1 = gg[degree-1].f; /* * Find Root called for the first time. * Calculate all derivatives of the polynomial. * and initialize static structure to hold them. */ if (flag==-1) { register int i,j; /* Initialize valid flags to !VALID */ for(i=0;i<8;i++) for(j=0;j<8;j++) gg[i].valid[j]=0; /* Initialize and calculate coefficients for all * derivatives of ff() */ for(i=0;i<=degree;i++) f0[i] = ff[i]; for(i=degree-1;i>=2;i--) for(j=0;j<=i;j++) gg[i].f[j] = (Float)(j+1) * gg[i+1].f[j+1]; flag = degree; /* initialize variable that flags runaway recursion. Again * I don't believe it will occur, but I won't rule it out. */ r_loop=1; } else /* * Not the 1st time through. Check to see if we've already calculated * a valid root for this derivative in this interval from a previous * recursive call. */ { register int ix; r_loop++; ix = 0; /* Check the valid flags and if valid, check to see if it * is in the current interval */ while( (gg[degree].valid[ix]==VALID) ) { if ( (gg[degree].roots[ix] > a ) && (gg[degree].roots[ix] < b ) ) { *root = gg[degree].roots[ix]; return TRUE; } ix++; /* * If all possible roots at this degree have been * found and none of them are in the interval [a,b] * then return FALSE. */ if (ix >= degree) return FALSE; } /* end of search for valid roots */ } /* end of not first time that FindRoot was called */ /* * If the leading coeff is very close to zero, approximate by reducing * the degree and recalling FindRoot. Prevents explosions. * * blob.c currently could do this since it constructs polynomials from * sums of other polynomials and not necessarily end up with degree 4. * * Check to see if this is ever called. TEST */ if ( FRIsZero3(gg[degree].f[degree]) ) return( FindRoot(f0,(degree-1),a,b,root,flag) ); /****************************************************************** * DEGREE==1 * * If ff() is of degree 1 then it's a simple line. */ if (degree==1) { if ( FRIsZero(f0[1]) ) return FALSE; else *root = -f0[0]/f0[1]; gg[degree].valid[0] = 1; gg[degree].roots[0] = *root; if ( (*root > a) && (*root < b) ) return TRUE; else return FALSE; } /****************************************************************** * DEGREE==2 * * If ff() is of degree 2 then use Quadratic Equation to find roots. */ if (degree==2) { register Float t,d0,d1; d0 = f0[1]; t = d0*d0 - 4.0 * f0[2]*f0[0]; /* Return FALSE if roots are imaginary. */ if (FR_SGN(t) == -1) return FALSE; /* There is only one root. */ if (FR_SGN(t) == 0.0) { /* Note f0[2] was checked for zero above * as gg[degree].f[degree]. */ d1 = -d0/(2.0 * f0[2]); /* Since for this set of coeff we found * all possible roots, set up static structure * to indicate this. */ gg[degree].valid[0] = VALID; gg[degree].valid[1] = VALID; gg[degree].roots[0] = d1; gg[degree].roots[1] = d1; if ( (d1>a) && (d1 0 which means there are two roots */ t = sqrt(t); /* If the smallest root has already been found, then it * was not in the current interval [a,b]. So skip this * and calculate the larger quadratic root. * NOTE: if both of the roots are valid then we wouldn't * be here, we would've either returned which ever one was * the smallest in the interval [a,b] or returned false. */ if (gg[degree].valid[0]!=VALID) { /* Calculate smallest Quadratic root */ if (f0[2]>0.0) d1= (-d0-t)/(2.0*f0[2]); else d1= (-d0+t)/(2.0*f0[2]); /* put it in the static structure */ gg[degree].valid[0] = VALID; gg[degree].roots[0] = d1; if ((d1>a) && (d10.0) d1= (-d0+t)/(2.0*f0[2]); else d1= (-d0-t)/(2.0*f0[2]); /* put it in the static structure */ gg[degree].valid[1] = VALID; gg[degree].roots[1] = d1; if ((d1>a) && (d1=3 * * Now for the fun stuff. * * The algorithm being find the smallest root of f'() in the interval * [a,b]. This will either be a min or a max(we don't care which). Call * this point c. Now evaluate f(a) and f(c). If they are of different * signs, the f() has a root in the interal [a,c], use Newton's method * to solve. If they are of the same sign, there is no root in [a,c] * so move the interval in question to [c,b] and recall FindRoot. * Oh, and if there is no min/max then the sign tests apply with a and b * in the same way as the applied to a and c. * * There's some tricky parts, like what happens if there's a root at a? * Since this is raytracing, we assume the ray was fired off from a * surface and we didn't quite clear it. So we walk a forward in * increments of FR_EPS2 which is smaller than FR_EPS. If after an * arbitrary amount, f(a) is still zero, then we return the root at a+. * Another situation where this occurs, is when we move the interval * from [a,b] to [c,b] and then search for a min/max along the new * interval, f'(c) is going to zero and we want the next root, not * the one we've already found. So we need to walk the start of the * interval in this case as well. * Walking might be an overkill, it might only be necessary to jump * one arbitrary amount(EPS or 2*EPS). I still need to experimentally * decide which is more efficient. * */ { Float y0,y1,c; int i,diff_ret; /* find the next root of f'() and store in c */ diff_ret = FindRoot(f1,(degree-1),a,b,&c,flag); /* Evaluate the polynomial at the beginning of the interval * and walk it forward if it is zero */ FR_EVALPOLY(f0,degree,a,y0); loop=0; while( (FRIsZero(y0)) && (loop<=10) ) { a+=FR_EPS_2; FR_EVALPOLY(f0,degree,a,y0); loop++; } /* * No roots of ff'() in [a,b] so there is no min/max */ if (!diff_ret) { /* Evaluate polynomial at point b */ FR_EVALPOLY(f0,degree,b,y1); /* if ff(a) and ff(b) have the same sign return false * else use Newton */ if (FR_SGN(y0) == FR_SGN(y1)) return FALSE; return( Newton(f0,f1,degree,a,b,root) ); } /* * ff'() does have a root at c where a < c < b. */ FR_EVALPOLY(f0,degree,c,y1); /* * if same sign move interval to [c,b] and recall FindRoot * else there is a root in [a,c] use Newton to solve. */ if (FR_SGN(y0) == FR_SGN(y1)) return(FindRoot(f0,degree,c,b,root,flag)); else return( Newton(f0,f1,degree,a,c,root) ); } /* end of degree >= 3 */ } /* end of routine */ /* * This routine numerically tries to find a root of f0() in the interval * [a,b]. */ int Newton(f0,f1,degree,a,b,root) Float *f0; /* f0() is the polynomial */ Float *f1; /* f1() is the 1st derivative of the polynomial */ int degree; /* degree is the degree of the polynomial f0() */ Float a,b; /* this is the interval [a,b] to find the roots on */ Float *root; /* if a root is found, put here */ { register int i,loop_num; register Float y0,y1,x0; /* loop an arbitrary number of times, if no root is found then * give up. In experiments with cubicspin and blob code, The root * was either resolved in <18 iterations or it was never(50+) * resolved. */ loop_num = 50; /* do a bisection to start with in case slope is ~0 at a. * this is likely since the interval moves to the min/max's * of the derivatives. */ x0 = 0.5 * (a + b); do { FR_EVALPOLY(f0,degree,x0,y0); /* y0 = f(x0) */ /* check to see if we are close enough to zero to * call it a success */ if ( FRIsZero(y0) ) { /* is it within range? This might be redundant. */ if ( (x0>(a+FR_EPS_2)) && (x0 b) tmp_x=(x0+b)/2.0; else if (tmp_x < a) tmp_x=(x0+a)/2.0; x0 = tmp_x; } loop_num--; } while(loop_num); return FALSE; } /* end of Newton */