1 /* Steven Andrews, 11/10/90.
2 See documentation called Rn_doc.doc.
3 Copyright 2003-2007 by Steven Andrews. This work is distributed under the terms
4 of the Gnu Lesser General Public License (LGPL). */
5
6 #include "Rn.h"
7 #include "random2.h"
8 #include <math.h>
9 #include <string.h>
10 #include "math2.h"
11
12 float detpart(float *a,int n,char s[],int i);
13
14
makeV(float * c,int n,char * s)15 int makeV(float *c,int n,char *s) {
16 int i,as;
17
18 as=0;
19 for(i=0;i<n;i++) {
20 if(sscanf(s,"%f",&c[i])) as++;
21 else c[i]=0;
22 s=strchr(s,' ');
23 if(s) s++; }
24 return as; }
25
26
makeM(float * c,int m,int n,char * s)27 int makeM(float *c,int m,int n,char *s) {
28 return makeV(c,m*n,s); }
29
30
setstdV(float * c,int n,int k)31 float *setstdV(float *c,int n,int k) {
32 int i;
33
34 if(k==0) for(i=0;i<n;i++) c[i]=0;
35 else if(k==1) for(i=0;i<n;i++) c[i]=1;
36 else if(k<0) {
37 for(i=0;i<n;i++) c[i]=0;
38 c[-k]=1; }
39 else if(k==2) for(i=0;i<n;i++) c[i]=i;
40 else if(k==3) for(i=0;i<n;i++) c[i]=randCCF();
41 return c; }
42
43
setstdM(float * c,int m,int n,int k)44 float *setstdM(float *c,int m,int n,int k) {
45 int i,j;
46
47 if(k==0)
48 for(i=0;i<m;i++)
49 for(j=0;j<n;j++) c[n*i+j]=0;
50 else if(k==1)
51 for(i=0;i<m;i++)
52 for(j=0;j<n;j++) c[n*i+j]=1.0*(i==j);
53 else if(k==2)
54 for(i=0;i<m;i++)
55 for(j=0;j<n;j++) c[n*i+j]=1;
56 else if(k==3)
57 for(i=0;i<m;i++)
58 for(j=0;j<n;j++) c[n*i+j]=randCCF();
59 return c; }
60
61
DirCosM(float * c,float theta,float phi,float chi)62 float *DirCosM(float *c,float theta,float phi,float chi) {
63 float cp,ct,cc,sp,st,sc;
64
65 cp=cos(phi);
66 ct=cos(theta);
67 cc=cos(chi);
68 sp=sin(phi);
69 st=sin(theta);
70 sc=sin(chi);
71 c[0]=cp*ct*cc-sp*sc;
72 c[1]=sp*ct*cc+cp*sc;
73 c[2]=-st*cc;
74 c[3]=-cp*ct*sc-sp*cc;
75 c[4]=-sp*ct*sc+cp*cc;
76 c[5]=st*sc;
77 c[6]=cp*st;
78 c[7]=sp*st;
79 c[8]=ct;
80 return c; }
81
82
DirCosMD(double * c,double theta,double phi,double chi)83 double *DirCosMD(double *c,double theta,double phi,double chi) {
84 double cp,ct,cc,sp,st,sc;
85
86 cp=cos(phi);
87 ct=cos(theta);
88 cc=cos(chi);
89 sp=sin(phi);
90 st=sin(theta);
91 sc=sin(chi);
92 c[0]=cp*ct*cc-sp*sc;
93 c[1]=sp*ct*cc+cp*sc;
94 c[2]=-st*cc;
95 c[3]=-cp*ct*sc-sp*cc;
96 c[4]=-sp*ct*sc+cp*cc;
97 c[5]=st*sc;
98 c[6]=cp*st;
99 c[7]=sp*st;
100 c[8]=ct;
101 return c; }
102
103
DirCosM2(float * c,float theta)104 float *DirCosM2(float *c,float theta) {
105 c[0]=c[3]=cos(theta);
106 c[2]=-(c[1]=sin(theta));
107 return c; }
108
109
DirCosM2D(double * c,double theta)110 double *DirCosM2D(double *c,double theta) {
111 c[0]=c[3]=cos(theta);
112 c[2]=-(c[1]=sin(theta));
113 return c; }
114
115
printV(float * a,int n)116 float *printV(float *a,int n) {
117 int i,er=0;
118
119 if(!a) return NULL;
120 if(n) if(printf("%f",a[0])<0) er=1;
121 for(i=1;i<n;i++) if(printf(" %f",a[i])<0) er=1;
122 if(printf("\n")<0) er=1;
123 return er?NULL:a; }
124
125
printVD(double * a,int n)126 double *printVD(double *a,int n) {
127 int i,er=0;
128
129 if(!a) return NULL;
130 if(n) if(printf("%g",a[0])<0) er=1;
131 for(i=1;i<n;i++) if(printf(" %g",a[i])<0) er=1;
132 if(printf("\n")<0) er=1;
133 return er?NULL:a; }
134
135
fprintV(FILE * stream,float * a,int n)136 float *fprintV(FILE *stream,float *a,int n) {
137 int i,er=0;
138
139 if(!a) return NULL;
140 for(i=0;i<n;i++) if(fprintf(stream,"%f ",a[i])<0) er=1;
141 if(fprintf(stream,"\n")<0) er=1;
142 return er?NULL:a; }
143
144
fprintVD(FILE * stream,double * a,int n)145 double *fprintVD(FILE *stream,double *a,int n) {
146 int i,er=0;
147
148 if(!a) return NULL;
149 for(i=0;i<n;i++) if(fprintf(stream,"%g ",a[i])<0) er=1;
150 if(fprintf(stream,"\n")<0) er=1;
151 return er?NULL:a; }
152
153
printM(float * a,int m,int n,char s[])154 float *printM(float *a,int m,int n,char s[]) {
155 int i,j,er=0;
156 char dft[]="%f ";
157
158 if(!a) return NULL;
159 if(!s) s=dft;
160 else if(s[0]=='\0') s=dft;
161 for(i=0;i<m;i++) {
162 for(j=0;j<n;j++)
163 if(printf(s,a[n*i+j])<0) er=1;
164 if(printf("\n")<0) er=1; }
165 return er?NULL:a; }
166
167
sprintM(float * a,int m,int n,char s[],char * t,int tn)168 float *sprintM(float *a,int m,int n,char s[],char *t,int tn) {
169 int i,j,ti;
170 char dft[]="%f ";
171 char s2[255];
172
173 if(!a) return NULL;
174 if(!s) s=dft;
175 else if(s[0]=='\0') s=dft;
176 ti=0;
177 for(i=0;i<m;i++) {
178 for(j=0;j<n;j++)
179 if(snprintf(s2,255,s,a[n*i+j])<tn-ti-1)
180 ti+=sprintf(t+ti,s,a[n*i+j]);
181 if(tn-ti>1) ti+=sprintf(t+ti,"\n"); }
182 return a; }
183
184
minV(float * a,int n)185 float minV(float *a,int n) {
186 float min;
187 int i;
188
189 min=a[0];
190 for(i=1;i<n;i++) if(a[i]<min) min=a[i];
191 return min; }
192
193
maxV(float * a,int n)194 float maxV(float *a,int n) {
195 float max;
196 int i;
197
198 max=a[0];
199 for(i=1;i<n;i++) if(a[i]>max) max=a[i];
200 return max; }
201
202
maxVD(double * a,int n,int * indx)203 double maxVD(double *a,int n,int *indx) {
204 double max;
205 int i,imax;
206
207 max=a[imax=0];
208 for(i=1;i<n;i++) if(a[i]>max) max=a[imax=i];
209 if(indx) *indx=imax;
210 return max; }
211
212
minVD(double * a,int n,int * indx)213 double minVD(double *a,int n,int *indx) {
214 double min;
215 int i,imin;
216
217 min=a[imin=0];
218 for(i=1;i<n;i++) if(a[i]<min) min=a[imin=i];
219 if(indx) *indx=imin;
220 return min; }
221
222
equalV(float * a,float * b,int n)223 int equalV(float *a,float *b,int n) {
224 int i;
225
226 for(i=0;i<n;i++) if(a[i]!=b[i]) return 0;
227 return 1; }
228
229
isevenspV(float * a,int n,float tol)230 int isevenspV(float *a,int n,float tol) {
231 int i;
232 float sp;
233
234 if(n<2) return 0;
235 sp=(a[n-1]-a[0])/(n-1);
236 tol*=fabs(sp);
237 for(i=1;i<n;i++)
238 if(fabs(a[i]-a[i-1]-sp)>tol) return 0;
239 return 1; }
240
241
detpart(float * a,int n,char s[],int i)242 float detpart(float *a,int n,char s[],int i) {
243 int j,p;
244 float sum;
245
246 if(i==n-1) {
247 for(j=0;s[j];j++)
248 ;
249 return a[n*i+j]; }
250 else {
251 sum=0;
252 p=1;
253 for(j=0;j<n;j++)
254 if(!s[j]) {
255 s[j]=1;
256 sum+=p*a[n*i+j]*detpart(a,n,s,i+1);
257 p=-p;
258 s[j]=0; }
259 return sum; }}
260
261
detM(float * a,int n)262 float detM(float *a,int n) {
263 int j;
264 char *s;
265 float det;
266
267 s=(char *) calloc(n,sizeof(char));
268 if(!s) return 0;
269 for(j=0;j<n;j++) s[j]=0;
270 det=detpart(a,n,s,0);
271 free(s);
272 return det; }
273
274
traceM(float * a,int n)275 float traceM(float *a,int n) {
276 int j;
277 float ans;
278
279 ans=0;
280 for(j=0;j<n;j++) ans+=a[j*n+j];
281 return ans; }
282
283
traceMD(double * a,int n)284 double traceMD(double *a,int n) {
285 int j;
286 double ans;
287
288 ans=0;
289 for(j=0;j<n;j++) ans+=a[j*n+j];
290 return ans; }
291
292
issymmetricMD(double * a,int n)293 int issymmetricMD(double *a,int n) {
294 int i,j,ans;
295
296 ans=1;
297 for(i=1;i<n&&ans;i++)
298 for(j=0;j<i&&ans;j++)
299 if(a[n*i+j]!=a[n*j+i]) ans=0;
300 return ans; }
301
302
columnM(float * a,float * c,int m,int n,int j)303 float *columnM(float *a,float *c,int m,int n,int j) {
304 int i;
305
306 for(i=0;i<m;i++)
307 c[i]=a[n*i+j];
308 return c; }
309
310
copyV(float * a,float * c,int n)311 float *copyV(float *a,float *c,int n) {
312 int i;
313
314 for(i=0;i<n;i++) c[i]=a[i];
315 return c; }
316
317
copyVD(double * a,double * c,int n)318 double *copyVD(double *a,double *c,int n) {
319 int i;
320
321 for(i=0;i<n;i++) c[i]=a[i];
322 return c; }
323
324
copyM(float * a,float * c,int m,int n)325 float *copyM(float *a,float *c,int m,int n) {
326 int i;
327
328 n*=m;
329 for(i=0;i<n;i++) c[i]=a[i];
330 return c; }
331
332
leftrotV(float * a,float * c,int n,int k)333 float *leftrotV(float *a,float *c,int n,int k) {
334 int i,i2,gcd;
335 float temp;
336
337 if(k<0) k+=((-k)/n+1)*n;
338 else k-=(k/n)*n;
339 if(k==0) return copyV(a,c,n);
340 gcd=gcomdiv(n,k);
341 for(i2=0;i2<gcd;i2++) {
342 temp=a[i2];
343 for(i=i2;(i+k)%n!=i2;i=(i+k)%n)
344 c[i]=a[(i+k)%n];
345 c[i]=temp; }
346 return c; }
347
348
sumV(float ax,float * a,float bx,float * b,float * c,int n)349 float *sumV(float ax,float *a,float bx,float *b,float *c,int n) {
350 int i;
351
352 for(i=0;i<n;i++) c[i]=ax*a[i]+bx*b[i];
353 return c; }
354
355
sumVD(double ax,double * a,double bx,double * b,double * c,int n)356 double *sumVD(double ax,double *a,double bx,double *b,double *c,int n) {
357 int i;
358
359 for(i=0;i<n;i++) c[i]=ax*a[i]+bx*b[i];
360 return c; }
361
362
sumM(float ax,float * a,float bx,float * b,float * c,int m,int n)363 float *sumM(float ax,float *a,float bx,float *b,float *c,int m,int n) {
364 int i;
365
366 n*=m;
367 for(i=0;i<n;i++) c[i]=ax*a[i]+bx*b[i];
368 return c; }
369
370
addKV(float k,float * a,float * c,int n)371 float *addKV(float k,float *a,float *c,int n) {
372 int i;
373
374 for(i=0;i<n;i++) c[i]=k+a[i];
375 return c; }
376
377
multKV(float k,float * a,float * c,int n)378 float *multKV(float k,float *a,float *c,int n) {
379 int i;
380
381 for(i=0;i<n;i++) c[i]=k*a[i];
382 return c; }
383
384
divKV(float k,float * a,float * c,int n)385 float *divKV(float k,float *a,float *c,int n) {
386 int i;
387
388 for(i=0;i<n;i++) c[i]=k/a[i];
389 return c; }
390
391
multV(float * a,float * b,float * c,int n)392 float *multV(float *a,float *b,float *c,int n) {
393 int i;
394
395 for(i=0;i<n;i++) c[i]=a[i]*b[i];
396 return c; }
397
398
multM(float * a,float * b,float * c,int m,int n)399 float *multM(float *a,float *b,float *c,int m,int n) {
400 int i;
401
402 n*=m;
403 for(i=0;i<n;i++) c[i]=a[i]*b[i];
404 return c; }
405
406
divV(float * a,float * b,float * c,int n)407 float *divV(float *a,float *b,float *c,int n) {
408 int i;
409
410 for(i=0;i<n;i++) c[i]=a[i]/b[i];
411 return c; }
412
413
divM(float * a,float * b,float * c,int m,int n)414 float *divM(float *a,float *b,float *c,int m,int n) {
415 int i;
416
417 n*=m;
418 for(i=0;i<n;i++) c[i]=a[i]/b[i];
419 return c; }
420
421
transM(float * a,float * c,int m,int n)422 float *transM(float *a,float *c,int m,int n) {
423 int i,j;
424
425 for(i=0;i<m;i++) for(j=0;j<n;j++) c[m*j+i]=a[n*i+j];
426 return c; }
427
428
dotVV(float * a,float * b,int n)429 float dotVV(float *a,float *b,int n) {
430 float x=0;
431 int i;
432
433 for(i=0;i<n;i++) x+=a[i]*b[i];
434 return x; }
435
436
dotVVD(double * a,double * b,int n)437 double dotVVD(double *a,double *b,int n) {
438 double x=0;
439 int i;
440
441 for(i=0;i<n;i++) x+=a[i]*b[i];
442 return x; }
443
444
crossVV(float * a,float * b,float * c)445 void crossVV(float *a,float *b,float *c) {
446 c[0]=a[2]*b[3]-a[3]*b[2];
447 c[1]=a[3]*b[1]-a[1]*b[3];
448 c[2]=a[1]*b[2]-a[2]*b[1];
449 return; }
450
451
distanceVV(float * a,float * b,int n)452 float distanceVV(float *a,float *b,int n) {
453 int i;
454 double x=0;
455
456 for(i=0;i<n;i++) x+=(a[i]-b[i])*(a[i]-b[i]);
457 return sqrt(x); }
458
459
distanceVVD(double * a,double * b,int n)460 double distanceVVD(double *a,double *b,int n) {
461 int i;
462 double x=0;
463
464 for(i=0;i<n;i++) x+=(a[i]-b[i])*(a[i]-b[i]);
465 return sqrt(x); }
466
467
normalizeV(float * a,int n)468 float normalizeV(float *a,int n) {
469 int i;
470 float x=0;
471
472 for(i=0;i<n;i++) x+=a[i]*a[i];
473 if(!x) return 0;
474 x=sqrt(x);
475 for(i=0;i<n;i++) a[i]/=x;
476 return x; }
477
478
normalizeVD(double * a,int n)479 float normalizeVD(double *a,int n) {
480 int i;
481 double x=0;
482
483 for(i=0;i<n;i++) x+=a[i]*a[i];
484 if(!x) return 0;
485 x=sqrt(x);
486 for(i=0;i<n;i++) a[i]/=x;
487 return x; }
488
489
averageV(float * a,int n,int p)490 float averageV(float *a,int n,int p) {
491 int i;
492 double x=0;
493
494 if(p==1)
495 for(i=0;i<n;i++) x+=a[i];
496 else if(p==2)
497 for(i=0;i<n;i++) x+=a[i]*a[i];
498 else if(p==-1)
499 for(i=0;i<n;i++) x+=1.0/a[i];
500 else if(p==0) x=n;
501 else for(i=0;i<n;i++) x+=pow(a[i],p);
502 return x/n; }
503
504
dotVM(float * a,float * b,float * c,int m,int n)505 float *dotVM(float *a,float *b,float *c,int m,int n) {
506 int i,j;
507
508 for(j=0;j<n;j++) {
509 c[j]=0;
510 for(i=0;i<m;i++) c[j]+=a[i]*b[n*i+j]; }
511 return c; }
512
513
dotMV(float * a,float * b,float * c,int m,int n)514 float *dotMV(float *a,float *b,float *c,int m,int n) {
515 int i,j;
516
517 for(i=0;i<m;i++) {
518 c[i]=0;
519 for(j=0;j<n;j++) c[i]+=a[n*i+j]*b[j]; }
520 return c; }
521
522
dotMVD(double * a,double * b,double * c,int m,int n)523 double *dotMVD(double *a,double *b,double *c,int m,int n) {
524 int i,j;
525
526 for(i=0;i<m;i++) {
527 c[i]=0;
528 for(j=0;j<n;j++) c[i]+=a[n*i+j]*b[j]; }
529 return c; }
530
531
dotMM(float * a,float * b,float * c,int ma,int na,int nb)532 float *dotMM(float *a,float *b,float *c,int ma,int na,int nb) {
533 int i,j,k;
534
535 for(i=0;i<ma;i++)
536 for(k=0;k<nb;k++) {
537 c[nb*i+k]=0;
538 for(j=0;j<na;j++) c[nb*i+k]+=a[na*i+j]*b[nb*j+k]; }
539 return c; }
540
541
dotMMD(double * a,double * b,double * c,int ma,int na,int nb)542 double *dotMMD(double *a,double *b,double *c,int ma,int na,int nb) {
543 int i,j,k;
544
545 for(i=0;i<ma;i++)
546 for(k=0;k<nb;k++) {
547 c[nb*i+k]=0;
548 for(j=0;j<na;j++) c[nb*i+k]+=a[na*i+j]*b[nb*j+k]; }
549 return c; }
550
551
minorM(float * a,int n,char * row,char * col)552 float minorM(float *a,int n,char *row,char *col) {
553 int i,j,sgn;
554 float sum;
555
556 for(i=0;row[i]&&i<n;i++)
557 ;
558 if(i==n) return 1;
559 row[i]=1;
560 sum=0;
561 sgn=1;
562 for(j=0;j<n;j++)
563 if(!col[j]) {
564 col[j]=1;
565 sum+=sgn*a[n*i+j]*minorM(a,n,row,col);
566 sgn=-sgn;
567 col[j]=0; }
568 row[i]=0;
569 return sum; }
570
571
invM(float * a,float * c,int n)572 float invM(float *a,float *c,int n) {
573 int i,j;
574 char *row,*col;
575 float det;
576
577 det=detM(a,n);
578 if(det) {
579 row=(char *) calloc(n,sizeof(char));
580 if(!row) return 0;
581 col=(char *) calloc(n,sizeof(char));
582 if(!col) return 0;
583 for(j=0;j<n;j++) row[j]=col[j]=0;
584 for(i=0;i<n;i++) {
585 row[i]=1;
586 for(j=0;j<n;j++) {
587 col[j]=1;
588 c[n*j+i]=minus1to(i+j)*minorM(a,n,row,col)/det;
589 col[j]=0; }
590 row[i]=0; }}
591 return det; }
592
593
deriv1V(float * a,float * c,int n)594 float *deriv1V(float *a,float *c,int n) {
595 int i;
596
597 if(n==1) c[0]=0;
598 else if(n==2) c[0]=c[1]=a[1]-a[0];
599 else {
600 c[0]=-1.5*a[0]+2.0*a[1]-0.5*a[2];
601 for(i=1;i<n-1;i++)
602 c[i]=(a[i+1]-a[i-1])/2.0;
603 c[n-1]=0.5*a[n-3]-2.0*a[n-2]+1.5*a[n-1]; }
604 return c; }
605
606
deriv2V(float * a,float * c,int n)607 float *deriv2V(float *a,float *c,int n) {
608 int i;
609
610 if(n==1) c[0]=0;
611 else if(n==2) c[0]=c[1]=0;
612 else {
613 c[0]=a[0]+a[2]-2.0*a[1];
614 for(i=1;i<n-1;i++)
615 c[i]=a[i-1]+a[i+1]-2.0*a[i];
616 c[n-1]=a[n-3]+a[n-1]-2.0*a[n-2]; }
617 return c; }
618
619
integV(float * a,float * c,int n)620 float *integV(float *a,float *c,int n) {
621 int i;
622
623 c[0]=a[0]/2;
624 for(i=1;i<n;i++)
625 c[i]=c[i-1]+(a[i-1]+a[i])/2;
626 return c; }
627
628
smoothV(float * a,float * c,int n,int k)629 float *smoothV(float *a,float *c,int n,int k) {
630 int i,j;
631 float *smo,norm;
632
633 if(k<0) return NULL;
634 smo=allocV(2*k+1);
635 if(!smo) return NULL;
636 smo+=k;
637 for(j=-k;j<=k;j++) smo[j]=choose(2*k,j+k);
638 setstdV(c,n,0);
639
640 for(i=0;i<n;i++) {
641 norm=0;
642 for(j=-k;j<=k;j++)
643 if(i+j>=0&&i+j<n) {
644 norm+=smo[j];
645 c[i]+=smo[j]*a[i+j]; }
646 c[i]/=norm; }
647 freeV(smo-k);
648 return c; }
649
650
convolveV(float * a,float * b,float * c,int na,int nb,int nc,int bz,float left,float right)651 float *convolveV(float *a,float *b,float *c,int na,int nb,int nc,int bz,float left,float right) {
652 int i,j,blo,bhi;
653 double sum;
654
655 blo=-bz;
656 bhi=nb-bz;
657 b+=bz;
658 for(i=0;i<nc;i++) {
659 sum=0;
660 for(j=blo;i-j>=na&&j<bhi;j++) sum+=right*b[j];
661 for(;i-j>=0&&j<bhi;j++) sum+=a[i-j]*b[j];
662 for(;j<bhi;j++) sum+=left*b[j];
663 c[i]=sum; }
664 return c; }
665
666
correlateV(float * a,float * b,float * c,int na,int nb,int nc,int bz,float left,float right)667 float *correlateV(float *a,float *b,float *c,int na,int nb,int nc,int bz,float left,float right) {
668 int i,j,blo,bhi;
669 double sum;
670
671 blo=-bz;
672 bhi=nb-bz;
673 b+=bz;
674 for(i=0;i<nc;i++) {
675 sum=0;
676 for(j=blo;i+j<0&&j<bhi;j++) sum+=left*b[j];
677 for(;i+j<na&&j<bhi;j++) sum+=a[i+j]*b[j];
678 for(;j<bhi;j++) sum+=right*b[j];
679 c[i]=sum; }
680 return c; }
681
682
histogramV(float * a,float * h,float lo,float hi,int na,int nh)683 int histogramV(float *a,float *h,float lo,float hi,int na,int nh) {
684 int asgn,i,j;
685 float delta;
686
687 delta=(hi-lo)/(nh-1);
688 asgn=na;
689 for(j=0;j<nh;j++) h[j]=0;
690 for(i=0;i<na;i++) {
691 j=(int)floor((a[i]-lo)/delta)+1;
692 if(j<=0) h[0]+=1.0;
693 else if(j<nh) h[j]+=1.0;
694 else asgn--; }
695 return asgn; }
696
697
histogramVdbl(double * a,double * h,double lo,double hi,int na,int nh)698 int histogramVdbl(double *a,double *h,double lo,double hi,int na,int nh) {
699 int asgn,i,j;
700 double delta;
701
702 delta=(hi-lo)/(nh-1);
703 asgn=na;
704 for(j=0;j<nh;j++) h[j]=0;
705 for(i=0;i<na;i++) {
706 j=(int)floor((a[i]-lo)/delta)+1;
707 if(j<=0) h[0]+=1.0;
708 else if(j<nh) h[j]+=1.0;
709 else asgn--; }
710 return asgn; }
711
712
713
714
715