summaryrefslogtreecommitdiff
path: root/tests/RRTest.c
blob: 9c09b28c97ddbe41ff6b4afa6cacad2c589f582a (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
/* Example file to test mpfr addition against NTL.
   Usage:
   0) compile this file with NTL
   1) compile tadd.c with -DCHECK_EXTERNAL
   2) ./tadd | egrep -v 'Seed|Inf|NaN' > /tmp/log
      (Warning, this produces a large file.)
   3) ./RRTest < /tmp/add.log
*/

#include <NTL/mat_RR.h>

NTL_CLIENT

void
ReadBinary (ZZ &x)
{
  int s = 1;
  ZZ b, y;

  cin >> b;

  if (b < 0)
    {
      s = -1;
      b = -b;
    }

  x = 0;
  y = 1;
  while (b != 0)
    {
      x += (b % 10) * y;
      y *= 2;
      b /= 10;
    }
  if (s < 0)
    x = -x;
}

long
ReadRR (RR &a)
{
  long p;
  ZZ x;
  long e;

  cin >> p;
  ReadBinary (x);
  cin >> e;
  MakeRRPrec (a, x, e, p);
  return p;
}

void
Output (RR a, long p)
{
  cout << a.mantissa() << "*2^(" << a.exponent() << ") [" << p << "]" << endl;
}

// ulp difference between a and b
long
ulp (RR a, RR b, long p)
{
  ZZ ma, mb;
  long ea, eb;

  ma = a.x;
  ea = a.e;
  while (NumBits (ma) < p)
    {
      ma *= 2;
      ea --;
    }
  mb = b.x;
  eb = b.e;
  while (NumBits (mb) < p)
    {
      mb *= 2;
      eb --;
    }
  if (ea != eb) abort ();
  return to_long (ma - mb);
}

// #define TWO_ARGS /* for functions of two arguments like add, sub, pow */

int
main (void)
{
  RR a, b, c, d;
  long line = 0, errors = 0;
  long pa, pb, pc;

  while (!feof(stdin))
    {
      if (++line % 10 == 0)
        cout << "line " << line << endl;
      pb = ReadRR (b);
#ifdef TWO_ARGS
      pc = ReadRR (c);
#endif
      pa = ReadRR (a);
      RR::SetPrecision(pa);
      cos (d, b
#ifdef TWO_ARGS
           , c
#endif
           );
      if (d != a)
        {
          cerr << "error at line " << line << endl;
          cerr << "b="; Output(b, pb);
#ifdef TWO_ARGS
          cerr << " c="; Output(c, pc);
#endif
          cerr << "expected "; Output(a, pa);
          cerr << "got      "; Output(d, pa);
          cerr << "difference is " << ulp (a, d, pa) << " ulps" << endl;
          cerr << ++errors << " errors" << endl;
        }
    }
}