/* file: zeros.c */
/* does: finds zeros in interval; puts them into linked list zeroLink. */
/************************* HEADER ***************************************/
#include <stdio.h>
#include <math.h>
/*** ZEROLINK *** */
struct zerolink{
double val;
double divergingUpperVal;
int isZero;
struct zerolink *next;
struct zerolink *prev;
};
typedef struct zerolink zeroLink;
zeroLink* newZeroLink(); /* creates link element, memory alloc. */
zeroLink* weedZeroLink(zeroLink*,int); /* can be used in conjunction with
getPointsOpenInterval(..) in
order to sort qualified zeros.*/
zeroLink* firstLink(zeroLink*);
zeroLink* lastLink(zeroLink*);
zeroLink* trashLink(zeroLink*);
/* DEFINE a DATATYPE FOR double-double FORTRAN type FUNCTION */
typedef double ddfortran(double*);
/* GETPOINTSINTERVAL (FOR USER) */
zeroLink* getPointsInterval(ddfortran*, double, double, double, double);
/* GETPOINTSINTERVAL RELIES ON THE FOLLOWING BELOW: */
/* (not really intended for user use) */
/* SLATEC DFZERO */
void dfzero(ddfortran*, double*, double*, double*,double*,double*,int*);
/* MODIFIED DFZERO */
int dfZero( ddfortran*, double*, double*, double, double, int );
/* FURTHER MODIFIED DFZERO (returns link-lst-elt) */
zeroLink* getPoint(ddfortran*,double,double,double,double);
/* USES GETPOINT(..) */
void getPointsOpenInterval(ddfortran*,zeroLink*,double,double,double,double);
/* EPSILONZERO */
double epsilonZero(int, ddfortran*, zeroLink*, double, double);
/* suggests a small distance about calculated zero, on the order of a few RE's
(relative errors), which is used to create an small neighborhood around
calculated zero in which program ASSUMES no other zero to occur in. The
boundary of this interval is picked so this point is known not to qualify as
a calculated zero by having the value of the function at this point>=AE.
*/
/************************************************************************/
#define JDB if (DEBUG_MODE>0)
int DEBUG_MODE=0;
/**************** FUNCTION IMPLEMETATIONS *******************************/
zeroLink* newZeroLink() {
zeroLink* ldp;
ldp=(zeroLink*) malloc(sizeof(zeroLink));
ldp->next=0;
ldp->prev=0;
ldp->val=0.;
ldp->isZero=0;
return(ldp);
}
zeroLink* weedZeroLink(zeroLink* k, int ae_per_billion) {
if (k==0) return(k);
for(k=firstLink(k); k->next ; k=k->next)
if (k->next->next) {
if ((k->next->isZero==0)||(k->next->isZero < -ae_per_billion))
trashLink(k->next);
}
if ((k->isZero==0)||(k->isZero < -ae_per_billion)) return(trashLink(k));
else return(k);
}
zeroLink* firstLink(zeroLink* k)
{ if (k) while(k->prev) k=k->prev; return(k); }
zeroLink* lastLink(zeroLink* k)
{ if (k) while (k->next) k=k->next; return(k); }
zeroLink* trashLink(zeroLink* k)
{ zeroLink* l;
l=k;
if (k==0) return(k);
if (k->prev) k->prev->next=k->next;
if (k->next) {k->next->prev=k->prev; k=k->next;}
else k=k->prev;
free(l);
return(k);
}
/***/
double epsilonZero(int plusminus, ddfortran* pf, zeroLink* z, double ae, double re)
{
int TIMEOUT=0;
double x,y,dx;
x=y=z->val;
dx=plusminus*re;
for (x+=dx; (fabs((*pf)(&x))<=20*ae) ;x+=dx) {
TIMEOUT++;
if (TIMEOUT==100) puts(
"Warning:epsilonZero(). Resolution within relative error may not be possible.\n"
);
if (TIMEOUT%100==0) dx+=(TIMEOUT/100)*re;
if (TIMEOUT>1000)
{
puts("Serious ERROR. Too many contiguous points qualify as zeros.\n");
break;
}
}
JDB printf("E%d ",TIMEOUT);
return(x);
}
/***/
int dfZero(ddfortran* pfun, double* a, double* b, double re, double ae, int persist) {
int flag,i,jk,ji,jn,k;
double guess, A,B,dx;
guess=*a; A=*a; B=*b; jn=5;
if (persist>0)
{
persist=2*persist+1;
for(k=3;k<=persist;k+=2)
{
dx=(B-A)/k;
for (i=0;i<k;i++)
{ JDB printf(",");
*a=A+i*dx;
*b=A+(i+1)*dx;
guess=(*a)+dx/3.;
dfzero(pfun, a, b, &guess, &re, &ae, &flag);
if (flag<=jn) {jn=flag; jk=k;ji=i;}
if (flag<=1) return(flag);
}
}
dx=(B-A)/jk; *a=A+ji*dx; *b=A+(ji+1)*dx; guess=(*a)+dx/3.;
dfzero(pfun, a, b, &guess, &re, &ae, &flag);
if(flag!=jn) puts("dfZero:intel inside?\n");
return(flag);
}
else
{
dfzero(pfun, a, b, &guess, &re, &ae, &flag);
if (flag==5) {
puts("a dfzero has crashed while evaluating the interval ");
printf("(,%g,%g).\n Any zeros in that interval could be ignored.",A,B);
}
}
return(flag);
}
/***/
#define PERSIST 9
/* _PERSIST_ IS A NUMBER THAT IS USED IN THE CASE WHEN NO ZERO IS FOUND BY
DFZERO. TO ENSURE THAT THERE REALLY IS NO ZERO IN THE INTERVAL, WE RETEST BY
CALLING DFZERO REPETITIVELY, FOR NUMEROUS INTERVALS THAT MAKE UP THE
INTERVAL IN QUESTION */
zeroLink* getPoint(ddfortran* pf,double a,double b,double ae,double re) {
int i,persist=PERSIST;
zeroLink* z;
if (a>b) return(0);
JDB printf("Call<%g|%g>",a,b);
if (a==b) i=4;
else
i=dfZero(pf, &a, &b, ae,re,persist); JDB printf(" %d\n",i);
if (i==4) if (fabs((*pf)(&a))>=ae) i=6;
/*** Let flag=6 stand for: "not really a zero" ***/
/* We don t know how to handle flag=5 (out of cpu time in dfzero yet)
This can be implemented later; right now if flag=5 assume no zeros exist
*/if (i==5) {i=6; printf("dfzero.flag=5 dfzero(%g,%g) gave up",a,b);}
if (i==6) return((zeroLink*) 0);
/* IMPROVE ACCURACY */
dfZero(pf,&a,&b,ae/64,re/16,0);
z=newZeroLink();
z->val=a;
z->divergingUpperVal=b;
/************************* CLASSIFICATION **********************/
/* isZero==0 => point of diverging f; not really a zero. */
/* isZero==1 =>healthy zero. (simple zero or odd multiplicity) */
/* isZero <0 =>possible zero even multiplicity. */
/****************************************************************/
if (i==4) z->isZero= -1-((int) (1000000000*(*pf)(&a)/ae));
/* CLASSIFY case as a (questionable/possible) zero of EVEN MULTIPLICITY.
isZero is set negative; the magnitude of isZero indicates the error,
in parts per billion of abs err. */
if (i==3) z->isZero= 0; /* CLASSIFY this as a point where f DIVERGES */
if (i==2) z->isZero= 2; /* CLASSIFY this as an exact zero (-> val=0) */
if (i==1) z->isZero= 1; /* CLASSIFY this as a HEALTHY zero */
return z;
}
/***/
#define PP 3.1415927
/* PP IS ONLY USED FOR DEBUGGING PURPOSES */
/* GIVEN A ZERO WITIN AN INTERVAL OF TWO SPECIFIED ENDPOINTS, BOTH NOT ZERO
(NOR WITHIN ABS ERROR OF ZERO), FIND (RECURSIVELY) ALL ZEROS AND POSSIBLY
DIVERGING POINTS IN THAT INTERVAL. POINTS WILL BE STORED IN A LINKED LIST
LINKED TO Z. FOR DIVERGING POINTS zeroLink.isZero IS SET ZERO. */
void getPointsOpenInterval(ddfortran* pf,zeroLink* z,
double a,double b,double ae,double re)
{
double zminus,zplus, dx,a1,b1,a2,b2;
int i,j;
zeroLink *z1, *z2;
a1=a; b2=b; dx=re;
if (z->isZero) { zminus=epsilonZero(-2,pf,z,ae,re); /* resolution=2*re */
zplus=epsilonZero(+2,pf,z,ae,re);
} else { while(z->val==z->divergingUpperVal)/*handle eps.div.pts*/
{z->val-=dx;z->divergingUpperVal+=dx; dx+=re; }
zminus=z->val;
zplus=z->divergingUpperVal;
} /* search for zeros in the two intervals below and above z->val */
b1=zminus;
a2=zplus;
z1=getPoint(pf,a1,b1,ae,re);
z2=getPoint(pf,a2,b2,ae,re);
JDB if (z1) printf(" (z1[%g]%d) ",z1->val/PP,z1->isZero);
JDB if (z2) printf(" (z2[%g]%d) ",z2->val/PP,z2->isZero);
if (z1)
{
z1->prev=z->prev;
z1->next=z;
if (z->prev) z->prev->next=z1;
z->prev=z1;
}
if (z2)
{
z2->next=z->next;
z2->prev=z;
if (z->next) z->next->prev=z2;
z->next=z2;
}
JDB if (z1) if (z1->prev) if (z1->next)
printf(" <%g<%g>%g> ",z1->prev->val/PP,z1->val/PP,z1->next->val/PP);
JDB if (z2) if (z2->prev) if (z2->next)
printf(" <%g<%g>%g> ",z2->prev->val/PP,z2->val/PP,z2->next->val/PP);
JDB if (z) if (z->prev) if (z->next)
printf("------<%g<%g>%g> ",z->prev->val/PP,z->val/PP,z->next->val/PP);
/* now start BRANCHING. Calling recursively if not finished: */
if (z1) getPointsOpenInterval(pf, z1, a, zminus,ae,re);
if (z2) getPointsOpenInterval(pf, z2, zplus, b, ae,re);
JDB if (z) if (z->prev) if (z->next)
printf("\n<%g<%g>%g>\n",z->prev->val/PP,z->val/PP,z->next->val/PP);
}
/***/
/* GETPOINTSINTERVAL */
zeroLink* getPointsInterval(ddfortran* pf,double a,double b,double absE,double relE) {
double c,d;
zeroLink *middle, *enda, *endb, *tmp;
tmp=newZeroLink();
tmp->val=a;
c=epsilonZero(+2,pf,tmp,absE,relE);
enda=getPoint(pf, a, c, absE,relE);
tmp->val=b;
d=epsilonZero(-2,pf,tmp,absE,relE);
endb=getPoint(pf, d, b, absE,relE);
middle=getPoint(pf, c,d,absE,relE);
JDB if (middle) printf("SEED IS %g. (%d)", middle->val, middle->isZero);
if (middle) getPointsOpenInterval(pf,middle, c,d,absE,relE);
middle=firstLink(middle);
if ((enda==0)&&am(endb==0)) { free(tmp); JDB puts("REGular\n"); return(middle);}
else;
if (enda==0) enda=tmp;
if (endb==0) endb=tmp;
if (middle==0) { enda->next=endb; endb->prev=middle=enda; }
else {
enda->next=middle;
if (enda!=tmp) middle->prev=enda;
endb->prev=middle=lastLink(middle);
if (endb!=tmp) middle->next=endb;
}
middle=firstLink(middle);
free(tmp);
return(middle);
}