DGtal 1.3.0
Loading...
Searching...
No Matches
testIntegerComputer.cpp
Go to the documentation of this file.
1
31#include <cstdlib>
32#include <iostream>
33#include "DGtal/base/Common.h"
34#include "DGtal/arithmetic/IntegerComputer.h"
36
37using namespace std;
38using namespace DGtal;
39
41// Functions for testing class IntegerComputer.
43
44template <typename Integer>
45bool testGCD( const IntegerComputer<Integer> & ic )
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
81template <typename Integer>
82bool testCFrac( const IntegerComputer<Integer> & ic )
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
121template <typename Integer>
122bool testCeilFloorDiv( const IntegerComputer<Integer> & ic )
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
153template <typename Integer>
154bool testExtendedEuclid( const IntegerComputer<Integer> & ic )
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
175template <typename Integer>
176bool testCoefficientIntersection( const IntegerComputer<Integer> & ic )
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
209template <typename Integer>
210bool testValidBezout( const IntegerComputer<Integer> & ic )
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
264bool testIntegerComputer()
265{
266 unsigned int nbtests = 50;
267 unsigned int nbok = 0;
268 unsigned int nb = 0;
269 typedef BigInteger Integer;
270 IntegerComputer<Integer> ic;
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
331int 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// //
DGtal::IntegerComputer
Aim: This class gathers several types and methods to make computation with integers.
Definition: IntegerComputer.h:83
DGtal::IntegerComputer::getCFrac
Integer getCFrac(std::vector< Integer > &quotients, IntegerParamType a, IntegerParamType b) const
DGtal::IntegerComputer::getGcd
void getGcd(Integer &g, IntegerParamType a, IntegerParamType b) const
DGtal::IntegerComputer::crossProduct
Integer crossProduct(const Vector2I &u, const Vector2I &v) const
DGtal::IntegerComputer::isZero
static bool isZero(IntegerParamType a)
DGtal::IntegerComputer::abs
static Integer abs(IntegerParamType a)
DGtal::IntegerComputer::ceilDiv
Integer ceilDiv(IntegerParamType na, IntegerParamType nb) const
DGtal::IntegerComputer::reduce
void reduce(Vector2I &p) const
DGtal::IntegerComputer::Vector2I
SpaceND< 2, Integer >::Vector Vector2I
Definition: IntegerComputer.h:94
DGtal::IntegerComputer::dotProduct
Integer dotProduct(const Vector2I &u, const Vector2I &v) const
DGtal::IntegerComputer::Point2I
SpaceND< 2, Integer >::Point Point2I
Definition: IntegerComputer.h:93
DGtal::IntegerComputer::convergent
Point2I convergent(const std::vector< Integer > &quotients, unsigned int k) const
DGtal::IntegerComputer::extendedEuclid
Vector2I extendedEuclid(IntegerParamType a, IntegerParamType b, IntegerParamType c) const
DGtal::IntegerComputer::getCoefficientIntersection
void getCoefficientIntersection(Integer &fl, Integer &ce, const Vector2I &p, const Vector2I &u, const Vector2I &N, IntegerParamType c) const
DGtal::IntegerComputer::getValidBezout
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
DGtal::IntegerComputer::getFloorCeilDiv
void getFloorCeilDiv(Integer &fl, Integer &ce, IntegerParamType na, IntegerParamType nb) const
DGtal::IntegerComputer::gcd
Integer gcd(IntegerParamType a, IntegerParamType b) const
DGtal::IntegerComputer::floorDiv
Integer floorDiv(IntegerParamType na, IntegerParamType nb) const
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Trace::emphase
std::ostream & emphase()
DGtal::Trace::info
std::ostream & info()
DGtal::Trace::endBlock
double endBlock()
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
Definition: ClosedIntegerHalfPlane.h:49
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::BigInteger
mpz_class BigInteger
Multi-precision integer with GMP implementation.
Definition: BasicTypes.h:79
std
STL namespace.
DGtal::NumberTraits
Aim: The traits class for all models of Cinteger.
Definition: NumberTraits.h:564
main
int main()
Definition: testBits.cpp:56
testGCD
bool testGCD(const IntegerComputer< Integer > &ic)
Definition: testIntegerComputer.cpp:45
testCeilFloorDiv
bool testCeilFloorDiv(const IntegerComputer< Integer > &ic)
Definition: testIntegerComputer.cpp:122
testValidBezout
bool testValidBezout(const IntegerComputer< Integer > &ic)
Definition: testIntegerComputer.cpp:210
testCoefficientIntersection
bool testCoefficientIntersection(const IntegerComputer< Integer > &ic)
Definition: testIntegerComputer.cpp:176
testCFrac
bool testCFrac(const IntegerComputer< Integer > &ic)
Definition: testIntegerComputer.cpp:82
testIntegerComputer
bool testIntegerComputer()
Definition: testIntegerComputer.cpp:264
testExtendedEuclid
bool testExtendedEuclid(const IntegerComputer< Integer > &ic)
Definition: testIntegerComputer.cpp:154
Generated on Sun Nov 27 2022 15:15:53 for DGtal by doxygen 1.9.5