Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2019 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 :
26 : namespace PLMD {
27 :
28 : using namespace std;
29 :
30 : const double epsilon=1e-14;
31 :
32 24666 : void LatticeReduction::sort(Vector v[3]) {
33 : const double onePlusEpsilon=(1.0+epsilon);
34 172662 : for(int i=0; i<3; i++) for(int j=i+1; j<3; j++) if(modulo2(v[i])>modulo2(v[j])) {
35 9957 : Vector x=v[i]; v[i]=v[j]; v[j]=x;
36 : }
37 73998 : for(int i=0; i<2; i++) plumed_assert(modulo2(v[i])<=modulo2(v[i+1])*onePlusEpsilon);
38 24666 : }
39 :
40 14670 : void LatticeReduction::reduce(Vector&a,Vector&b) {
41 : const double onePlusEpsilon=(1.0+epsilon);
42 14670 : double ma=modulo2(a);
43 14670 : double mb=modulo2(b);
44 : unsigned counter=0;
45 : while(true) {
46 15247 : if(mb>ma) {
47 9029 : Vector t(a); a=b; b=t;
48 : double mt(ma); ma=mb; mb=mt;
49 : }
50 15247 : a-=b*floor(dotProduct(a,b)/modulo2(b)+0.5);
51 15247 : ma=modulo2(a);
52 15247 : if(mb<=ma*onePlusEpsilon) break;
53 577 : counter++;
54 577 : if(counter%10000==0) fprintf(stderr,"WARNING: LatticeReduction::reduce stuck after %u iterations\n",counter);
55 : }
56 :
57 14670 : Vector t(a); a=b; b=t;
58 14670 : }
59 :
60 0 : void LatticeReduction::reduce2(Vector&a,Vector&b,Vector&c) {
61 0 : Vector v[3];
62 0 : v[0]=a; v[1]=b; v[2]=c;
63 : int iter=0;
64 : int ok=0;
65 0 : while(ok<3) {
66 : int i,j;
67 0 : if(iter%3==0) {
68 : i=0; j=1;
69 0 : } else if(iter%3==1) {
70 : i=0; j=2;
71 : } else {
72 : i=1; j=2;
73 : }
74 0 : if(isReduced(v[i],v[j])) ok++;
75 : else {
76 0 : reduce(v[i],v[j]);
77 : ok=1;
78 : }
79 0 : iter++;
80 : }
81 0 : a=v[0]; b=v[1]; c=v[2];
82 0 : }
83 :
84 0 : bool LatticeReduction::isReduced(const Vector&a,const Vector&b) {
85 : const int cut=5;
86 0 : for(int i=-cut; i<=cut; i++) {
87 0 : if(modulo2(b+i*a)<modulo2(b)) return false;
88 : }
89 0 : return modulo2(a)<=modulo2(b) && 2.0*dotProduct(a,b)<=modulo2(a);
90 : }
91 :
92 0 : void LatticeReduction::reduce2(Tensor&t) {
93 0 : Vector a=t.getRow(0);
94 0 : Vector b=t.getRow(1);
95 0 : Vector c=t.getRow(2);
96 0 : reduce2(a,b,c);
97 0 : t.setRow(0,a);
98 0 : t.setRow(1,b);
99 0 : t.setRow(2,c);
100 0 : }
101 :
102 9996 : void LatticeReduction::reduce(Tensor&t) {
103 9996 : reduceFast(t);
104 9996 : }
105 :
106 9996 : void LatticeReduction::reduceFast(Tensor&t) {
107 : const double onePlusEpsilon=(1.0+epsilon);
108 39984 : Vector v[3];
109 9996 : v[0]=t.getRow(0);
110 9996 : v[1]=t.getRow(1);
111 9996 : v[2]=t.getRow(2);
112 : unsigned counter=0;
113 : while(true) {
114 14670 : sort(v);
115 14670 : reduce(v[0],v[1]);
116 14670 : double b11=modulo2(v[0]);
117 14670 : double b22=modulo2(v[1]);
118 14670 : double b12=dotProduct(v[0],v[1]);
119 14670 : double b13=dotProduct(v[0],v[2]);
120 14670 : double b23=dotProduct(v[1],v[2]);
121 14670 : double z=b11*b22-b12*b12;
122 14670 : double y2=-(b11*b23-b12*b13)/z;
123 14670 : double y1=-(b22*b13-b12*b23)/z;
124 14670 : int x1min=floor(y1);
125 14670 : int x1max=x1min+1;
126 14670 : int x2min=floor(y2);
127 14670 : int x2max=x2min+1;
128 : bool first=true;
129 : double mbest,mtrial;
130 14670 : Vector trial,best;
131 73350 : for(int x1=x1min; x1<=x1max; x1++)
132 146700 : for(int x2=x2min; x2<=x2max; x2++) {
133 58680 : trial=v[2]+x2*v[1]+x1*v[0];
134 58680 : mtrial=modulo2(trial);
135 58680 : if(first || mtrial<mbest) {
136 : mbest=mtrial;
137 21588 : best=trial;
138 : first=false;
139 : }
140 : }
141 14670 : if(modulo2(best)*onePlusEpsilon>=modulo2(v[2])) break;
142 4674 : counter++;
143 4674 : if(counter%10000==0) fprintf(stderr,"WARNING: LatticeReduction::reduceFast stuck after %u iterations\n",counter);
144 4674 : v[2]=best;
145 4674 : }
146 9996 : sort(v);
147 9996 : t.setRow(0,v[0]);
148 9996 : t.setRow(1,v[1]);
149 9996 : t.setRow(2,v[2]);
150 9996 : }
151 :
152 :
153 0 : void LatticeReduction::reduceSlow(Tensor&t) {
154 0 : Vector v[3];
155 0 : v[0]=t.getRow(0);
156 0 : v[1]=t.getRow(1);
157 0 : v[2]=t.getRow(2);
158 0 : reduce2(v[0],v[1],v[2]);
159 0 : double e01=dotProduct(v[0],v[1]);
160 0 : double e02=dotProduct(v[0],v[2]);
161 0 : double e12=dotProduct(v[1],v[2]);
162 0 : if(e01*e02*e12<0) {
163 0 : int eps01=0; if(e01>0.0) eps01=1; else if(e01<0.0) eps01=-1;
164 0 : int eps02=0; if(e02>0.0) eps02=1; else if(e02<0.0) eps02=-1;
165 0 : Vector n=v[0]-eps01*v[1]-eps02*v[2];
166 0 : int i=0; double mx=modulo2(v[i]);
167 0 : for(int j=1; j<3; j++) {
168 0 : double f=modulo2(v[j]);
169 0 : if(f>mx) {
170 : i=j;
171 : mx=f;
172 : }
173 : }
174 0 : if(modulo2(n)<mx) v[i]=n;
175 : }
176 0 : sort(v);
177 0 : t.setRow(0,v[0]);
178 0 : t.setRow(1,v[1]);
179 0 : t.setRow(2,v[2]);
180 0 : }
181 :
182 0 : bool LatticeReduction::isReduced2(const Vector&a,const Vector&b,const Vector &c) {
183 0 : return isReduced(a,b) && isReduced(a,b) && isReduced(b,c);
184 : }
185 :
186 0 : bool LatticeReduction::isReduced(const Tensor&t) {
187 0 : Vector v[3];
188 : double m[3];
189 0 : v[0]=t.getRow(0);
190 0 : v[1]=t.getRow(1);
191 0 : v[2]=t.getRow(2);
192 0 : for(int i=0; i<3; i++) m[i]=modulo2(v[i]);
193 0 : if(!((m[0]<=m[1]) && m[1]<=m[2])) return false;
194 : const int cut=5;
195 0 : for(int i=-cut; i<=cut; i++) {
196 0 : double mm=modulo2(v[1]+i*v[0]);
197 0 : if(mm<m[1]) return false;
198 0 : for(int j=-cut; j<=cut; j++) {
199 0 : double mx=modulo2(v[2]+i*v[1]+j*v[0]);
200 0 : if(mx<m[2])return false;
201 : }
202 : }
203 : return true;
204 : }
205 :
206 : }
|
