DGtal  0.9.4.1
testIntegerComputer.cpp
Go to the documentation of this file.
1 
30 #include <cstdlib>
32 #include <iostream>
33 #include "DGtal/base/Common.h"
34 #include "DGtal/arithmetic/IntegerComputer.h"
36 
37 using namespace std;
38 using namespace DGtal;
39 
41 // Functions for testing class IntegerComputer.
43 
44 template <typename Integer>
46 {
47  unsigned int nbok = 0;
48  unsigned int nb = 0;
49  Integer a = rand();
50  Integer b = rand();
51  Integer g = ic.gcd( a, b );
52  trace.info() << "GCD(" << a << "," << b << ")"
53  << " = " << g << std::endl;
54  Integer ra = a % g;
55  Integer rb = b % g;
56  nbok += ic.isZero( ra ) ? 1 : 0;
57  nb++;
58  trace.info() << "(" << nbok << "/" << nb << ") "
59  << a << " % " << g << " == 0" << std::endl;
60  nbok += ic.isZero( rb ) ? 1 : 0;
61  nb++;
62  trace.info() << "(" << nbok << "/" << nb << ") "
63  << b << " % " << g << " == 0" << std::endl;
64  a /= g; b /= g;
65  g = ic.gcd( a, b );
66  nbok += g == Integer( 1 ) ? 1 : 0;
67  nb++;
68  trace.info() << "(" << nbok << "/" << nb << ") "
69  << "GCD(" << a << "," << b << ") == 1" << std::endl;
70  Integer c = rand();
71  ++c; // avoids zero.
72  a *= c; b *= c;
73  ic.getGcd( g, a, b );
74  nbok += g == c ? 1 : 0;
75  nb++;
76  trace.info() << "(" << nbok << "/" << nb << ") "
77  << "GCD(" << a << "," << b << ") == " << c << std::endl;
78  return nbok == nb;
79 }
80 
81 template <typename Integer>
83 {
84  unsigned int nbok = 0;
85  unsigned int nb = 0;
86  Integer a = rand();
87  Integer b = rand();
88  Integer g = ic.gcd( a, b );
89  trace.info() << "a / b = " << a << " / " << b << std::endl;
90  std::vector<Integer> quotients;
91  Integer g2 = ic.getCFrac( quotients, a, b );
92  nbok += g == g2 ? 1 : 0;
93  nb++;
94  trace.info() << "(" << nbok << "/" << nb << ") "
95  << g << " == " << g2 << std::endl;
96  trace.info() << a << " / " << b << " = ";
97  for ( typename std::vector<Integer>::const_iterator it = quotients.begin(),
98  it_end = quotients.end(); it != it_end; ++it )
99  trace.info() << *it;
100  trace.info() << std::endl;
101  double da = NumberTraits<Integer>::castToDouble( a );
102  double db = NumberTraits<Integer>::castToDouble( b );
103  double q = floor( da / db );
104  nbok += Integer( (int) q ) == quotients[ 0 ] ? 1 : 0;
105  nb++;
106  trace.info() << "(" << nbok << "/" << nb << ") "
107  << q << " == " << quotients[ 0 ] << std::endl;
108  typedef typename IntegerComputer<Integer>::Point2I Point2I;
109  Point2I p = ic.convergent( quotients, quotients.size() );
110  nbok += p[ 0 ] == ( a / g ) ? 1 : 0;
111  nb++;
112  trace.info() << "(" << nbok << "/" << nb << ") "
113  << "convergent p[ 0 ] " << p[ 0 ] << " == a / g " << std::endl;
114  nbok += p[ 1 ] == ( b / g ) ? 1 : 0;
115  nb++;
116  trace.info() << "(" << nbok << "/" << nb << ") "
117  << "convergent p[ 1 ] " << p[ 1 ] << " == b / g " << std::endl;
118  return nbok == nb;
119 }
120 
121 template <typename Integer>
123 {
124  unsigned int nbok = 0;
125  unsigned int nb = 0;
126  Integer a = rand();
127  a -= rand();
128  Integer b = rand();
129  b -= rand();
130  if ( ic.isZero( b ) ) ++b;
131  trace.info() << "- a / b = " << a << " / " << b << std::endl;
132  Integer fl = ic.floorDiv( a, b );
133  Integer ce = ic.ceilDiv( a, b );
134  Integer fl2, ce2;
135  ic.getFloorCeilDiv( fl2, ce2, a, b );
136  nbok += fl == fl2 ? 1 : 0;
137  nb++;
138  trace.info() << "(" << nbok << "/" << nb << ") "
139  << "fl == fl2 " << fl2 << std::endl;
140  nbok += ce == ce2 ? 1 : 0;
141  nb++;
142  trace.info() << "(" << nbok << "/" << nb << ") "
143  << "ce == ce2 " << ce2 << std::endl;
144  Integer m = a % b;
145  nbok += ( ( m == 0 ) && ( fl == ce ) ) || ( fl + 1 == ce ) ? 1 : 0;
146  nb++;
147  trace.info() << "(" << nbok << "/" << nb << ") "
148  << "( ( m == 0 ) && ( fl == ce ) ) || ( fl+1 == ce )"
149  << std::endl;
150  return nbok == nb;
151 }
152 
153 template <typename Integer>
155 {
156  typedef typename IntegerComputer<Integer>::Point2I Point2I;
157  unsigned int nbok = 0;
158  unsigned int nb = 0;
159  Integer a = rand();
160  Integer b = rand();
161  Integer g = ic.gcd( a, b );
162  trace.info() << "a / b = " << a << " / " << b
163  << " gcd=" << g << std::endl;
164  Point2I v = ic.extendedEuclid( a, b, g );
165  trace.info() << "Bezout = " << v[ 0 ] << "," << v[ 1 ] << std::endl;
166  Integer rem = a * v[ 0 ] + b * v[ 1 ];
167  nbok += rem == g ? 1 : 0;
168  nb++;
169  trace.info() << "(" << nbok << "/" << nb << ") "
170  << "a * v[ 0 ] + b * v[ 1 ] == g "
171  << "(" << rem << " == " << g << ")" << std::endl;
172  return nbok == nb;
173 }
174 
175 template <typename Integer>
177 {
178  typedef typename IntegerComputer<Integer>::Point2I Point2I;
179  unsigned int nbok = 0;
180  unsigned int nb = 0;
181  Point2I p, N;
182  Integer c;
183  p = Point2I( rand(), rand() );
184  N = Point2I( rand() , rand() );
185  c = rand() * rand();
186  Point2I u( rand() / 100, rand() / 100);
187  trace.info() << "p = " << p << std::endl;
188  trace.info() << "u = " << u << std::endl;
189  trace.info() << "N = " << N << std::endl;
190  trace.info() << "c = " << c << std::endl;
191  Integer fl, ce;
192  ic.getCoefficientIntersection( fl, ce, p, u, N, c );
193  trace.info() << "fl = " << fl << ", ce = " << ce << std::endl;
194  Point2I p1 = p + (u * fl);
195  Point2I p2 = p + (u * ce);
196  Integer c1 = ic.dotProduct( p1, N );
197  Integer c2 = ic.dotProduct( p2, N );
198  trace.info() << "c1 = " << c1
199  << " <= c = " << c
200  << " <= c2 = " << c2 << std::endl;
201  nbok += ( ( c1 == c2 ) && ( c == c1 ) )
202  || ( ( c1 <= c ) && ( c < c2 ) ) ? 1 : 0;
203  nb++;
204  trace.info() << "(" << nbok << "/" << nb << ") "
205  << "( ( c1 == c2 ) && ( c == c1 ) ) || ( ( c1 <= c ) && ( c < c2 ) )" << std::endl;
206  return nbok == nb;
207 }
208 
209 template <typename Integer>
211 {
212  typedef typename IntegerComputer<Integer>::Point2I Point2I;
213  typedef typename IntegerComputer<Integer>::Vector2I Vector2I;
214  unsigned int nbok = 0;
215  unsigned int nb = 0;
216  Vector2I v;
217  Point2I A( rand(), rand() );
218  Vector2I u( rand() / 100, rand() / 100 );
219  ic.reduce( u );
220  Vector2I N( rand(), rand() );
221  Vector2I N2( rand(), rand() );
222  Integer c = rand() * rand();
223  Integer c2 = rand() * rand();
224  trace.info() << "A = " << A << std::endl;
225  trace.info() << "u = " << u << std::endl;
226  trace.info() << "N = " << N << std::endl;
227  trace.info() << "c = " << c << std::endl;
228  trace.info() << "N2 = " << N2 << std::endl;
229  trace.info() << "c2 = " << c2 << std::endl;
230  ic.getValidBezout ( v,
231  A, u, N, c, N2, c2, true );
232  trace.info() << "-> v = " << v << std::endl;
233  Integer a0 = ic.crossProduct(u,v);
234  nbok += ( ic.abs( a0 ) == 1 ) ? 1 : 0;
235  ++nb;
236  trace.info() << "(" << nbok << "/" << nb << ") "
237  << " v^u = " << a0 << " == +/- 1 " << std::endl;
238  // (A+v).N <= c , (A+v).N2 <= c2, (A+v+u).N2 > c2.
239  // Integer a1 = ic.dotProduct( A+v, N );
240  Integer a2 = ic.dotProduct( A+v, N2 );
241  Integer a3 = ic.dotProduct( A+v+u, N2 );
242  // nbok += a1 <= c ? 1 : 0;
243  // ++nb;
244  // trace.info() << "(" << nbok << "/" << nb << ") "
245  // << "(A+v).N = " << a1
246  // << " <= c = " << c << std::endl;
247  nbok += a2 <= c2 ? 1 : 0;
248  ++nb;
249  trace.info() << "(" << nbok << "/" << nb << ") "
250  << "(A+v).N2 = " << a2
251  << " <= c = " << c2 << std::endl;
252  nbok += a3 > c2 ? 1 : 0;
253  ++nb;
254  trace.info() << "(" << nbok << "/" << nb << ") "
255  << "(A+v+u).N2 = " << a2
256  << " > c2 = " << c2 << std::endl;
257  return nbok == nb;
258 }
259 
265 {
266  unsigned int nbtests = 50;
267  unsigned int nbok = 0;
268  unsigned int nb = 0;
269  typedef BigInteger Integer;
271  trace.beginBlock ( "Testing block: multiple random gcd." );
272  for ( unsigned int i = 0; i < nbtests; ++i )
273  {
274  nbok += testGCD<Integer>( ic ) ? 1 : 0;
275  nb++;
276  }
277  trace.info() << "(" << nbok << "/" << nb << ") gcd tests." << std::endl;
278  trace.endBlock();
279 
280  trace.beginBlock ( "Testing block: multiple random cfrac." );
281  for ( unsigned int i = 0; i < nbtests; ++i )
282  {
283  nbok += testCFrac<Integer>( ic ) ? 1 : 0;
284  nb++;
285  }
286  trace.info() << "(" << nbok << "/" << nb << ") cfrac tests." << std::endl;
287  trace.endBlock();
288 
289  trace.beginBlock ( "Testing block: multiple ceil floor division." );
290  for ( unsigned int i = 0; i < nbtests; ++i )
291  {
292  nbok += testCeilFloorDiv<Integer>( ic ) ? 1 : 0;
293  nb++;
294  }
295  trace.info() << "(" << nbok << "/" << nb << ") ceil floor division." << std::endl;
296  trace.endBlock();
297 
298  trace.beginBlock ( "Testing block: multiple Bezout / extended euclid: ax+by= gcd(a,b)." );
299  for ( unsigned int i = 0; i < nbtests; ++i )
300  {
301  nbok += testExtendedEuclid<Integer>( ic ) ? 1 : 0;
302  nb++;
303  }
304  trace.info() << "(" << nbok << "/" << nb << ") Bezout / extended euclid." << std::endl;
305  trace.endBlock();
306 
307  trace.beginBlock ( "Testing block: multiple coefficient intersection." );
308  for ( unsigned int i = 0; i < nbtests; ++i )
309  {
310  nbok += testCoefficientIntersection( ic ) ? 1 : 0;
311  nb++;
312  }
313  trace.info() << "(" << nbok << "/" << nb << ") coefficient intersection." << std::endl;
314  trace.endBlock();
315 
316  trace.beginBlock ( "Testing block: multiple valid bezout." );
317  for ( unsigned int i = 0; i < nbtests; ++i )
318  {
319  nbok += testValidBezout( ic ) ? 1 : 0;
320  nb++;
321  }
322  trace.info() << "(" << nbok << "/" << nb << ") valid bezout." << std::endl;
323  trace.endBlock();
324 
325  return nbok == nb;
326 }
327 
329 // Standard services - public :
330 
331 int main( int /*argc*/, char** /*argv*/ )
332 {
333  trace.beginBlock ( "Testing class IntegerComputer" );
334  bool res = testIntegerComputer(); // && ... other tests
335  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
336  trace.endBlock();
337  return res ? 0 : 1;
338 }
339 // //
void beginBlock(const std::string &keyword="")
Integer dotProduct(const Vector2I &u, const Vector2I &v) const
void getGcd(Integer &g, IntegerParamType a, IntegerParamType b) const
Trace trace
Definition: Common.h:137
void getValidBezout(Vector2I &v, const Point2I &A, const Vector2I &u, const Vector2I &N, IntegerParamType c, const Vector2I &N2, IntegerParamType c2, bool compute_v=true) const
bool testIntegerComputer()
STL namespace.
double endBlock()
Vector2I extendedEuclid(IntegerParamType a, IntegerParamType b, IntegerParamType c) const
bool testExtendedEuclid(const IntegerComputer< Integer > &ic)
Integer floorDiv(IntegerParamType na, IntegerParamType nb) const
mpz_class BigInteger
Multi-precision integer with GMP implementation.
Definition: BasicTypes.h:79
Integer crossProduct(const Vector2I &u, const Vector2I &v) const
static bool isZero(IntegerParamType a)
static Integer abs(IntegerParamType a)
int main(int, char **)
bool testCoefficientIntersection(const IntegerComputer< Integer > &ic)
bool testValidBezout(const IntegerComputer< Integer > &ic)
std::ostream & emphase()
Point2I convergent(const std::vector< Integer > &quotients, unsigned int k) const
Aim: The traits class for all models of Cinteger.
Definition: NumberTraits.h:69
DGtal is the top-level namespace which contains all DGtal functions and types.
bool testGCD(const IntegerComputer< Integer > &ic)
std::ostream & info()
void getFloorCeilDiv(Integer &fl, Integer &ce, IntegerParamType na, IntegerParamType nb) const
Integer getCFrac(std::vector< Integer > &quotients, IntegerParamType a, IntegerParamType b) const
void reduce(Vector2I &p) const
Integer ceilDiv(IntegerParamType na, IntegerParamType nb) const
bool testCFrac(const IntegerComputer< Integer > &ic)
void getCoefficientIntersection(Integer &fl, Integer &ce, const Vector2I &p, const Vector2I &u, const Vector2I &N, IntegerParamType c) const
Integer gcd(IntegerParamType a, IntegerParamType b) const
bool testCeilFloorDiv(const IntegerComputer< Integer > &ic)