1 /*
2  * Copyright (c) 1994 Cygnus Support.
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms are permitted
6  * provided that the above copyright notice and this paragraph are
7  * duplicated in all such forms and that any documentation,
8  * and/or other materials related to such
9  * distribution and use acknowledge that the software was developed
10  * at Cygnus Support, Inc.  Cygnus Support, Inc. may not be used to
11  * endorse or promote products derived from this software without
12  * specific prior written permission.
13  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
14  * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
15  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
16  */
17 
18 #include "test.h"
19 #include <errno.h>
20 
21 
22 int
randi(void)23 randi (void)
24 {
25   static int32_t next;
26   next = (next * 1103515245) + 12345;
27   return ((next >> 16) & 0xffff);
28 }
29 
30 #ifndef __FLOAT_WORD_ORDER__
31 #define __FLOAT_WORD_ORDER__ __BYTE_ORDER__
32 #endif
33 
randx(void)34 double randx (void)
35 {
36   double res;
37 
38   do
39   {
40     union {
41 	short parts[4];
42 	double res;
43       } u;
44 
45 #if __FLOAT_WORD_ORDER__ == __ORDER_LITTLE_ENDIAN__
46     u.parts[0] = randi();
47     u.parts[1] = randi();
48     u.parts[2] = randi();
49     u.parts[3] = randi();
50 #else
51     u.parts[3] = randi();
52     u.parts[2] = randi();
53     u.parts[1] = randi();
54     u.parts[0] = randi();
55 #endif
56     res = u.res;
57 
58   } while (!finite(res));
59 
60   return res ;
61 }
62 
63 /* Return a random double, but bias for numbers closer to 0 */
randy(void)64 double randy (void)
65 {
66   int pow;
67   double r= randx();
68   r = frexp(r, &pow);
69   return ldexp(r, randi() & 0x1f);
70 }
71 
72 void
test_frexp(void)73 test_frexp (void)
74 {
75   int i;
76   double r;
77   int t;
78 
79   float xf;
80   double gives;
81 
82   int pow;
83 
84 
85   /* Frexp of x return a and n, where a * 2**n == x, so test this with a
86      set of random numbers */
87   for (t = 0; t < 2; t++)
88   {
89     for (i = 0; i < 1000; i++)
90     {
91 
92       double x = randx();
93       line(i);
94       switch (t)
95       {
96       case 0:
97 	newfunc("frexp/ldexp");
98 	r = frexp(x, &pow);
99 	if (r > 1.0 || r < -1.0)
100 	{
101 	  /* Answer can never be > 1 or < 1 */
102 	  test_iok(0,1);
103 	}
104 
105 	gives = ldexp(r ,pow);
106 	test_mok(gives,x,62);
107 	break;
108       case 1:
109 	newfunc("frexpf/ldexpf");
110 	if (x > (double) FLT_MIN && x < (double) FLT_MAX)
111 	{
112 	  /* test floats too, but they have a smaller range so make sure x
113 	     isn't too big. Also x can get smaller than a float can
114 	     represent to make sure that doesn't happen too */
115 	  xf = x;
116 	  r = (double) frexpf(xf, &pow);
117 	  if (r > 1.0 || r < -1.0)
118 	  {
119 	    /* Answer can never be > 1 or < -1 */
120 	    test_iok(0,1);
121 	  }
122 
123 	  gives = (double) ldexpf(r ,pow);
124 	  test_mok(gives,x, 32);
125 
126 	}
127       }
128 
129     }
130 
131   }
132 
133   /* test a few numbers manually to make sure frexp/ldexp are not
134      testing as ok because both are broken */
135 
136   r = frexp(64.0, &i);
137 
138   test_mok(r, 0.5,64);
139   test_iok(i, 7);
140 
141   r = frexp(96.0, &i);
142 
143   test_mok(r, 0.75, 64);
144   test_iok(i, 7);
145 
146 }
147 
148 /* Test mod - this is given a real hammering by the strtod type
149    routines, here are some more tests.
150 
151    By definition
152 
153    modf = func(value, &iptr)
154 
155       (*iptr + modf) == value
156 
157    we test this
158 
159 */
160 void
test_mod(void)161 test_mod (void)
162 {
163   int i;
164 
165   newfunc("modf");
166 
167 
168   for (i = 0; i < 1000; i++)
169   {
170     double intpart;
171     double n;
172     line(i);
173     n  = randx();
174     if (finite(n) && n != 0.0 )
175     {
176       volatile double r = modf(n, &intpart);
177       line(i);
178       test_mok(intpart + r, n, 63);
179     }
180 
181   }
182   newfunc("modff");
183 
184   for (i = 0; i < 1000; i++)
185   {
186     float intpart;
187     double nd;
188     line(i);
189     nd  = randx() ;
190     if (fabs(nd) < (double) FLT_MAX && finitef(nd) && nd != 0.0)
191     {
192       volatile float n = nd;
193       volatile double r = (double) modff(n, &intpart);
194       line(i);
195       test_mok((double) intpart + r, (double) n, 32);
196     }
197   }
198 
199 
200 }
201 
202 /*
203 Test pow by multiplying logs
204 */
205 void
test_pow(void)206 test_pow (void)
207 {
208   unsigned int i;
209   newfunc("pow");
210 
211   for (i = 0; i < 1000; i++)
212   {
213     double n1;
214     double n2;
215     double res;
216     double shouldbe;
217 
218     line(i);
219     n1 = fabs(randy());
220     n2 = fabs(randy()/100.0);
221     res = pow(n1, n2);
222     shouldbe = exp(log(n1) * n2);
223     test_mok(shouldbe, res,55);
224   }
225 
226   newfunc("powf");
227 
228   for (i = 0; i < 1000; i++)
229   {
230     float n1;
231     float n2;
232     float res;
233     float shouldbe;
234 
235     errno = 0;
236 
237     line(i);
238     n1 = fabs(randy());
239     n2 = fabs(randy()/100.0);
240     res = powf(n1, n2);
241     shouldbe = expf(logf(n1) * n2);
242     if (!errno)
243       test_mok((double) shouldbe, (double) res,28);
244   }
245 
246 
247 
248 
249 }
250 
251 
252 
253 void
test_math2(void)254 test_math2 (void)
255 {
256   test_mod();
257   test_frexp();
258   test_pow();
259 }
260