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 31782 : void LatticeReduction::sort(Vector v[3]) {
32 : const double onePlusEpsilon=(1.0+epsilon);
33 : double m[3];
34 31782 : m[0]=modulo2(v[0]);
35 31782 : m[1]=modulo2(v[1]);
36 31782 : m[2]=modulo2(v[2]);
37 222474 : for(int i=0; i<3; i++) for(int j=i+1; j<3; j++) if(m[i]>m[j]) {
38 11780 : std::swap(v[i],v[j]);
39 : std::swap(m[i],m[j]);
40 : }
41 31782 : plumed_assert(m[0]<=m[1]*onePlusEpsilon);
42 31782 : plumed_assert(m[1]<=m[2]*onePlusEpsilon);
43 31782 : }
44 :
45 19233 : void LatticeReduction::reduce(Vector&a,Vector&b) {
46 : const double onePlusEpsilon=(1.0+epsilon);
47 19233 : double ma=modulo2(a);
48 19233 : double mb=modulo2(b);
49 : unsigned counter=0;
50 : while(true) {
51 19799 : if(mb>ma) {
52 : std::swap(a,b);
53 : std::swap(ma,mb);
54 : }
55 19799 : a-=b*floor(dotProduct(a,b)/mb+0.5);
56 19799 : ma=modulo2(a);
57 19799 : if(mb<=ma*onePlusEpsilon) break;
58 566 : counter++;
59 566 : if(counter%100==0) { // only test rarely since this might be expensive
60 0 : plumed_assert(!std::isnan(ma));
61 0 : plumed_assert(!std::isnan(mb));
62 : }
63 566 : if(counter%10000==0) std::fprintf(stderr,"WARNING: LatticeReduction::reduce stuck after %u iterations\n",counter);
64 : }
65 :
66 : std::swap(a,b);
67 19233 : }
68 :
69 1 : void LatticeReduction::reduce2(Vector&a,Vector&b,Vector&c) {
70 4 : Vector v[3];
71 1 : v[0]=a; v[1]=b; v[2]=c;
72 : int iter=0;
73 : int ok=0;
74 12 : while(ok<3) {
75 : int i,j;
76 11 : if(iter%3==0) {
77 : i=0; j=1;
78 7 : } else if(iter%3==1) {
79 : i=0; j=2;
80 : } else {
81 : i=1; j=2;
82 : }
83 11 : if(isReduced(v[i],v[j])) ok++;
84 : else {
85 7 : reduce(v[i],v[j]);
86 : ok=1;
87 : }
88 11 : iter++;
89 : }
90 1 : a=v[0]; b=v[1]; c=v[2];
91 1 : }
92 :
93 11 : bool LatticeReduction::isReduced(const Vector&a,const Vector&b) {
94 : const int cut=5;
95 84 : for(int i=-cut; i<=cut; i++) {
96 79 : if(modulo2(b+i*a)<modulo2(b)) return false;
97 : }
98 5 : return modulo2(a)<=modulo2(b) && 2.0*dotProduct(a,b)<=modulo2(a);
99 : }
100 :
101 0 : void LatticeReduction::reduce2(Tensor&t) {
102 0 : Vector a=t.getRow(0);
103 0 : Vector b=t.getRow(1);
104 0 : Vector c=t.getRow(2);
105 0 : reduce2(a,b,c);
106 0 : t.setRow(0,a);
107 0 : t.setRow(1,b);
108 0 : t.setRow(2,c);
109 0 : }
110 :
111 12554 : void LatticeReduction::reduce(Tensor&t) {
112 12554 : reduceFast(t);
113 12554 : }
114 :
115 12555 : void LatticeReduction::reduceFast(Tensor&t) {
116 : const double onePlusEpsilon=(1.0+epsilon);
117 50220 : Vector v[3];
118 12555 : v[0]=t.getRow(0);
119 12555 : v[1]=t.getRow(1);
120 12555 : v[2]=t.getRow(2);
121 : unsigned counter=0;
122 : while(true) {
123 19226 : sort(v);
124 19226 : reduce(v[0],v[1]);
125 19226 : double b11=modulo2(v[0]);
126 19226 : double b22=modulo2(v[1]);
127 19226 : double b12=dotProduct(v[0],v[1]);
128 19226 : double b13=dotProduct(v[0],v[2]);
129 19226 : double b23=dotProduct(v[1],v[2]);
130 19226 : double z=b11*b22-b12*b12;
131 19226 : double y2=-(b11*b23-b12*b13)/z;
132 19226 : double y1=-(b22*b13-b12*b23)/z;
133 19226 : int x1min=floor(y1);
134 19226 : int x1max=x1min+1;
135 19226 : int x2min=floor(y2);
136 19226 : int x2max=x2min+1;
137 : bool first=true;
138 : double mbest,mtrial;
139 19226 : Vector trial,best;
140 57678 : for(int x1=x1min; x1<=x1max; x1++)
141 115356 : for(int x2=x2min; x2<=x2max; x2++) {
142 76904 : trial=v[2]+x2*v[1]+x1*v[0];
143 76904 : mtrial=modulo2(trial);
144 76904 : if(first || mtrial<mbest) {
145 : mbest=mtrial;
146 26770 : best=trial;
147 : first=false;
148 : }
149 : }
150 19226 : if(modulo2(best)*onePlusEpsilon>=modulo2(v[2])) break;
151 6671 : counter++;
152 6671 : if(counter%10000==0) std::fprintf(stderr,"WARNING: LatticeReduction::reduceFast stuck after %u iterations\n",counter);
153 6671 : v[2]=best;
154 6671 : }
155 12555 : sort(v);
156 12555 : t.setRow(0,v[0]);
157 12555 : t.setRow(1,v[1]);
158 12555 : t.setRow(2,v[2]);
159 12555 : }
160 :
161 :
162 1 : void LatticeReduction::reduceSlow(Tensor&t) {
163 4 : Vector v[3];
164 1 : v[0]=t.getRow(0);
165 1 : v[1]=t.getRow(1);
166 1 : v[2]=t.getRow(2);
167 1 : reduce2(v[0],v[1],v[2]);
168 1 : double e01=dotProduct(v[0],v[1]);
169 1 : double e02=dotProduct(v[0],v[2]);
170 1 : double e12=dotProduct(v[1],v[2]);
171 1 : if(e01*e02*e12<0) {
172 0 : int eps01=0; if(e01>0.0) eps01=1; else if(e01<0.0) eps01=-1;
173 0 : int eps02=0; if(e02>0.0) eps02=1; else if(e02<0.0) eps02=-1;
174 0 : Vector n=v[0]-eps01*v[1]-eps02*v[2];
175 0 : int i=0; double mx=modulo2(v[i]);
176 0 : for(int j=1; j<3; j++) {
177 0 : double f=modulo2(v[j]);
178 0 : if(f>mx) {
179 : i=j;
180 : mx=f;
181 : }
182 : }
183 0 : if(modulo2(n)<mx) v[i]=n;
184 : }
185 1 : sort(v);
186 1 : t.setRow(0,v[0]);
187 1 : t.setRow(1,v[1]);
188 1 : t.setRow(2,v[2]);
189 1 : }
190 :
191 0 : bool LatticeReduction::isReduced2(const Vector&a,const Vector&b,const Vector &c) {
192 0 : return isReduced(a,b) && isReduced(a,b) && isReduced(b,c);
193 : }
194 :
195 2 : bool LatticeReduction::isReduced(const Tensor&t) {
196 8 : Vector v[3];
197 : double m[3];
198 2 : v[0]=t.getRow(0);
199 2 : v[1]=t.getRow(1);
200 2 : v[2]=t.getRow(2);
201 8 : for(int i=0; i<3; i++) m[i]=modulo2(v[i]);
202 2 : if(!((m[0]<=m[1]) && m[1]<=m[2])) return false;
203 : const int cut=5;
204 13 : for(int i=-cut; i<=cut; i++) {
205 12 : double mm=modulo2(v[1]+i*v[0]);
206 12 : if(mm<m[1]) return false;
207 142 : for(int j=-cut; j<=cut; j++) {
208 131 : double mx=modulo2(v[2]+i*v[1]+j*v[0]);
209 131 : if(mx<m[2])return false;
210 : }
211 : }
212 : return true;
213 : }
214 :
215 : }
|
