1/*	$NetBSD: systfloat.c,v 1.7 2004/04/15 19:01:57 matt Exp $	*/
2
3/* This is a derivative work. */
4
5/*-
6 * Copyright (c) 2001 The NetBSD Foundation, Inc.
7 * All rights reserved.
8 *
9 * This code is derived from software contributed to The NetBSD Foundation
10 * by Ross Harvey.
11 *
12 * Redistribution and use in source and binary forms, with or without
13 * modification, are permitted provided that the following conditions
14 * are met:
15 * 1. Redistributions of source code must retain the above copyright
16 *    notice, this list of conditions and the following disclaimer.
17 * 2. Redistributions in binary form must reproduce the above copyright
18 *    notice, this list of conditions and the following disclaimer in the
19 *    documentation and/or other materials provided with the distribution.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
22 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
23 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
24 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
25 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 */
33
34/*
35===============================================================================
36
37This C source file is part of TestFloat, Release 2a, a package of programs
38for testing the correctness of floating-point arithmetic complying to the
39IEC/IEEE Standard for Floating-Point.
40
41Written by John R. Hauser.  More information is available through the Web
42page `http://HTTP.CS.Berkeley.EDU/~jhauser/arithmetic/TestFloat.html'.
43
44THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort
45has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT
46TIMES RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO
47PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY
48AND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE.
49
50Derivative works are acceptable, even for commercial purposes, so long as
51(1) they include prominent notice that the work is derivative, and (2) they
52include prominent notice akin to these four paragraphs for those parts of
53this code that are retained.
54
55===============================================================================
56*/
57
58#include <sys/cdefs.h>
59#ifndef __lint
60__RCSID("$NetBSD: systfloat.c,v 1.7 2004/04/15 19:01:57 matt Exp $");
61#endif
62
63#include <math.h>
64#include <ieeefp.h>
65#include "milieu.h"
66#include "softfloat.h"
67#include "systfloat.h"
68#include "systflags.h"
69#include "systmodes.h"
70
71typedef union {
72    float32 f32;
73    float f;
74} union32;
75typedef union {
76    float64 f64;
77    double d;
78} union64;
79#if defined( FLOATX80 ) && defined( LONG_DOUBLE_IS_FLOATX80 )
80typedef union {
81    floatx80 fx80;
82    long double ld;
83} unionx80;
84#endif
85#if defined( FLOAT128 ) && defined( LONG_DOUBLE_IS_FLOAT128 )
86typedef union {
87    float128 f128;
88    long double ld;
89} union128;
90#endif
91
92fp_except
93syst_float_flags_clear(void)
94{
95    return fpsetsticky(0)
96	& (FP_X_IMP | FP_X_UFL | FP_X_OFL | FP_X_DZ | FP_X_INV);
97}
98
99void
100syst_float_set_rounding_mode(fp_rnd direction)
101{
102    fpsetround(direction);
103    fpsetmask(0);
104}
105
106float32 syst_int32_to_float32( int32 a )
107{
108    const union32 uz = { .f = a };
109
110    return uz.f32;
111
112}
113
114float64 syst_int32_to_float64( int32 a )
115{
116    const union64 uz = { .d = a };
117
118    return uz.f64;
119
120}
121
122#if defined( FLOATX80 ) && defined( LONG_DOUBLE_IS_FLOATX80 )
123
124floatx80 syst_int32_to_floatx80( int32 a )
125{
126    const unionx80 uz = { .ld = a };
127
128    return uz.fx80;
129
130}
131
132#endif
133
134#if defined( FLOAT128 ) && defined( LONG_DOUBLE_IS_FLOAT128 )
135
136float128 syst_int32_to_float128( int32 a )
137{
138    const union128 uz = { .ld = a };
139
140    return uz.f128;
141
142}
143
144#endif
145
146#ifdef BITS64
147
148float32 syst_int64_to_float32( int64 a )
149{
150    const union32 uz = { .f = a };
151
152    return uz.f32;
153}
154
155float64 syst_int64_to_float64( int64 a )
156{
157    const union64 uz = { .d = a };
158
159    return uz.f64;
160}
161
162#if defined( FLOATX80 ) && defined( LONG_DOUBLE_IS_FLOATX80 )
163
164floatx80 syst_int64_to_floatx80( int64 a )
165{
166    const unionx80 uz = { .ld = a };
167
168    return uz.fx80;
169}
170
171#endif
172
173#if defined( FLOAT128 ) && defined( LONG_DOUBLE_IS_FLOAT128 )
174
175float128 syst_int64_to_float128( int64 a )
176{
177    const union128 uz = { .ld = a };
178
179    return uz.f128;
180}
181
182#endif
183
184#endif
185
186int32 syst_float32_to_int32_round_to_zero( float32 a )
187{
188    const union32 uz = { .f32 = a };
189
190    return uz.f;
191
192}
193
194#ifdef BITS64
195
196int64 syst_float32_to_int64_round_to_zero( float32 a )
197{
198    const union32 uz = { .f32 = a };
199
200    return uz.f;
201}
202
203#endif
204
205float64 syst_float32_to_float64( float32 a )
206{
207    const union32 ua = { .f32 = a };
208    union64 uz;
209
210    uz.d = ua.f;
211    return uz.f64;
212
213}
214
215#if defined( FLOATX80 ) && defined( LONG_DOUBLE_IS_FLOATX80 )
216
217floatx80 syst_float32_to_floatx80( float32 a )
218{
219    const union32 ua = { .f32 = a };
220    unionx80 uz;
221
222    uz.ld = ua.f;
223    return uz.fx80;
224}
225
226#endif
227
228#if defined( FLOAT128 ) && defined( LONG_DOUBLE_IS_FLOAT128 )
229
230float128 syst_float32_to_float128( float32 a )
231{
232    const union32 ua = { .f32 = a };
233    union128 ub;
234
235    ub.ld = ua.f;
236    return ub.f128;
237}
238
239#endif
240
241float32 syst_float32_add( float32 a, float32 b )
242{
243    const union32 ua = { .f32 = a }, ub = { .f32 = b };
244    union32 uz;
245
246    uz.f = ua.f + ub.f;
247    return uz.f32;
248}
249
250float32 syst_float32_sub( float32 a, float32 b )
251{
252    const union32 ua = { .f32 = a }, ub = { .f32 = b };
253    union32 uz;
254
255    uz.f = ua.f - ub.f;
256    return uz.f32;
257}
258
259float32 syst_float32_mul( float32 a, float32 b )
260{
261    const union32 ua = { .f32 = a }, ub = { .f32 = b };
262    union32 uz;
263
264    uz.f = ua.f * ub.f;
265    return uz.f32;
266}
267
268float32 syst_float32_div( float32 a, float32 b )
269{
270    const union32 ua = { .f32 = a }, ub = { .f32 = b };
271    union32 uz;
272
273    uz.f = ua.f / ub.f;
274    return uz.f32;
275}
276
277flag syst_float32_eq( float32 a, float32 b )
278{
279    const union32 ua = { .f32 = a }, ub = { .f32 = b };
280
281    return ua.f == ub.f;
282}
283
284flag syst_float32_le( float32 a, float32 b )
285{
286    const union32 ua = { .f32 = a }, ub = { .f32 = b };
287
288    return ua.f <= ub.f;
289}
290
291flag syst_float32_lt( float32 a, float32 b )
292{
293    const union32 ua = { .f32 = a }, ub = { .f32 = b };
294
295    return ua.f < ub.f;
296}
297
298int32 syst_float64_to_int32_round_to_zero( float64 a )
299{
300    const union64 uz = { .f64 = a };
301
302    return uz.d;
303}
304
305#ifdef BITS64
306
307int64 syst_float64_to_int64_round_to_zero( float64 a )
308{
309    const union64 uz = { .f64 = a };
310
311    return uz.d;
312}
313
314#endif
315
316float32 syst_float64_to_float32( float64 a )
317{
318    const union64 ua = { .f64 = a };
319    union32 uz;
320
321    uz.f = ua.d;
322    return uz.f32;
323}
324
325#if defined( FLOATX80 ) && defined( LONG_DOUBLE_IS_FLOATX80 )
326
327floatx80 syst_float64_to_floatx80( float64 a )
328{
329    const union64 ua = { .f64 = a };
330    unionx80 u;
331
332    u.ld = ua.d;
333    return u.fx80;
334}
335
336#endif
337
338#if defined( FLOAT128 ) && defined( LONG_DOUBLE_IS_FLOAT128 )
339
340float128 syst_float64_to_float128( float64 a )
341{
342    const union64 ua = { .f64 = a };
343    union128 uz;
344
345    uz.ld = ua.d;
346    return uz.f128;
347}
348
349#endif
350
351float64 syst_float64_add( float64 a, float64 b )
352{
353    const union64 ua = { .f64 = a }, ub = { .f64 = b };
354    union64 uz;
355
356    uz.d = ua.d + ub.d;
357    return uz.f64;
358}
359
360float64 syst_float64_sub( float64 a, float64 b )
361{
362    const union64 ua = { .f64 = a }, ub = { .f64 = b };
363    union64 uz;
364
365    uz.d = ua.d - ub.d;
366    return uz.f64;
367}
368
369float64 syst_float64_mul( float64 a, float64 b )
370{
371    const union64 ua = { .f64 = a }, ub = { .f64 = b };
372    union64 uz;
373
374    uz.d = ua.d * ub.d;
375    return uz.f64;
376}
377
378float64 syst_float64_div( float64 a, float64 b )
379{
380    const union64 ua = { .f64 = a }, ub = { .f64 = b };
381    union64 uz;
382
383    uz.d = ua.d / ub.d;
384    return uz.f64;
385}
386
387float64 syst_float64_sqrt( float64 a )
388{
389    const union64 ua = { .f64 = a };
390    union64 uz;
391
392    uz.d = sqrt(ua.d);
393    return uz.f64;
394}
395
396flag syst_float64_eq( float64 a, float64 b )
397{
398    const union64 ua = { .f64 = a }, ub = { .f64 = b };
399
400    return ua.d == ub.d;
401}
402
403flag syst_float64_le( float64 a, float64 b )
404{
405    const union64 ua = { .f64 = a }, ub = { .f64 = b };
406
407    return ua.d <= ub.d;
408}
409
410flag syst_float64_lt( float64 a, float64 b )
411{
412    const union64 ua = { .f64 = a }, ub = { .f64 = b };
413
414    return ua.d < ub.d;
415}
416
417#if defined( FLOATX80 ) && defined( LONG_DOUBLE_IS_FLOATX80 )
418
419int32 syst_floatx80_to_int32_round_to_zero( floatx80 a )
420{
421    const unionx80 uz = { .fx80 = a };
422
423    return uz.ld;
424}
425
426#ifdef BITS64
427
428int64 syst_floatx80_to_int64_round_to_zero( floatx80 a )
429{
430    const unionx80 uz = { .fx80 = a };
431
432    return uz.ld;
433}
434
435#endif
436
437float32 syst_floatx80_to_float32( floatx80 a )
438{
439    const unionx80 ua = { .fx80 = a };
440    union32 uz;
441
442    uz.f = ua.ld;
443    return uz.f32;
444}
445
446float64 syst_floatx80_to_float64( floatx80 a )
447{
448    const unionx80 ua = { .fx80 = a };
449    union64 uz;
450
451    uz.d = ua.ld;
452    return uz.f64;
453}
454
455floatx80 syst_floatx80_add( floatx80 a, floatx80 b )
456{
457    const unionx80 ua = { .fx80 = a }, ub = { .fx80 = b };
458    unionx80 uz;
459
460    uz.ld = ua.ld + ub.ld;
461    return uz.fx80;
462}
463
464floatx80 syst_floatx80_sub( floatx80 a, floatx80 b )
465{
466    const unionx80 ua = { .fx80 = a }, ub = { .fx80 = b };
467    unionx80 uz;
468
469    uz.ld = ua.ld - ub.ld;
470    return uz.fx80;
471}
472
473floatx80 syst_floatx80_mul( floatx80 a, floatx80 b )
474{
475    const unionx80 ua = { .fx80 = a }, ub = { .fx80 = b };
476    unionx80 uz;
477
478    uz.ld = ua.ld * ub.ld;
479    return uz.fx80;
480}
481
482floatx80 syst_floatx80_div( floatx80 a, floatx80 b )
483{
484    const unionx80 ua = { .fx80 = a }, ub = { .fx80 = b };
485    unionx80 uz;
486
487    uz.ld = ua.ld / ub.ld;
488    return uz.fx80;
489}
490
491flag syst_floatx80_eq( floatx80 a, floatx80 b )
492{
493    const unionx80 ua = { .fx80 = a }, ub = { .fx80 = b };
494
495    return ua.ld == ub.ld;
496}
497
498flag syst_floatx80_le( floatx80 a, floatx80 b )
499{
500    const unionx80 ua = { .fx80 = a }, ub = { .fx80 = b };
501
502    return ua.ld <= ub.ld;
503}
504
505flag syst_floatx80_lt( floatx80 a, floatx80 b )
506{
507    const unionx80 ua = { .fx80 = a }, ub = { .fx80 = b };
508
509    return ua.ld < ub.ld;
510}
511
512#endif
513
514#if defined( FLOAT128 ) && defined( LONG_DOUBLE_IS_FLOAT128 )
515
516int32 syst_float128_to_int32_round_to_zero( float128 a )
517{
518    const union128 ua = { .f128 = a };
519
520    return ua.ld;
521}
522
523#ifdef BITS64
524
525int64 syst_float128_to_int64_round_to_zero( float128 a )
526{
527    const union128 ua = { .f128 = a };
528
529    return ua.ld;
530}
531
532#endif
533
534float32 syst_float128_to_float32( float128 a )
535{
536    const union128 ua = { .f128 = a };
537    union32 uz;
538
539    uz.f = ua.ld;
540    return uz.f32;
541
542}
543
544float64 syst_float128_to_float64( float128 a )
545{
546    const union128 ua = { .f128 = a };
547    union64 uz;
548
549    uz.d = ua.ld;
550    return uz.f64;
551}
552
553float128 syst_float128_add( float128 a, float128 b )
554{
555    const union128 ua = { .f128 = a }, ub = { .f128 = b };
556    union128 uz;
557
558    uz.ld = ua.ld + ub.ld;
559    return uz.f128;
560
561}
562
563float128 syst_float128_sub( float128 a, float128 b )
564{
565    const union128 ua = { .f128 = a }, ub = { .f128 = b };
566    union128 uz;
567
568    uz.ld = ua.ld - ub.ld;
569    return uz.f128;
570}
571
572float128 syst_float128_mul( float128 a, float128 b )
573{
574    const union128 ua = { .f128 = a }, ub = { .f128 = b };
575    union128 uz;
576
577    uz.ld = ua.ld * ub.ld;
578    return uz.f128;
579}
580
581float128 syst_float128_div( float128 a, float128 b )
582{
583    const union128 ua = { .f128 = a }, ub = { .f128 = b };
584    union128 uz;
585
586    uz.ld = ua.ld / ub.ld;
587    return uz.f128;
588}
589
590flag syst_float128_eq( float128 a, float128 b )
591{
592    const union128 ua = { .f128 = a }, ub = { .f128 = b };
593
594    return ua.ld == ub.ld;
595}
596
597flag syst_float128_le( float128 a, float128 b )
598{
599    const union128 ua = { .f128 = a }, ub = { .f128 = b };
600
601    return ua.ld <= ub.ld;
602}
603
604flag syst_float128_lt( float128 a, float128 b )
605{
606    const union128 ua = { .f128 = a }, ub = { .f128 = b };
607
608    return ua.ld < ub.ld;
609}
610
611#endif
612