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 :
27 : namespace PLMD {
28 :
29 : const double epsilon=1e-14;
30 :
31 55460 : void LatticeReduction::sort(Vector v[3]) {
32 : const double onePlusEpsilon=(1.0+epsilon);
33 : double m[3];
34 55460 : m[0]=modulo2(v[0]);
35 55460 : m[1]=modulo2(v[1]);
36 55460 : m[2]=modulo2(v[2]);
37 221840 : for(int i=0; i<3; i++)
38 332760 : for(int j=i+1; j<3; j++)
39 166380 : if(m[i]>m[j]) {
40 12236 : std::swap(v[i],v[j]);
41 : std::swap(m[i],m[j]);
42 : }
43 55460 : plumed_assert(m[0]<=m[1]*onePlusEpsilon);
44 55460 : plumed_assert(m[1]<=m[2]*onePlusEpsilon);
45 55460 : }
46 :
47 35223 : void LatticeReduction::reduce(Vector&a,Vector&b) {
48 : const double onePlusEpsilon=(1.0+epsilon);
49 35223 : double ma=modulo2(a);
50 35223 : double mb=modulo2(b);
51 : unsigned counter=0;
52 : while(true) {
53 35802 : if(mb>ma) {
54 : std::swap(a,b);
55 : std::swap(ma,mb);
56 : }
57 35802 : a-=b*floor(dotProduct(a,b)/mb+0.5);
58 35802 : ma=modulo2(a);
59 35802 : if(mb<=ma*onePlusEpsilon) {
60 : break;
61 : }
62 579 : counter++;
63 579 : if(counter%100==0) { // only test rarely since this might be expensive
64 0 : plumed_assert(!std::isnan(ma));
65 0 : plumed_assert(!std::isnan(mb));
66 : }
67 579 : if(counter%10000==0) {
68 0 : std::fprintf(stderr,"WARNING: LatticeReduction::reduce stuck after %u iterations\n",counter);
69 : }
70 : }
71 :
72 : std::swap(a,b);
73 35223 : }
74 :
75 1 : void LatticeReduction::reduce2(Vector&a,Vector&b,Vector&c) {
76 4 : Vector v[3];
77 1 : v[0]=a;
78 1 : v[1]=b;
79 1 : v[2]=c;
80 : int iter=0;
81 : int ok=0;
82 12 : while(ok<3) {
83 : int i,j;
84 11 : if(iter%3==0) {
85 : i=0;
86 : j=1;
87 7 : } else if(iter%3==1) {
88 : i=0;
89 : j=2;
90 : } else {
91 : i=1;
92 : j=2;
93 : }
94 11 : if(isReduced(v[i],v[j])) {
95 4 : ok++;
96 : } else {
97 7 : reduce(v[i],v[j]);
98 : ok=1;
99 : }
100 11 : iter++;
101 : }
102 1 : a=v[0];
103 1 : b=v[1];
104 1 : c=v[2];
105 1 : }
106 :
107 11 : bool LatticeReduction::isReduced(const Vector&a,const Vector&b) {
108 : const int cut=5;
109 84 : for(int i=-cut; i<=cut; i++) {
110 79 : if(modulo2(b+i*a)<modulo2(b)) {
111 : return false;
112 : }
113 : }
114 5 : return modulo2(a)<=modulo2(b) && 2.0*dotProduct(a,b)<=modulo2(a);
115 : }
116 :
117 0 : void LatticeReduction::reduce2(Tensor&t) {
118 0 : Vector a=t.getRow(0);
119 0 : Vector b=t.getRow(1);
120 0 : Vector c=t.getRow(2);
121 0 : reduce2(a,b,c);
122 0 : t.setRow(0,a);
123 0 : t.setRow(1,b);
124 0 : t.setRow(2,c);
125 0 : }
126 :
127 20242 : void LatticeReduction::reduce(Tensor&t) {
128 20242 : reduceFast(t);
129 20242 : }
130 :
131 20243 : void LatticeReduction::reduceFast(Tensor&t) {
132 : const double onePlusEpsilon=(1.0+epsilon);
133 80972 : Vector v[3];
134 20243 : v[0]=t.getRow(0);
135 20243 : v[1]=t.getRow(1);
136 20243 : v[2]=t.getRow(2);
137 : unsigned counter=0;
138 : while(true) {
139 35216 : sort(v);
140 35216 : reduce(v[0],v[1]);
141 35216 : double b11=modulo2(v[0]);
142 35216 : double b22=modulo2(v[1]);
143 35216 : double b12=dotProduct(v[0],v[1]);
144 35216 : double b13=dotProduct(v[0],v[2]);
145 35216 : double b23=dotProduct(v[1],v[2]);
146 35216 : double z=b11*b22-b12*b12;
147 35216 : double y2=-(b11*b23-b12*b13)/z;
148 35216 : double y1=-(b22*b13-b12*b23)/z;
149 35216 : int x1min=floor(y1);
150 35216 : int x1max=x1min+1;
151 35216 : int x2min=floor(y2);
152 35216 : int x2max=x2min+1;
153 : bool first=true;
154 : double mbest,mtrial;
155 35216 : Vector trial,best;
156 105648 : for(int x1=x1min; x1<=x1max; x1++)
157 211296 : for(int x2=x2min; x2<=x2max; x2++) {
158 140864 : trial=v[2]+x2*v[1]+x1*v[0];
159 140864 : mtrial=modulo2(trial);
160 140864 : if(first || mtrial<mbest) {
161 : mbest=mtrial;
162 42378 : best=trial;
163 : first=false;
164 : }
165 : }
166 35216 : if(modulo2(best)*onePlusEpsilon>=modulo2(v[2])) {
167 : break;
168 : }
169 14973 : counter++;
170 14973 : if(counter%10000==0) {
171 0 : std::fprintf(stderr,"WARNING: LatticeReduction::reduceFast stuck after %u iterations\n",counter);
172 : }
173 14973 : v[2]=best;
174 14973 : }
175 20243 : sort(v);
176 20243 : t.setRow(0,v[0]);
177 20243 : t.setRow(1,v[1]);
178 20243 : t.setRow(2,v[2]);
179 20243 : }
180 :
181 :
182 1 : void LatticeReduction::reduceSlow(Tensor&t) {
183 4 : Vector v[3];
184 1 : v[0]=t.getRow(0);
185 1 : v[1]=t.getRow(1);
186 1 : v[2]=t.getRow(2);
187 1 : reduce2(v[0],v[1],v[2]);
188 1 : double e01=dotProduct(v[0],v[1]);
189 1 : double e02=dotProduct(v[0],v[2]);
190 1 : double e12=dotProduct(v[1],v[2]);
191 1 : if(e01*e02*e12<0) {
192 : int eps01=0;
193 0 : if(e01>0.0) {
194 : eps01=1;
195 0 : } else if(e01<0.0) {
196 : eps01=-1;
197 : }
198 : int eps02=0;
199 0 : if(e02>0.0) {
200 : eps02=1;
201 0 : } else if(e02<0.0) {
202 : eps02=-1;
203 : }
204 0 : Vector n=v[0]-eps01*v[1]-eps02*v[2];
205 : int i=0;
206 0 : double mx=modulo2(v[i]);
207 0 : for(int j=1; j<3; j++) {
208 0 : double f=modulo2(v[j]);
209 0 : if(f>mx) {
210 : i=j;
211 : mx=f;
212 : }
213 : }
214 0 : if(modulo2(n)<mx) {
215 0 : v[i]=n;
216 : }
217 : }
218 1 : sort(v);
219 1 : t.setRow(0,v[0]);
220 1 : t.setRow(1,v[1]);
221 1 : t.setRow(2,v[2]);
222 1 : }
223 :
224 0 : bool LatticeReduction::isReduced2(const Vector&a,const Vector&b,const Vector &c) {
225 0 : return isReduced(a,b) && isReduced(a,b) && isReduced(b,c);
226 : }
227 :
228 2 : bool LatticeReduction::isReduced(const Tensor&t) {
229 8 : Vector v[3];
230 : double m[3];
231 2 : v[0]=t.getRow(0);
232 2 : v[1]=t.getRow(1);
233 2 : v[2]=t.getRow(2);
234 8 : for(int i=0; i<3; i++) {
235 6 : m[i]=modulo2(v[i]);
236 : }
237 2 : if(!((m[0]<=m[1]) && m[1]<=m[2])) {
238 : return false;
239 : }
240 : const int cut=5;
241 13 : for(int i=-cut; i<=cut; i++) {
242 12 : double mm=modulo2(v[1]+i*v[0]);
243 12 : if(mm<m[1]) {
244 : return false;
245 : }
246 142 : for(int j=-cut; j<=cut; j++) {
247 131 : double mx=modulo2(v[2]+i*v[1]+j*v[0]);
248 131 : if(mx<m[2]) {
249 : return false;
250 : }
251 : }
252 : }
253 : return true;
254 : }
255 :
256 : }
|
