1 // This file is part of the FXT library.
2 // Copyright (C) 2010 Joerg Arndt
3 // License: GNU General Public License version 3 or later,
4 // see the file COPYING.txt in the main directory.
5
6 #include <cmath> // fabs()
7
8 #include "fxttypes.h"
9 #include "aux0/sumdiff.h"
10 #include "aux1/arith1.h"
11 //#include "jjassert.h"
12
13 #include "aux1/copy.h"
14
15
16 int
symm_q(const double * f,ulong n,double eps)17 symm_q(const double *f, ulong n, double eps/*=1e-5*/)
18 //
19 // check symmetry k <--> n-1-k (k=0..n-1)
20 // return:
21 // 0: no symmetry
22 // 1: symmetry
23 // -1: neg. symm.
24 //
25 {
26 ulong k, nh = n/2;
27 for (k=0; k<nh; ++k) if ( fabs(f[k]-f[n-1-k]) > eps ) break;
28 if ( k==nh ) return +1;
29 else
30 {
31 for (k=0; k<nh; ++k) if ( fabs(f[k]+f[n-1-k]) > eps ) break;
32 if ( k==nh ) return -1;
33 }
34 return 0;
35 }
36 // -------------------------
37
38
39 int
symm_0_q(const double * f,ulong n,double eps)40 symm_0_q(const double *f, ulong n, double eps/*=1e-5*/)
41 // check symmetry k <--> n-k (k=1..n/2-1)
42 {
43 ulong k, nh = n/2;
44 for (k=1; k<nh; ++k) if ( fabs(f[k]-f[n-k]) > eps ) break;
45
46 if ( k==nh ) return +1;
47
48 for (k=1; k<nh; ++k) if ( fabs(f[k]+f[n-k]) > eps ) break;
49
50 if ( k==nh )
51 {
52 if ( (fabs(f[nh])<eps) && (fabs(f[0])<eps) ) return -1;
53 }
54
55 return 0;
56 }
57 // -------------------------
58
59
60 int
any_symm_q(const double * f,ulong n,double eps)61 any_symm_q(const double *f, ulong n, double eps/*=1e-5*/)
62 {
63 return symm_q(f, n, eps) + 2 * symm_0_q(f, n, eps);
64 }
65 // -------------------------
66
67
68 void
symm_0(double * f,ulong n,int s=+1)69 symm_0(double *f, ulong n, int s=+1)
70 // make symmetric k <--> n-k (k=1..n/2-1)
71 {
72 ulong nh = n/2;
73
74 for (ulong k=1, j=n-1; k<j; ++k, --j) sumdiff(f[k],f[j]);
75
76 if ( s>0 )
77 {
78 copy_reverse(f+1, f+nh+1, nh-1);
79 f[0] *= 2;
80 if ( 0==(n&1) ) f[nh] *= 2; // only for n even
81 }
82 else
83 {
84 copy_reverse(f+nh+1, f+1, nh-1);
85 negate(f+1, nh-1);
86 if ( 0==(n&1) ) f[nh] = 0.0; // only for n even
87 f[0] = 0.0;
88 }
89
90 // jjassert( s == symm_0_q(f, n) );
91 }
92 // -------------------------
93
94
95 void
symm(double * f,ulong n,int s=+1)96 symm(double *f, ulong n, int s=+1)
97 // make symmetric k <--> n-1-k (k=1..n/2-1)
98 {
99 ulong nh = n/2;
100
101 for (ulong k=0, j=n-1; k<j; ++k, --j) sumdiff(f[k],f[j]);
102
103 if ( s>0 )
104 {
105 copy_reverse(f, f+nh, nh);
106 if ( n & 1 ) f[nh] *= 2.0; // only for n odd
107 }
108 else
109 {
110 copy_reverse(f+nh, f, nh);
111 negate(f, nh);
112 if ( n & 1 ) f[nh] = 0.0; // only for n odd
113 }
114
115 // jjassert( s == symm_q(f, n) );
116 }
117 // -------------------------
118