5
#include "rheolef/config.h"
6
#include "predicates_tst.h"
7
using namespace rheolef;
12
#ifdef _RHEOLEF_HAVE_GMPXX_H
14
typedef mpf_class real;
15
typedef mpz_class integer;
16
int sign(integer x) { return (x > 0) ? 1 : ((x < 0) ? -1 : 0); }
17
int sign(real x) { return (x > 0) ? 1 : ((x < 0) ? -1 : 0); }
19
#endif // _RHEOLEF_HAVE_GMPXX_H
22
#include <LEDA/real.h>
23
#include <LEDA/integer.h>
24
int sign(real x) { x.sign(1); }
29
extern "C" int getpid();
30
extern "C" long clock();
32
void exact_out(ostream& s, double a)
40
// random value generation
41
// return in a double a signed integer with i bits
44
return (((int) random()%2)?1:-1) *
45
((double) ((int) random()% ((int)exp2(i)) ));
47
return (((int) random()%2)?1:-1) *
48
(((double) ((int) random()% ((int)exp2(i-25)) )) * exp2(25) +
49
((double) ((int) random()% ((int)exp2(25)) )));
53
// Few 2x2 determinant sign computation
54
int quadruple_det2x2(double a, double b, double c, double d)
56
long double DET=(long double)a*(long double)d-(long double)b*(long double)c;
57
if (DET>0.0) return 1;
58
if (DET<0.0) return -1;
62
int double_det2x2(double a, double b, double c, double d)
63
{ return ( ((a*d)<(b*c)) ? -1 : 1) ; }
66
int leda_integer_det2x2(double a, double b, double c, double d)
68
integer D = integer(a)*integer(d)-integer(b)*integer(c);
71
int leda_real_det2x2(double a, double b, double c, double d)
73
real D = real(a)*real(d)-real(b)*real(c);
79
int double_det3x3(double a, double b, double c,
80
double d, double e, double f, double g, double h, double i)
81
{ if ( a*(e*i-f*h) + b*(f*g-d*i) + c*(d*h-e*g) < 0 ) return -1 ;
85
// Few 3x3 determinant sign computation
86
int leda_integer_det3x3(double a, double b, double c,
87
double d, double e, double f,
88
double g, double h, double i)
90
integer D = integer(a)* (integer(e)* integer(i) -integer(f)* integer(h) )
91
+ integer(b)* (integer(f)* integer(g) -integer(d)* integer(i) )
92
+ integer(c)* (integer(d)* integer(h) -integer(e)* integer(g) );
95
int leda_real_det3x3(double a, double b, double c,
96
double d, double e, double f,
97
double g, double h, double i)
99
real D = real(a)* (real(e)* real(i) - real(f)* real(h) )
100
+ real(b)* (real(f)* real(g) - real(d)* real(i) )
101
+ real(c)* (real(d)* real(h) - real(e)* real(g) );
111
void generate2(double& a,double& b,double& c,double& d, char type)
140
case 'f' : // y = floor(u*x) u in [0,1]
143
b = floor(a*cxx/exp2(52));
145
d = floor(c*cxx/exp2(52));
148
case 'c' : // all entries equal
149
a = b = c = d = rnd(53);
153
void generate3(double& a,double& b,double& c,double& d,double& e,
154
double& f,double& g,double& h,double& i, char type)
170
case 'q' : // aU bV cU+dV
175
c = cxx * a + cyy * b;
178
f = cxx * d + cyy * e;
181
i = cxx * g + cyy * h;
192
case 'p' : // x+y+z=0
204
case 'f' : // z = floor( t*x + u*y ) t,u in [0,1]
209
c = floor((-a*cxx-b*cyy)/exp2(52));
212
f = floor((-d*cxx-e*cyy)/exp2(52));
215
i = floor((-g*cxx-h*cyy)/exp2(52));
218
case 'c' : // all entries equal
219
a = b = c = d = e = f = g = h = i = rnd(51);
223
void perturbation2(double& a,double& b,double& c,double& d)
224
{ a += rnd(2); b += rnd(2); c += rnd(2); d += rnd(2);}
226
void perturbation3(double& a,double& b,double& c,double& d,double& e,
227
double& f,double& g,double& h,double& i)
228
{ a += rnd(2); b += rnd(2); c += rnd(2); d += rnd(2); e += rnd(2); f += rnd(2);
229
g += rnd(2); h += rnd(2); i += rnd(2); }
231
void transpose2(double& a,double& b,double& c,double& d)
232
{ double swap; (void) a; (void) d; swap = b; b = c; c = swap; }
234
void transpose3(double& a,double& b,double& c,double& d,double& e,
235
double& f,double& g,double& h,double& i)
236
{ double swap; swap = b; b = d; d = swap; swap = c; c = g; g = swap;
237
swap = f; f = h; h = swap; (void) a; (void) e; (void) i; }
246
double a,b,c,d,e,f,g,h,i;
250
while (1) if (dim == 2) {
251
cout << "a = "; cin >> a;
252
cout << "b = "; cin >> b;
253
cout << "c = "; cin >> c;
254
cout << "d = "; cin >> d;
256
cout << " | " ; exact_out(cout,a);
257
cout << " " ; exact_out(cout,c);
258
cout << " |" << endl;
259
cout << " | " ; exact_out(cout,b);
260
cout << " " ; exact_out(cout,d);
261
cout << " |" << endl;
262
cout << "double :" << double_det2x2(a,b,c,d) << endl;
263
cout << "quadruple :" << quadruple_det2x2(a,b,c,d) << endl;
265
cout << "integer :" << leda_integer_det2x2(a,b,c,d) << endl;
266
cout << "real:" << leda_real_det2x2(a,b,c,d) << endl;
268
cout << "ABDPY not lazy:" << ABDPY__not_lazy_det2x2(a,b,c,d) << endl;
269
cout << "ABDPY lazy:" << ABDPY_det2x2(a,b,c,d) << endl;
270
cout << "ABDPY secure:" << ABDPY_secure_det2x2(a,b,c,d) << endl;
271
} else if (dim == 3) {
272
cout << "a = "; cin >> a;
273
cout << "b = "; cin >> b;
274
cout << "c = "; cin >> c;
275
cout << "d = "; cin >> d;
276
cout << "e = "; cin >> e;
277
cout << "f = "; cin >> f;
278
cout << "g = "; cin >> g;
279
cout << "h = "; cin >> h;
280
cout << "i = "; cin >> i;
281
cout << " | " ; exact_out(cout,a);
282
cout << " " ; exact_out(cout,d);
283
cout << " " ; exact_out(cout,g);
284
cout << " |" << endl;
285
cout << " | " ; exact_out(cout,b);
286
cout << " " ; exact_out(cout,e);
287
cout << " " ; exact_out(cout,h);
288
cout << " |" << endl;
289
cout << " | " ; exact_out(cout,c);
290
cout << " " ; exact_out(cout,f);
291
cout << " " ; exact_out(cout,i);
292
cout << " |" << endl;
293
cout << "double :" << double_det3x3(a,b,c,d,e,f,g,h,i) << endl;
295
cout << "integer :" << leda_integer_det3x3(a,b,c,d,e,f,g,h,i) << endl;
296
cout << "real :" << leda_real_det3x3(a,b,c,d,e,f,g,h,i) << endl;
298
cout << "ABDPY not lazy:"<<ABDPY__not_lazy_det3x3(a,b,c,d,e,f,g,h,i)<<endl;
299
cout << "ABDPY lazy:" << ABDPY_det3x3(a,b,c,d,e,f,g,h,i) << endl;
300
cout << "ABDPY secure:" << ABDPY_secure_det3x3(a,b,c,d,e,f,g,h,i) << endl;
307
int main(int argc, char **argv)
310
double epsilon_mach = exactinit();
311
#endif // USE_SHEWCHUK
314
int error_2d [7] = {0, 0, 0, 0, 0, 0};
315
int error_3d [7] = {0, 0, 0, 0, 0, 0};
316
double a,b,c,d,e,f,g,h,i;
319
int input, method, trial, j;
322
int (*FF2)(double, double, double, double);
323
int (*FF3)(double, double, double, double, double, double, double, double, double);
328
cerr <<"init random with pid ="<<n<<endl;
333
if ( (**argv>='1') && (**argv<='9') ) n = atoi(*argv);
339
cout <<"\\documentstyle{article}" << endl;
340
cout <<"\\voffset=-2.5cm" << endl;
341
cout <<"\\hoffset=-2.5cm" << endl;
342
cout <<"\\textwidth 17.5cm" << endl;
343
cout <<"\\textheight 24.2cm" << endl;
344
cout <<"\\begin{document}" << endl;
345
cout <<"\\title{Test of determinant sign algorithms}" << endl;
346
cout <<"\\author{O. Devillers}" << endl;
347
cout <<"\\maketitle" << endl;
349
cout <<"Few determinant sign algorithms" << endl;
350
cout <<"are tested with various kinds of input." << endl;
351
cout <<"For each kind of input, "<<n<<" differents trials are done." << endl;
352
cout <<"For each, the determinant is computed "<<N<<" times," << endl;
353
cout <<"and the time are obtained using the {\\tt clock} function." << endl;
355
cout <<"Thus the precision on the given running time is about " << endl;
356
cout <<16000/(double)N << "$\\mu$s." << endl;
358
cout <<"The program verify that all methods" << endl;
359
cout <<"produces the same result, except" << endl;
360
cout <<"the direct floating point calculation" << endl;
361
cout <<"which can be wrong due to" << endl;
362
cout <<"rounding error." << endl;
364
cout <<"\\section*{Platform}" << endl;
365
cout <<"% ADD HERE YOUR PLATFORM" << endl;
366
cout <<"Computer :\\\\" << endl;
367
cout <<"Compiler :\\\\" << endl;
368
cout <<"Compiler options : " << endl;
370
cout <<"\\section*{Algorithms}" << endl;
371
cout <<"Here are the different algorithms used," << endl;
372
cout <<"the original data are of {\\tt double} type." << endl;
374
cout <<"\\subsection*{Direct computation}" << endl;
375
cout <<"Direct computation of the determinant" << endl;
376
cout <<"using $ad-bc$ formula in 2D and" << endl;
377
cout <<"$a(ei-fh) + b(fg-di) + c(dh-eg)$ in 3D." << endl;
378
cout <<"Plus an evaluation of the sign." << endl;
379
cout <<"\\subsubsection*{double}" << endl;
380
cout <<"If the computation is done using {\\tt double} arithmetic," << endl;
381
cout <<"it is fast, but there is rounding errors." << endl;
382
cout <<"\\subsubsection*{quadruple}" << endl;
383
cout <<"If the computation is done using quadruple precision" << endl;
384
cout <<"({\\tt long double} type) the result is exact in 2D." << endl;
385
cout <<"\\subsubsection*{leda-integer}" << endl;
386
cout <<"LEDA provides exact computation on integer of arbitrary" << endl;
387
cout <<"length. The result is exact \\cite{leda}." << endl;
388
cout <<"\\subsubsection*{leda-real}" << endl;
389
cout <<"LEDA provides a floating point filter which invokes" << endl;
390
cout <<"an exact computation only when the rounded computation" << endl;
391
cout <<"does not allow a conclusion \\cite{leda}." << endl;
392
cout <<"The result is exact." << endl;
393
cout <<"\\subsection*{Iterative method}" << endl;
394
cout <<"Avnaim, Boissonnat, Devillers, Preparata and Yvinec" << endl;
395
cout <<"provides an iterative method which transform the " << endl;
396
cout <<"determinant in another one with smaller entries" << endl;
397
cout <<"up the entries have sign allowing a conclusion" << endl;
398
cout <<"\\cite{iter1,iter2}." << endl;
399
cout <<"The algorithm uses double arithmetic to perform" << endl;
400
cout <<"fast operations on 53 bits integers, thus there is precondition" << endl;
401
cout <<"that the entries must represent an integer and their" << endl;
402
cout <<"absolute value must be smaller than $2^{53}$ in 2D" << endl;
403
cout <<"and $2^{51}$ in 3D." << endl;
404
cout <<"\\subsubsection*{not-lazy}" << endl;
405
cout <<"The algorithm is iterated, up to the comparison" << endl;
406
cout <<"of the entries and computation of sign of minors" << endl;
407
cout <<"in 3D allows a conclusion." << endl;
408
cout <<"\\subsubsection*{lazy}" << endl;
409
cout <<"The algorithm is combinated with a floating point filter." << endl;
410
cout <<"\\subsubsection*{secure}" << endl;
411
cout <<"The algorithm is combinated with a floating point filter," << endl;
412
cout <<"and the entries are verified to match the preconditions." << endl;
414
cout <<"\\small" << endl;
415
cout <<"\\section*{Input 2D}" << endl;
416
cout <<"$a$ is random on $b$ bits means a is evenly distributed" << endl;
417
cout <<"among the integers between $-2^b+1$ and $2^b-1$." << endl;
418
cout <<"\\subsection*{random}" << endl;
419
cout <<"$\\left| \\begin{array}{cc} a& c\\\\ b & d \\end{array} \\right|" << endl;
420
cout <<"\\hspace{1cm} a,b,c,d \\mbox{ random on 53 bits} $" << endl;
421
cout <<"\\subsection*{$x=-y$}" << endl;
422
cout <<"$\\left| \\begin{array}{cc} a& c\\\\ -a & -c \\end{array} \\right|" << endl;
423
cout <<"\\hspace{1cm} a,c \\mbox{ random on 53 bits} $" << endl;
424
cout <<"\\subsection*{$x=-y+\\varepsilon$}" << endl;
425
cout <<"$\\left| \\begin{array}{cc} a+e_a& c+e_c\\\\ -a+e_b & -c+e_d \\end{array} \\right|" << endl;
426
cout <<"\\hspace{1cm} a,c \\mbox{ random on 53 bits, } e_a,e_b,e_c,e_d \\mbox{ on 2 bits}$" << endl;
427
cout <<"\\subsection*{$x=-y^t$}" << endl;
428
cout <<"$\\left| \\begin{array}{cc} a& -a\\\\ b & -b \\end{array} \\right|" << endl;
429
cout <<"\\hspace{1cm} a,c \\mbox{ random on 53 bits} $" << endl;
430
cout <<"\\subsection*{$x=-y+\\varepsilon^t$}" << endl;
431
cout <<"$\\left| \\begin{array}{cc} a+e_a& -a+e_c\\\\ b+e_b & -b+e_d \\end{array} \\right|" << endl;
432
cout <<"\\hspace{1cm} a,b \\mbox{ random on 53 bits, } e_a,e_b,e_c,e_d \\mbox{ on 2 bits}$" << endl;
433
cout <<"\\subsection*{$kU,lU$}" << endl;
434
cout <<"$\\left| \\begin{array}{cc} ka& la\\\\ kb & lb \\end{array} \\right|" << endl;
435
cout <<"\\hspace{1cm} a,b \\mbox{ random on 26 bits, } k,l \\mbox{ on 27 bits}$" << endl;
436
cout <<"\\subsection*{$kU,lU+\\varepsilon$}" << endl;
437
cout <<"$\\left| \\begin{array}{cc} ka+e_a& la+e_c\\\\ kb+e_b & lb+e_d \\end{array} \\right|" << endl;
438
cout <<"\\hspace{1cm} a,b \\mbox{ random on 53 bits, } k,l \\mbox{ on 27 bits, } e_a,e_b,e_c,e_d \\mbox{ on 2 bits}$" << endl;
439
cout <<"\\subsection*{$U,\\lfloor\\alpha U\\rfloor$}" << endl;
440
cout <<"$\\left| \\begin{array}{cc} a & c\\\\ \\lfloor \\alpha a\\rfloor & \\lfloor \\alpha c\\rfloor \\end{array} \\right|" << endl;
441
cout <<"\\hspace{1cm} a,c \\mbox{ random on 53 bits, } \\alpha \\mbox{ random in [-1,1]}$" << endl;
442
cout <<"\\subsection*{$U,\\lfloor\\alpha U\\rfloor^t$}" << endl;
443
cout <<"$\\left| \\begin{array}{cc} a& \\lfloor \\alpha a\\rfloor\\\\ b & \\lfloor \\alpha b\\rfloor \\end{array} \\right|" << endl;
444
cout <<"\\hspace{1cm} a,b \\mbox{ random on 53 bits, } \\alpha \\mbox{ random in [-1,1]}$" << endl;
445
cout <<"\\subsection*{$=$}" << endl;
446
cout <<"$\\left| \\begin{array}{cc} a& a\\\\ a & a \\end{array} \\right|" << endl;
447
cout <<"\\hspace{1cm} a \\mbox{ random on 53 bits} $" << endl;
448
cout <<"\\subsection*{$=+\\varepsilon$}" << endl;
449
cout <<"$\\left| \\begin{array}{cc} a+e_a& a+e_c\\\\ a+e_b & a+e_d \\end{array} \\right|" << endl;
450
cout <<"\\hspace{1cm} a \\mbox{ random on 53 bits, } e_a,e_b,e_c,e_d \\mbox{ on 2 bits}$" << endl;
452
cout <<"\\subsection*{Running time in $\\mu$s.}" << endl;
453
cout <<"\\begin{tabular}{|l||r|r|r|r|r|r|r|}\\hline" << endl;
454
cout <<"input&double&quadruple&leda-integer&leda-real¬-lazy&lazy&secure\\\\\\hline" << endl;
457
for ( input = 0 ; input < 11 ; input++ )
460
case 0: cout<< "random"; break;
461
case 1: cout<< "$x=-y$"; break;
462
case 2: cout<< "$x=-y+\\varepsilon$"; break;
463
case 3: cout<< "$x=-y^t$"; break;
464
case 4: cout<< "$x=-y+\\varepsilon^t$"; break;
465
case 5: cout<< "$kU,lU$"; break;
466
case 6: cout<< "$kU,lU+\\varepsilon$"; break;
467
case 7: cout<< "$U,\\lfloor\\alpha U\\rfloor$"; break;
468
case 8: cout<< "$U,\\lfloor\\alpha U\\rfloor^t$"; break;
469
case 9: cout<< "$=$"; break;
470
case 10: cout<< "$=+\\varepsilon$"; break;
473
for ( method = 0 ; method < 7 ; method++) time[method] = 0.0;
474
for (trial=0; trial < n ; trial++){
476
case 0: generate2(a,b,c,d,'r'); break;
477
case 1: generate2(a,b,c,d,'p'); break;
478
case 2: generate2(a,b,c,d,'p'); perturbation2(a,b,c,d); break;
479
case 3: generate2(a,b,c,d,'p'); transpose2(a,b,c,d); break;
480
case 4: generate2(a,b,c,d,'p'); perturbation2(a,b,c,d);transpose2(a,b,c,d); break;
481
case 5: generate2(a,b,c,d,'q'); break;
482
case 6: generate2(a,b,c,d,'q'); perturbation2(a,b,c,d); break;
483
case 7: generate2(a,b,c,d,'f'); break;
484
case 8: generate2(a,b,c,d,'f'); transpose2(a,b,c,d); break;
485
case 9: generate2(a,b,c,d,'c'); break;
486
case 10: generate2(a,b,c,d,'c'); perturbation2(a,b,c,d); break;
489
sign = leda_integer_det2x2(a,b,c,d);
491
sign = ABDPY__not_lazy_det2x2(a,b,c,d);
493
for ( method = 0 ; method < 7 ; method++) {
495
case 0 : FF2 = double_det2x2; break;
496
case 1 : FF2 = quadruple_det2x2; break;
498
case 2 : FF2 = leda_integer_det2x2; break;
499
case 3 : FF2 = leda_real_det2x2; N/=100; break;
501
case 2 : FF2 = 0; break;
502
case 3 : FF2 = 0; N/=100; break;
504
case 4 : FF2 = ABDPY__not_lazy_det2x2; N*=100; break;
505
case 5 : FF2 = ABDPY_det2x2; break;
507
case 6 : FF2 = ABDPY_secure_det2x2; break;
509
case 6 : FF2 = 0; break;
510
#endif // HAVE_SECURE
513
int sign_method = FF2(a,b,c,d);
514
if (sign != sign_method) error_2d[method]++;
516
if (method >= 2) { // compare leda_integer with all except method=0,1
518
if (method == 4) { // compare leda_integer with ABDPY__not_lazy
519
#endif // USE_SHEWCHUK
520
if (sign != sign_method) {
521
cerr << "ERROR det2x2: " << endl
522
<< input << " " << trial << endl
523
<< " double_det2x2 " << double_det2x2(a,b,c,d) << endl
524
<< " quadruple_det2x2 " << quadruple_det2x2(a,b,c,d) << endl
526
<< " leda_integer " << leda_integer_det2x2(a,b,c,d) << endl
527
<< " leda_real_det2x2 " << leda_real_det2x2(a,b,c,d) << endl
529
<< " ABDPY__not_lazy_det2x2 " << ABDPY__not_lazy_det2x2(a,b,c,d) << endl
530
<< " ABDPY_det2x2 " << ABDPY_det2x2(a,b,c,d) << endl
531
<< " ABDPY_secure_det2x2 " << ABDPY_secure_det2x2(a,b,c,d) << endl
533
cerr <<"/* maxima test code */" << endl;
534
cerr <<"a : ";exact_out(cerr,a); cerr << ";" << endl;
535
cerr <<"b : ";exact_out(cerr,b); cerr << ";" << endl;
536
cerr <<"c : ";exact_out(cerr,c); cerr << ";" << endl;
537
cerr <<"d : ";exact_out(cerr,d); cerr << ";" << endl;
538
cerr <<"det : a*d-b*c;" << endl;
539
cerr <<"/* end of maxima test code */" << endl;
540
cerr<<endl<<"send the above lines to olivier.devillers@sophia.inria.fr"<<endl;
546
for (j=0; j<N; j++) (void) FF2(a,b,c,d);
549
time[method] += t2 / (double) N;
552
if(!argc)for ( method = 0 ; method < 7 ; method++)cout<<"&"<<time[method]/n;
553
if(!argc)cout << "\\\\\\hline";
559
cout << "errors & " << endl;
560
for (size_t method = 0; method < 7; method++) {
561
cout << error_3d[method];
562
if (method != 6) cout << " & ";
563
else cout << " \\\\\\hline" << endl;
565
cout <<"\\end{tabular}" << endl;
567
cout <<"\\section*{Input 3D}" << endl;
568
cout <<"\\subsection*{random}" << endl;
569
cout <<"$\\left| \\begin{array}{ccc} a&d&g\\\\b&e&h\\\\c&f&i \\end{array} \\right|" << endl;
570
cout <<"\\hspace{1cm} a,b,c,d,e,f,g,h,i \\mbox{ random on 51 bits} $" << endl;
571
cout <<"\\subsection*{$x+y+z=0$}" << endl;
572
cout <<"$\\left| \\begin{array}{ccc} a&d&g\\\\b&e&h\\\\-a-b&-d-e&-g-h \\end{array} \\right|" << endl;
573
cout <<"\\hspace{1cm} a,b,d,e,g,h \\mbox{ random on 51 bits} $" << endl;
574
cout <<"\\subsection*{$x+y+z=0+\\varepsilon$}" << endl;
575
cout <<"$\\left| \\begin{array}{ccc} a+e_a&d+e_d&g+e_g\\\\b+e_b&e+e_e&h+e_h\\\\-a-b+e_c&-d-e+e_f&-g-h+e_i \\end{array} \\right|" << endl;
576
cout <<"\\hspace{1cm} a,b,d,e,g,h \\mbox{ random on 51 bits, }" << endl;
577
cout <<"e_a,e_b,e_c,e_d,e_e,e_f,e_g,e_h,e_i \\mbox{ on 2 bits}$" << endl;
578
cout <<"\\subsection*{$x+y+z=0^t$}" << endl;
579
cout <<"$\\left| \\begin{array}{ccc} a&d&-a-d\\\\b&e&-b-e\\\\c&f&-c-f \\end{array} \\right|" << endl;
580
cout <<"\\hspace{1cm} a,b,c,d,e,f \\mbox{ random on 51 bits} $" << endl;
581
cout <<"\\subsection*{$x+y+z=0+\\varepsilon^t$}" << endl;
582
cout <<"$\\left| \\begin{array}{ccc} a+e_a&d+e_d&-a-d+e_g\\\\b+e_b&e+e_e&-b-e+e_h\\\\c+e_c&f+e_f&-c-f+e_i \\end{array} \\right|" << endl;
583
cout <<"\\hspace{1cm} a,b,c,d,e,f \\mbox{ random on 51 bits, }" << endl;
584
cout <<"e_a,e_b,e_c,e_d,e_e,e_f,e_g,e_h,e_i \\mbox{ on 2 bits}$" << endl;
585
cout <<"\\subsection*{$kU,lV,mU+nV$}" << endl;
586
cout <<"$\\left| \\begin{array}{ccc} ka&ld&ma+nd\\\\kb&le&mb+ne\\\\kc&lf&mc+lf \\end{array} \\right|" << endl;
587
cout <<"\\hspace{1cm} a,b,c,d,e,f,m,n \\mbox{ random on 25 bits, } k,l \\mbox{ on 26 bits}$" << endl;
588
cout <<"\\subsection*{$kU,lV,mU+nV+\\varepsilon$}" << endl;
589
cout <<"$\\left| \\begin{array}{ccc} ka+e_a&ld+e_d&ma+nd+e_g\\\\kb+e_b&le+e_e&mb+ne+e_h\\\\kc+e_c&lf+e_f&mc+lf+e_i \\end{array} \\right|" << endl;
590
cout <<"\\hspace{1cm} \\begin{array}{l}a,b,c,d,e,f,m,n \\mbox{ random on 25 bits, } k,l \\mbox{ on 26 bits,}\\\\" << endl;
591
cout <<"e_a,e_b,e_c,e_d,e_e,e_f,e_g,e_h,e_i \\mbox{ on 2 bits} \\end{array}$" << endl;
592
cout <<"\\subsection*{$U,V,\\lfloor\\alpha U+\\beta V\\rfloor$}" << endl;
593
cout <<"$\\left| \\begin{array}{ccc} a&d&g\\\\b&e&h\\\\" << endl;
594
cout <<"\\lfloor \\alpha a+\\beta b\\rfloor&\\lfloor \\alpha d+\\beta e\\rfloor&" << endl;
595
cout <<"\\lfloor \\alpha g+\\beta h\\rfloor \\end{array} \\right|" << endl;
596
cout <<"\\hspace{1cm} a,b,d,e,g,h \\mbox{ random on 51 bits, } \\alpha,\\beta \\mbox{ random in [-1/2,1/2]}$" << endl;
597
cout <<"\\subsection*{$U,V,\\lfloor\\alpha U+\\beta V\\rfloor^t$}" << endl;
598
cout <<"$\\left| \\begin{array}{ccc} a&d&\\lfloor \\alpha a+\\beta d\\rfloor\\\\" << endl;
599
cout <<"b&e&\\lfloor \\alpha b+\\beta e\\rfloor\\\\" << endl;
600
cout <<"c&f&\\lfloor \\alpha c+\\beta f\\rfloor \\end{array} \\right|" << endl;
601
cout <<"\\hspace{1cm} a,b,c,d,e,f \\mbox{ random on 51 bits, } \\alpha,\\beta \\mbox{ random in [-1/2,1/2]}$" << endl;
602
cout <<"\\subsection*{$=$}" << endl;
603
cout <<"$\\left| \\begin{array}{ccc} a&a&a\\\\a&a&a\\\\a&a&a \\end{array} \\right|" << endl;
604
cout <<"\\hspace{1cm} a \\mbox{ random on 53 bits} $" << endl;
605
cout <<"\\subsection*{$=+\\varepsilon$}" << endl;
606
cout <<"$\\left| \\begin{array}{ccc} a+e_a&d+e_d&g+e_g\\\\b+e_b&e+e_e&h+e_h\\\\-a-b+e_c&-d-e+" << endl;
607
cout <<"e_f&-g-h+e_i \\end{array} \\right|" << endl;
608
cout <<"\\hspace{1cm} a \\mbox{ random on 53 bits, }" << endl;
609
cout <<"e_a,e_b,e_c,e_d,e_e,e_f,e_g,e_h,e_i \\mbox{ on 2 bits}$" << endl;
611
cout <<"\\subsection*{Running time in $\\mu$s.}" << endl;
612
cout <<"\\begin{tabular}{|l||r|r|r|r|r|r|}\\hline" << endl;
613
cout <<"input&double&leda-integer&leda-real¬-lazy&lazy&secure\\\\\\hline" << endl;
616
for ( input = 0 ; input < 11 ; input++ )
619
case 0: cout<< "random"; break;
620
case 1: cout<< "$x+y+z=0$"; break;
621
case 2: cout<< "$x+y+z=0+\\varepsilon$"; break;
622
case 3: cout<< "$x+y+z=0^t$"; break;
623
case 4: cout<< "$x+y+z=0+\\varepsilon^t$"; break;
624
case 5: cout<< "$kU,lV,mU+nV$"; break;
625
case 6: cout<< "$kU,lV,mU+nV+\\varepsilon$"; break;
626
case 7: cout<< "$U,V,\\lfloor\\alpha U+\\beta V\\rfloor$"; break;
627
case 8: cout<< "$U,V,\\lfloor\\alpha U+\\beta V\\rfloor^t$"; break;
628
case 9: cout<< "$=$"; break;
629
case 10: cout<< "$=+\\varepsilon$"; break;
631
for ( method = 0 ; method < 6 ; method++) time[method] = 0.0;
632
for (trial=0; trial < n ; trial++){
634
case 0: generate3(a,b,c,d,e,f,g,h,i,'r'); break;
635
case 1: generate3(a,b,c,d,e,f,g,h,i,'p'); break;
636
case 2: generate3(a,b,c,d,e,f,g,h,i,'p'); perturbation3(a,b,c,d,e,f,g,h,i); break;
637
case 3: generate3(a,b,c,d,e,f,g,h,i,'p'); transpose3(a,b,c,d,e,f,g,h,i); break;
638
case 4: generate3(a,b,c,d,e,f,g,h,i,'p'); perturbation3(a,b,c,d,e,f,g,h,i);transpose3(a,b,c,d,e,f,g,h,i); break;
639
case 5: generate3(a,b,c,d,e,f,g,h,i,'q'); break;
640
case 6: generate3(a,b,c,d,e,f,g,h,i,'q'); perturbation3(a,b,c,d,e,f,g,h,i); break;
641
case 7: generate3(a,b,c,d,e,f,g,h,i,'f'); break;
642
case 8: generate3(a,b,c,d,e,f,g,h,i,'f'); transpose3(a,b,c,d,e,f,g,h,i); break;
643
case 9: generate3(a,b,c,d,e,f,g,h,i,'c'); break;
644
case 10: generate3(a,b,c,d,e,f,g,h,i,'c'); perturbation3(a,b,c,d,e,f,g,h,i); break;
647
sign = leda_integer_det3x3(a,b,c,d,e,f,g,h,i);
649
sign = ABDPY__not_lazy_det3x3(a,b,c,d,e,f,g,h,i);
651
for ( method = 0 ; method < 6 ; method++) {
653
case 0 : FF3 = double_det3x3; break;
655
case 1 : FF3 = leda_integer_det3x3; break;
656
case 2 : FF3 = leda_real_det3x3; N/=100; break;
658
case 1: FF3 = 0; break;
659
case 2: FF3 = 0; N/=100; break;
661
case 3 : FF3 = ABDPY__not_lazy_det3x3; N*=100; break;
662
case 4 : FF3 = ABDPY_det3x3; break;
664
case 5 : FF3 = ABDPY_secure_det3x3; break;
666
case 5 : FF3 = 0; break;
667
#endif // HAVE_SECURE
670
int sign_method = FF3(a,b,c,d,e,f,g,h,i);
671
if (sign != sign_method) error_3d[method]++;
673
if (method != 0 && method != 2) { // compare leda_integer with all except method=0,2
675
if (method == 3) { // compare leda_integer with ABDPY__not_lazy
676
#endif // USE_SHEWCHUK
677
if (sign != sign_method) {
678
cerr << "ERROR det3x3:"
679
<< input << " " <<trial << endl
680
<< " double_det3x3 " << double_det3x3(a,b,c,d,e,f,g,h,i) << endl
682
<< " leda_integer_det3x3 " << leda_integer_det3x3(a,b,c,d,e,f,g,h,i) << endl
683
<< " leda_real_det3x3 " << leda_real_det3x3(a,b,c,d,e,f,g,h,i) << endl
685
<< " ABDPY__not_lazy_det3x3 " << ABDPY__not_lazy_det3x3(a,b,c,d,e,f,g,h,i) << endl
686
<< " ABDPY_det3x3 " << ABDPY_det3x3(a,b,c,d,e,f,g,h,i) << endl
687
<< " ABDPY_secure_det3x3 " << ABDPY_secure_det3x3(a,b,c,d,e,f,g,h,i) << endl;
688
cerr <<"/* maxima test code */" << endl;
689
cerr<<"a : ";exact_out(cerr,a); cerr << ";" << endl;
690
cerr<<"b : ";exact_out(cerr,b); cerr << ";" << endl;
691
cerr<<"c : ";exact_out(cerr,c); cerr << ";" << endl;
692
cerr<<"d : ";exact_out(cerr,d); cerr << ";" << endl;
693
cerr<<"e : ";exact_out(cerr,e); cerr << ";" << endl;
694
cerr<<"f : ";exact_out(cerr,f); cerr << ";" << endl;
695
cerr<<"g : ";exact_out(cerr,g); cerr << ";" << endl;
696
cerr<<"h : ";exact_out(cerr,h); cerr << ";" << endl;
697
cerr<<"i : ";exact_out(cerr,i); cerr << ";" << endl;
698
cerr << "m : matrix([a,b,c],[d,e,f],[g,h,i]);" << endl;
699
cerr << "det : determinant(m);" << endl;
700
cerr <<"/* end of maxima test code */" << endl;
701
cerr<<endl<<"send the above lines to olivier.devillers@sophia.inria.fr"<<endl;
708
(void) FF3(a,b,c,d,e,f,g,h,i);
711
time[method] += t2 / (double) N;
714
if(!argc)for ( method = 0 ; method < 6 ; method++) cout <<"&" << time[method]/n;
715
if(!argc)cout << "\\\\\\hline";
721
cout << "errors & " << endl;
722
for (size_t method = 0; method < 6; method++) {
723
cout << error_3d[method];
724
if (method != 5) cout << " & ";
725
else cout << " \\\\\\hline" << endl;
727
cout <<"\\end{tabular}" << endl;
730
cout <<"\\newcommand{\\etalchar}[1]{$^{#1}$}" << endl;
731
cout <<"\\begin{thebibliography}{ABD{\\etalchar{+}}94}" << endl;
732
cout <<"\\bibitem[ABD{\\etalchar{+}}94]{iter1}" << endl;
733
cout <<"F.~Avnaim, J.-D. Boissonnat, O.~Devillers, F.~Preparata, and M.~Yvinec." << endl;
734
cout <<"\\newblock Evaluating signs of determinants using single-precision arithmetic." << endl;
735
cout <<"\\newblock Research {Report} 2306, INRIA, BP93, 06902 Sophia-Antipolis, France, 1994." << endl;
736
cout <<"\\bibitem[ABD{\\etalchar{+}}95]{iter2}" << endl;
737
cout <<"F.~Avnaim, J.-D. Boissonnat, O.~Devillers, F.~Preparata, and M.~Yvinec." << endl;
738
cout <<"\\newblock Evaluating signs of determinants using single-precision arithmetic." << endl;
739
cout <<"\\newblock In {\\em Proc. 11th Annu. ACM Sympos. Comput. Geom.},1995." << endl;
740
cout <<"{\\tt http://www.inria.fr:/prisme/personnel/devillers/anglais/determinant.html}" << endl;
741
cout <<"\\bibitem[MBK{\\etalchar{+}}95]{leda}" << endl;
742
cout <<"K.~Mehlhorn, C.~Burnikel, J.~K\\\"onnemann, S.~N\\\"aher, S.~Schirra, and C.~Uhrig." << endl;
743
cout <<"\\newblock Exact Geometric Computation in LEDA" << endl;
744
cout <<"\\newblock In {\\em Proc. 11th Annu. ACM Sympos. Comput. Geom.},1995." << endl;
745
cout <<"{\\tt http://www.mpi-sb.mpg.de/guide/staff/uhrig/leda.html}" << endl;
746
cout <<"\\end{thebibliography}" << endl;
747
cout <<"\\end{document}" << endl;
749
return (n_error_2d == 0 && n_error_3d == 0) ? 0 : 1;
added
removed
nfem/basis/predicates_tst.cc
