diff options
author | Norbert Preining <norbert@preining.info> | 2019-09-02 13:46:59 +0900 |
---|---|---|
committer | Norbert Preining <norbert@preining.info> | 2019-09-02 13:46:59 +0900 |
commit | e0c6872cf40896c7be36b11dcc744620f10adf1d (patch) | |
tree | 60335e10d2f4354b0674ec22d7b53f0f8abee672 /support/fastpictex/src/myfunctions.cc |
Initial commit
Diffstat (limited to 'support/fastpictex/src/myfunctions.cc')
-rw-r--r-- | support/fastpictex/src/myfunctions.cc | 211 |
1 files changed, 211 insertions, 0 deletions
diff --git a/support/fastpictex/src/myfunctions.cc b/support/fastpictex/src/myfunctions.cc new file mode 100644 index 0000000000..1bfb73b6b1 --- /dev/null +++ b/support/fastpictex/src/myfunctions.cc @@ -0,0 +1,211 @@ +#include <iostream> +#include <stdlib.h> +#include <string.h> +#include <math.h> + +int myatoi(char *s, int *i) { + float wert; + int rc=0; + char *endp; + + wert=strtod(s, &endp); + if (s==endp) { + rc=0; + } + else { + *i=(int)floor(wert); + rc=1; + } + return(rc); +} + +int myatof(char *s, float *f) { + double wert; + int rc=0; + char *endp; + + wert=strtod(s, &endp); + if (s==endp) { + rc=0; + } + else { + *f=(float)wert; + rc=1; + } + return(rc); +} + +char *getsig(char *s) { + char *a, *b, *rc; + /* allocates memory for a string in */ + /* doublequotes and returns a pointer to it */ + if (strchr(s, '"')) { + a=strchr(s, '"'); + b=strrchr(s, '"'); + if (a<b) { + rc=new char[b-a]; + strncpy(rc, (a+1), b-a-1); + rc[b-a-1]='\0'; + } + else { + rc=new char[1]; + strcpy(rc,"\0"); + } + } + else { + rc=new char[1]; + strcpy(rc,"\0"); + } + return rc; +} + +float prho(int n, float is, int *ifault) { +/* + Algorithm AS 89 Appl. Statist. (1975) Vol.24, No. 3, P377. + + To evaluate the probability of obtaining a value greater than or + equal to is, where is=(n**3-n)*(1-r)/6, r=Spearman's rho and n + must be greater than 1 + + n is the number of data pairs in the correlation + The calculated probability is for one-sided tests. + For two-sided tests multiply the probability value with 2. + + Auxiliary function required: ALNORM = algorithm AS66 +*/ + + float c1=0.2274, c2=0.2531, c3=0.1745, c4=0.0758, c5=0.1033, c6=0.3932; + float c7=0.0879, c8=0.0151, c9=0.0072, c10=0.0831, c11=0.0131, c12=0.00046; + float rc; + double b, x, y, u; + long i, js, ifr, nfac, m, n1, mt, nn, ise; + long l[7]; + + extern float alnorm(float x, int upper); + + /* Test admissibility of arguments and initialize */ + rc=1.0; + *ifault=1; + if (n<=1) return(rc); + *ifault=0; + if (is<=0) return(rc); + rc=0.0; + if (is > n*(n*n-1)/3) return(rc); + js=(long)is; + if (js!= 2*(js/2)) js=js+1; + /* if (n>6) goto mark6; */ + goto mark6; + + /* Exact evaluation of probability */ + nfac=1; + for (i=1; i<=n; i++) { + nfac*=i; + l[i]=i; + } +mark1: + rc=1.0/(float)(nfac); + if (js != n*(n*n-1)/3) return(rc); + ifr=0; + for (m=1; m<=nfac; m++) { + ise=0; + for (i=1; i<=n; i++) { + ise+=(i-l[i])*(i-l[i]); + } + mark2: + if (js <= ise) ifr++; + n1=n; + mark3: + mt=l[1]; + nn=n1-1; + for (i=1; i<=nn; i++) { + l[i]=l[i+1]; + } + mark4: + l[n1]=mt; + if ((l[n1]!=n1) || (n1==2)) goto mark5; + n1--; + if (m!=nfac) goto mark3; + } +mark5: + rc=(float)(ifr)/(float)(nfac); + return(rc); + /* Evaluation by Edgeworth series expansion */ +mark6: + b=1.0/(double)(n); + x=(6.0*((double)(js)-1.0)*b/(1.0/(b*b)-1.0)-1.0)*sqrt(1.0/b-1.0); + y=x*x; + u=x*b*(c1+b*(c2+c3*b)+y*(-c4+b*(c5+c6*b)-y*b*(c7+c8*b-y*(c9-c10*b+y*b*(c11-c12*y))))); + /* Call to algorithm AS 66 */ + rc=u/exp(y/2.0)+alnorm(x,1); + if (rc < 0.0) rc=0.0; + if (rc > 1.0) rc=1.0; + return(rc); +} + + +/* + Algorithm AS66 Applied Statistics (1973) vol22 no.3 + Evaluates the tail area of the standardised normal curve + from x to infinity if upper is .true. or + from minus infinity to x if upper is .false. +*/ +float stdnorm(float x) { + return (1.0/sqrt(2*3.141592653))*exp(-0.5*x*x); +} +float alnorm(float x, int upper) { + float q; + + float stdnorm(float x); + extern float qtrap(float (*func)(float x), float a, float b); + + if (x==0.0) q=0; + else q=qtrap(&stdnorm, 0, x); + + if (upper) q=1.0-q-0.5; + else q+=0.5; + + return(q); +} + + +#define EPS 1.0e-5 +#define JMAX 20 +float qtrap(float (*func)(float), float a, float b) +{ + float trapzd(float (*func)(float), float a, float b, int n); + void nrerror(char error_text[]); + int j; + float s,olds; + + olds = -1.0e30; + for (j=1;j<=JMAX;j++) { + s=trapzd(func,a,b,j); + if (fabs(s-olds) < EPS*fabs(olds)) return s; + olds=s; + } + nrerror("Too many steps in routine qtrap"); + return 0.0; +} +#undef EPS +#undef JMAX + +#define FUNC(x) ((*func)(x)) +float trapzd(float (*func)(float), float a, float b, int n) +{ + float x,tnm,sum,del; + static float s; + int it,j; + + if (n == 1) { + return (s=0.5*(b-a)*(FUNC(a)+FUNC(b))); + } else { + for (it=1,j=1;j<n-1;j++) it <<= 1; + tnm=it; + del=(b-a)/tnm; + x=a+0.5*del; + for (sum=0.0,j=1;j<=it;j++,x+=del) sum += FUNC(x); + s=0.5*(s+(b-a)*sum/tnm); + return s; + } +} +#undef FUNC |