1 /* Copyright (C) 1997-2013 Free Software Foundation, Inc.
2 This file is part of the GNU C Library.
3 Contributed by Geoffrey Keating <Geoff.Keating@anu.edu.au>, 1997.
5 The GNU C Library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
10 The GNU C Library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
15 You should have received a copy of the GNU Lesser General Public
16 License along with the GNU C Library; if not, see
17 <http://www.gnu.org/licenses/>. */
26 #define PRINT_ERRORS 0
32 #define mpbpl (CHAR_BIT * sizeof (mp_limb_t))
33 #define SZ (FRAC / mpbpl + 1)
34 typedef mp_limb_t mp1
[SZ
], mp2
[SZ
* 2];
36 /* These strings have exactly 100 hex digits in them. */
37 static const char sin1
[101] =
38 "d76aa47848677020c6e9e909c50f3c3289e511132f518b4def"
39 "b6ca5fd6c649bdfb0bd9ff1edcd4577655b5826a3d3b50c264";
40 static const char cos1
[101] =
41 "8a51407da8345c91c2466d976871bd29a2373a894f96c3b7f2"
42 "300240b760e6fa96a94430a52d0e9e43f3450e3b8ff99bc934";
43 static const char hexdig
[] = "0123456789abcdef";
46 print_mpn_hex (const mp_limb_t
*x
, unsigned size
)
50 const unsigned final
= (size
* 4 > SZ
* mpbpl
) ? SZ
* mpbpl
/ 4 : size
;
52 memset (value
, '0', size
);
54 for (i
= 0; i
< final
; i
++)
55 value
[size
-1-i
] = hexdig
[x
[i
* 4 / mpbpl
] >> (i
* 4) % mpbpl
& 0xf];
58 fputs (value
, stdout
);
62 sincosx_mpn (mp1 si
, mp1 co
, mp1 xx
, mp1 ix
)
71 memset (si
, 0, sizeof (mp1
));
72 memset (co
, 0, sizeof (mp1
));
74 memcpy (x
, xx
, sizeof (mp1
));
77 mpn_sub_n (x
, xx
, ix
, SZ
);
79 for (i
= 0; i
< 1 << N
; i
++)
81 #define add_shift_mulh(d,x,s1,s2,sh,n) \
82 /* d = (n ? -1 : 1) * (s1 + (s2>>sh)) * x / (1>>N); */ \
86 assert (sh < mpbpl); \
87 mpn_lshift (tmp, s1, SZ, sh); \
88 chk = (n ? mpn_sub_n : mpn_add_n)(tmp,tmp,s2+FRAC/mpbpl,SZ); \
90 chk = (n ? mpn_sub_n : mpn_add_n)(tmp,s1,s2+FRAC/mpbpl,SZ); \
91 /* assert(chk == 0); */ \
92 mpn_mul_n(d,tmp,x,SZ); \
94 mpn_mul_n(d,s1,x,SZ); \
95 /* assert(d[SZ*2-1] == 0); */ \
96 assert(N+sh < mpbpl); \
97 if (N+sh > 0) mpn_rshift(d,d,2*SZ,N+sh); \
99 #define summ(d,ss,s,n) \
100 /* d = ss +/- (s[0]+2*s[1]+2*s[2]+s[3])/6; */ \
102 chk = mpn_add_n(tmp,s[1]+FRAC/mpbpl,s[2]+FRAC/mpbpl,SZ); \
103 mpn_lshift(tmp,tmp,SZ,1); \
104 chk |= mpn_add_n(tmp,tmp,s[0]+FRAC/mpbpl,SZ); \
105 chk |= mpn_add_n(tmp,tmp,s[3]+FRAC/mpbpl,SZ); \
106 round = mpn_divmod_1(tmp,tmp,SZ,6); \
107 /* chk |= mpn_add_1(tmp,tmp,SZ, (round > 3) ); */ \
108 chk |= (n ? mpn_sub_n : mpn_add_n)(d,ss,tmp,SZ); \
109 /* assert(chk == 0); */ \
112 add_shift_mulh (s
[0], x
, co
, NULL
, 0, 0); /* s0 = h * c; */
113 add_shift_mulh (c
[0], x
, si
, NULL
, 0, 0); /* c0 = h * s; */
114 add_shift_mulh (s
[1], x
, co
, c
[0], 1, 1); /* s1 = h * (c - c0/2); */
115 add_shift_mulh (c
[1], x
, si
, s
[0], 1, 0); /* c1 = h * (s + s0/2); */
116 add_shift_mulh (s
[2], x
, co
, c
[1], 1, 1); /* s2 = h * (c - c1/2); */
117 add_shift_mulh (c
[2], x
, si
, s
[1], 1, 0); /* c2 = h * (s + s1/2); */
118 add_shift_mulh (s
[3], x
, co
, c
[2], 0, 1); /* s3 = h * (c - c2); */
119 add_shift_mulh (c
[3], x
, si
, s
[2], 0, 0); /* c3 = h * (s + s2); */
120 summ (si
, si
, s
, 0); /* s = s + (s0+2*s1+2*s2+s3)/6; */
121 summ (co
, co
, c
, 1); /* c = c - (c0+2*c1+2*c2+c3)/6; */
123 #undef add_shift_mulh
128 mpn_bitsize (const mp_limb_t
*SRC_PTR
, mp_size_t SIZE
)
131 for (i
= SIZE
- 1; i
> 0; i
--)
134 for (j
= mpbpl
- 1; j
>= 0; j
--)
135 if ((SRC_PTR
[i
] & (mp_limb_t
)1 << j
) != 0)
138 return i
* mpbpl
+ j
;
144 mp1 si
, co
, x
, ox
, xt
, s2
, c2
, s3
, c3
;
146 int sin_errors
= 0, cos_errors
= 0;
147 int sin_failures
= 0, cos_failures
= 0;
148 mp1 sin_maxerror
, cos_maxerror
;
149 int sin_maxerror_s
= 0, cos_maxerror_s
= 0;
150 const double sf
= pow (2, mpbpl
);
152 /* assert(mpbpl == mp_bits_per_limb); */
153 assert(FRAC
/ mpbpl
* mpbpl
== FRAC
);
155 memset (sin_maxerror
, 0, sizeof (mp1
));
156 memset (cos_maxerror
, 0, sizeof (mp1
));
157 memset (xt
, 0, sizeof (mp1
));
158 xt
[(FRAC
- N2
) / mpbpl
] = (mp_limb_t
)1 << (FRAC
- N2
) % mpbpl
;
160 for (i
= 0; i
< 1 << N2
; i
++)
162 int s2s
, s3s
, c2s
, c3s
, j
;
165 mpn_mul_1 (x
, xt
, SZ
, i
);
166 sincosx_mpn (si
, co
, x
, i
== 0 ? NULL
: ox
);
167 memcpy (ox
, x
, sizeof (mp1
));
168 ds2
= sin (i
/ (double) (1 << N2
));
169 dc2
= cos (i
/ (double) (1 << N2
));
170 for (j
= SZ
-1; j
>= 0; j
--)
172 s2
[j
] = (mp_limb_t
) ds2
;
173 ds2
= (ds2
- s2
[j
]) * sf
;
174 c2
[j
] = (mp_limb_t
) dc2
;
175 dc2
= (dc2
- c2
[j
]) * sf
;
177 if (mpn_cmp (si
, s2
, SZ
) >= 0)
178 mpn_sub_n (s3
, si
, s2
, SZ
);
180 mpn_sub_n (s3
, s2
, si
, SZ
);
181 if (mpn_cmp (co
, c2
, SZ
) >= 0)
182 mpn_sub_n (c3
, co
, c2
, SZ
);
184 mpn_sub_n (c3
, c2
, co
, SZ
);
186 s2s
= mpn_bitsize (s2
, SZ
);
187 s3s
= mpn_bitsize (s3
, SZ
);
188 c2s
= mpn_bitsize (c2
, SZ
);
189 c3s
= mpn_bitsize (c3
, SZ
);
190 if ((s3s
>= 0 && s2s
- s3s
< 54)
191 || (c3s
>= 0 && c2s
- c3s
< 54)
195 printf ("%06x ", i
* (0x100000 / (1 << N2
)));
196 print_mpn_hex(si
, (FRAC
/ 4) + 1);
198 print_mpn_hex (co
, (FRAC
/ 4) + 1);
201 print_mpn_hex (s2
, (FRAC
/ 4) + 1);
203 print_mpn_hex (c2
, (FRAC
/ 4) + 1);
206 s3s
>= 0 && s2s
-s3s
< 54 ? s2s
- s3s
== 53 ? 'e' : 'F' : 'P',
207 c3s
>= 0 && c2s
-c3s
< 54 ? c2s
- c3s
== 53 ? 'e' : 'F' : 'P');
208 print_mpn_hex (s3
, (FRAC
/ 4) + 1);
210 print_mpn_hex (c3
, (FRAC
/ 4) + 1);
213 sin_errors
+= s2s
- s3s
== 53;
214 cos_errors
+= c2s
- c3s
== 53;
215 sin_failures
+= s2s
- s3s
< 53;
216 cos_failures
+= c2s
- c3s
< 53;
218 if (s3s
>= sin_maxerror_s
219 && mpn_cmp (s3
, sin_maxerror
, SZ
) > 0)
221 memcpy (sin_maxerror
, s3
, sizeof (mp1
));
222 sin_maxerror_s
= s3s
;
224 if (c3s
>= cos_maxerror_s
225 && mpn_cmp (c3
, cos_maxerror
, SZ
) > 0)
227 memcpy (cos_maxerror
, c3
, sizeof (mp1
));
228 cos_maxerror_s
= c3s
;
232 /* Check Range-Kutta against precomputed values of sin(1) and cos(1). */
233 memset (x
, 0, sizeof (mp1
));
234 x
[FRAC
/ mpbpl
] = (mp_limb_t
)1 << FRAC
% mpbpl
;
235 sincosx_mpn (si
, co
, x
, ox
);
237 memset (s2
, 0, sizeof (mp1
));
238 memset (c2
, 0, sizeof (mp1
));
239 for (i
= 0; i
< 100 && i
< FRAC
/ 4; i
++)
241 s2
[(FRAC
- i
* 4 - 4) / mpbpl
] |= ((mp_limb_t
) (strchr (hexdig
, sin1
[i
])
243 << (FRAC
- i
* 4 - 4) % mpbpl
);
244 c2
[(FRAC
- i
* 4 - 4) / mpbpl
] |= ((mp_limb_t
) (strchr (hexdig
, cos1
[i
])
246 << (FRAC
- i
* 4 - 4) % mpbpl
);
249 if (mpn_cmp (si
, s2
, SZ
) >= 0)
250 mpn_sub_n (s3
, si
, s2
, SZ
);
252 mpn_sub_n (s3
, s2
, si
, SZ
);
253 if (mpn_cmp (co
, c2
, SZ
) >= 0)
254 mpn_sub_n (c3
, co
, c2
, SZ
);
256 mpn_sub_n (c3
, c2
, co
, SZ
);
259 printf ("%d failures; %d errors; error rate %0.2f%%\n",
260 sin_failures
, sin_errors
, sin_errors
* 100.0 / (double) (1 << N2
));
261 fputs ("maximum error: ", stdout
);
262 print_mpn_hex (sin_maxerror
, (FRAC
/ 4) + 1);
263 fputs ("\nerror in sin(1): ", stdout
);
264 print_mpn_hex (s3
, (FRAC
/ 4) + 1);
266 fputs ("\n\ncos:\n", stdout
);
267 printf ("%d failures; %d errors; error rate %0.2f%%\n",
268 cos_failures
, cos_errors
, cos_errors
* 100.0 / (double) (1 << N2
));
269 fputs ("maximum error: ", stdout
);
270 print_mpn_hex (cos_maxerror
, (FRAC
/ 4) + 1);
271 fputs ("\nerror in cos(1): ", stdout
);
272 print_mpn_hex (c3
, (FRAC
/ 4) + 1);
275 return (sin_failures
== 0 && cos_failures
== 0) ? 0 : 1;