summaryrefslogtreecommitdiff
path: root/support/fastpictex/src/myfunctions.cc
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
committerNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
commite0c6872cf40896c7be36b11dcc744620f10adf1d (patch)
tree60335e10d2f4354b0674ec22d7b53f0f8abee672 /support/fastpictex/src/myfunctions.cc
Initial commit
Diffstat (limited to 'support/fastpictex/src/myfunctions.cc')
-rw-r--r--support/fastpictex/src/myfunctions.cc211
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