Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "LatticeReduction.h"
23 : #include "Exception.h"
24 : #include <cstdio>
25 : #include <cmath>
26 : #include <iomanip>
27 :
28 : namespace PLMD {
29 :
30 : const double epsilon=1e-14;
31 :
32 48700 : void LatticeReduction::sort(Vector v[3]) {
33 : const double onePlusEpsilon=(1.0+epsilon);
34 : double m[3];
35 48700 : m[0]=modulo2(v[0]);
36 48700 : m[1]=modulo2(v[1]);
37 48700 : m[2]=modulo2(v[2]);
38 194800 : for(int i=0; i<3; i++) {
39 292200 : for(int j=i+1; j<3; j++) {
40 146100 : if(m[i]>m[j]) {
41 12265 : std::swap(v[i],v[j]);
42 : std::swap(m[i],m[j]);
43 : }
44 : }
45 : }
46 48700 : if(m[0]>m[1]*onePlusEpsilon) {
47 0 : plumed_error() << std::scientific << std::setprecision(17)
48 : << "v[0] " << v[0] << " " << "m[0] " << m[0] <<"\n"
49 : << "v[1] " << v[1] << " " << "m[1] " << m[1] <<"\n"
50 0 : << "v[2] " << v[2] << " " << "m[2] " << m[2] <<"\n";
51 : }
52 48700 : if(m[1]>m[2]*onePlusEpsilon) {
53 0 : plumed_error() << std::scientific << std::setprecision(17)
54 : << "v[0] " << v[0] << " " << "m[0] " << m[0] <<"\n"
55 : << "v[1] " << v[1] << " " << "m[1] " << m[1] <<"\n"
56 0 : << "v[2] " << v[2] << " " << "m[2] " << m[2] <<"\n";
57 : }
58 48700 : }
59 :
60 27811 : void LatticeReduction::reduce(Vector&a,Vector&b) {
61 : const double onePlusEpsilon=(1.0+epsilon);
62 27811 : double ma=modulo2(a);
63 27811 : double mb=modulo2(b);
64 : unsigned counter=0;
65 : while(true) {
66 28395 : if(mb>ma) {
67 : std::swap(a,b);
68 : std::swap(ma,mb);
69 : }
70 28395 : a-=b*floor(dotProduct(a,b)/mb+0.5);
71 28395 : ma=modulo2(a);
72 28395 : if(mb<=ma*onePlusEpsilon) {
73 : break;
74 : }
75 584 : counter++;
76 584 : if(counter%100==0) { // only test rarely since this might be expensive
77 0 : plumed_assert(!std::isnan(ma));
78 0 : plumed_assert(!std::isnan(mb));
79 : }
80 584 : if(counter%10000==0) {
81 0 : std::fprintf(stderr,"WARNING: LatticeReduction::reduce stuck after %u iterations\n",counter);
82 : }
83 : }
84 :
85 : std::swap(a,b);
86 27811 : }
87 :
88 1 : void LatticeReduction::reduce2(Vector&a,Vector&b,Vector&c) {
89 4 : Vector v[3];
90 1 : v[0]=a;
91 1 : v[1]=b;
92 1 : v[2]=c;
93 : int iter=0;
94 : int ok=0;
95 12 : while(ok<3) {
96 : int i,j;
97 11 : if(iter%3==0) {
98 : i=0;
99 : j=1;
100 7 : } else if(iter%3==1) {
101 : i=0;
102 : j=2;
103 : } else {
104 : i=1;
105 : j=2;
106 : }
107 11 : if(isReduced(v[i],v[j])) {
108 4 : ok++;
109 : } else {
110 7 : reduce(v[i],v[j]);
111 : ok=1;
112 : }
113 11 : iter++;
114 : }
115 1 : a=v[0];
116 1 : b=v[1];
117 1 : c=v[2];
118 1 : }
119 :
120 11 : bool LatticeReduction::isReduced(const Vector&a,const Vector&b) {
121 : const int cut=5;
122 84 : for(int i=-cut; i<=cut; i++) {
123 79 : if(modulo2(b+i*a)<modulo2(b)) {
124 : return false;
125 : }
126 : }
127 5 : return modulo2(a)<=modulo2(b) && 2.0*dotProduct(a,b)<=modulo2(a);
128 : }
129 :
130 0 : void LatticeReduction::reduce2(Tensor&t) {
131 0 : Vector a=t.getRow(0);
132 0 : Vector b=t.getRow(1);
133 0 : Vector c=t.getRow(2);
134 0 : reduce2(a,b,c);
135 0 : t.setRow(0,a);
136 0 : t.setRow(1,b);
137 0 : t.setRow(2,c);
138 0 : }
139 :
140 20894 : void LatticeReduction::reduce(Tensor&t) {
141 20894 : reduceFast(t);
142 20894 : }
143 :
144 20895 : void LatticeReduction::reduceFast(Tensor&t) {
145 : const double onePlusEpsilon=(1.0+epsilon);
146 83580 : Vector v[3];
147 20895 : v[0]=t.getRow(0);
148 20895 : v[1]=t.getRow(1);
149 20895 : v[2]=t.getRow(2);
150 : unsigned counter=0;
151 : while(true) {
152 27804 : sort(v);
153 27804 : reduce(v[0],v[1]);
154 27804 : double b11=modulo2(v[0]);
155 27804 : double b22=modulo2(v[1]);
156 27804 : double b12=dotProduct(v[0],v[1]);
157 27804 : double b13=dotProduct(v[0],v[2]);
158 27804 : double b23=dotProduct(v[1],v[2]);
159 27804 : double z=b11*b22-b12*b12;
160 27804 : double y2=-(b11*b23-b12*b13)/z;
161 27804 : double y1=-(b22*b13-b12*b23)/z;
162 27804 : int x1min=floor(y1);
163 27804 : int x1max=x1min+1;
164 27804 : int x2min=floor(y2);
165 27804 : int x2max=x2min+1;
166 : bool first=true;
167 : double mbest,mtrial;
168 27804 : Vector trial,best;
169 83412 : for(int x1=x1min; x1<=x1max; x1++)
170 166824 : for(int x2=x2min; x2<=x2max; x2++) {
171 111216 : trial=v[2]+x2*v[1]+x1*v[0];
172 111216 : mtrial=modulo2(trial);
173 111216 : if(first || mtrial<mbest) {
174 : mbest=mtrial;
175 35415 : best=trial;
176 : first=false;
177 : }
178 : }
179 27804 : if(modulo2(best)*onePlusEpsilon>=modulo2(v[2])) {
180 : break;
181 : }
182 6909 : counter++;
183 6909 : if(counter%10000==0) {
184 0 : std::fprintf(stderr,"WARNING: LatticeReduction::reduceFast stuck after %u iterations\n",counter);
185 : }
186 6909 : v[2]=best;
187 6909 : }
188 20895 : sort(v);
189 20895 : t.setRow(0,v[0]);
190 20895 : t.setRow(1,v[1]);
191 20895 : t.setRow(2,v[2]);
192 20895 : }
193 :
194 :
195 1 : void LatticeReduction::reduceSlow(Tensor&t) {
196 4 : Vector v[3];
197 1 : v[0]=t.getRow(0);
198 1 : v[1]=t.getRow(1);
199 1 : v[2]=t.getRow(2);
200 1 : reduce2(v[0],v[1],v[2]);
201 1 : double e01=dotProduct(v[0],v[1]);
202 1 : double e02=dotProduct(v[0],v[2]);
203 1 : double e12=dotProduct(v[1],v[2]);
204 1 : if(e01*e02*e12<0) {
205 : int eps01=0;
206 0 : if(e01>0.0) {
207 : eps01=1;
208 0 : } else if(e01<0.0) {
209 : eps01=-1;
210 : }
211 : int eps02=0;
212 0 : if(e02>0.0) {
213 : eps02=1;
214 0 : } else if(e02<0.0) {
215 : eps02=-1;
216 : }
217 0 : Vector n=v[0]-eps01*v[1]-eps02*v[2];
218 : int i=0;
219 0 : double mx=modulo2(v[i]);
220 0 : for(int j=1; j<3; j++) {
221 0 : double f=modulo2(v[j]);
222 0 : if(f>mx) {
223 : i=j;
224 : mx=f;
225 : }
226 : }
227 0 : if(modulo2(n)<mx) {
228 0 : v[i]=n;
229 : }
230 : }
231 1 : sort(v);
232 1 : t.setRow(0,v[0]);
233 1 : t.setRow(1,v[1]);
234 1 : t.setRow(2,v[2]);
235 1 : }
236 :
237 0 : bool LatticeReduction::isReduced2(const Vector&a,const Vector&b,const Vector &c) {
238 0 : return isReduced(a,b) && isReduced(a,b) && isReduced(b,c);
239 : }
240 :
241 2 : bool LatticeReduction::isReduced(const Tensor&t) {
242 8 : Vector v[3];
243 : double m[3];
244 2 : v[0]=t.getRow(0);
245 2 : v[1]=t.getRow(1);
246 2 : v[2]=t.getRow(2);
247 8 : for(int i=0; i<3; i++) {
248 6 : m[i]=modulo2(v[i]);
249 : }
250 2 : if(!((m[0]<=m[1]) && m[1]<=m[2])) {
251 : return false;
252 : }
253 : const int cut=5;
254 13 : for(int i=-cut; i<=cut; i++) {
255 12 : double mm=modulo2(v[1]+i*v[0]);
256 12 : if(mm<m[1]) {
257 : return false;
258 : }
259 142 : for(int j=-cut; j<=cut; j++) {
260 131 : double mx=modulo2(v[2]+i*v[1]+j*v[0]);
261 131 : if(mx<m[2]) {
262 : return false;
263 : }
264 : }
265 : }
266 : return true;
267 : }
268 :
269 : }
|
