1 | |
---|
2 | #include <NTL/ctools.h> |
---|
3 | |
---|
4 | |
---|
5 | #include <stdlib.h> |
---|
6 | #include <math.h> |
---|
7 | |
---|
8 | #if (defined(NTL_CXX_ONLY) && !defined(__cplusplus)) |
---|
9 | #error "CXX_ONLY flag set...must use C++ compiler" |
---|
10 | #endif |
---|
11 | |
---|
12 | /* |
---|
13 | * An IEEE double x is finite if and only if x - x == 0. |
---|
14 | * The function _ntl_IsFinite implements this logic; however, |
---|
15 | * it does not completely trust that an optimizing compiler |
---|
16 | * really implements this correctly, and so it goes out of its way to |
---|
17 | * confuse the compiler. For a good compiler that respects IEEE floating |
---|
18 | * point arithmetic, this is not be necessary, but it is better |
---|
19 | * to be a bit paranoid. |
---|
20 | * |
---|
21 | * Like the routine _ntl_ForceToMem below, this routine has the |
---|
22 | * side effect of forcing its argument into memory. |
---|
23 | */ |
---|
24 | |
---|
25 | double _ntl_IsFinite__local; |
---|
26 | double *_ntl_IsFinite__ptr1 = &_ntl_IsFinite__local; |
---|
27 | double *_ntl_IsFinite__ptr2 = &_ntl_IsFinite__local; |
---|
28 | double *_ntl_IsFinite__ptr3 = &_ntl_IsFinite__local; |
---|
29 | double *_ntl_IsFinite__ptr4 = &_ntl_IsFinite__local; |
---|
30 | |
---|
31 | long _ntl_IsFinite(double *p) |
---|
32 | { |
---|
33 | *_ntl_IsFinite__ptr1 = *p; |
---|
34 | *_ntl_IsFinite__ptr3 = (*_ntl_IsFinite__ptr2 - *p); |
---|
35 | if (*_ntl_IsFinite__ptr4 != 0.0) return 0; |
---|
36 | return 1; |
---|
37 | } |
---|
38 | |
---|
39 | |
---|
40 | /* |
---|
41 | * On machines with wide floating point registers, the routine _ntl_ForceToMem |
---|
42 | * is used to force a floating point double to a memory location. |
---|
43 | * This relies on "separate compilation" model, so that optimizing |
---|
44 | * compilers cannot "optimize away" the whole thing. |
---|
45 | */ |
---|
46 | |
---|
47 | |
---|
48 | void _ntl_ForceToMem(double *p) |
---|
49 | { } |
---|
50 | |
---|
51 | |
---|
52 | |
---|
53 | /* |
---|
54 | * The routine _ntl_ldexp(x, e) is like the standard ldexp(x, e) routine, |
---|
55 | * except that it takes a long exponent e, rather than an int exponenet. |
---|
56 | * Some care is taken to ensure reasonable overflow/undeflow behavior. |
---|
57 | * If the value of e does not fit into an int, then the result |
---|
58 | * is x*infinity or x*0, as appropriate. |
---|
59 | * Of course, this can only happen on platforms where long is wider |
---|
60 | * than int (e.g., most 64-bit platforms). |
---|
61 | * |
---|
62 | * We go out of our way to hide the fact that we are multiplying/dividing |
---|
63 | * by zero, so as to avoid unnecessary warnings, and to prevent |
---|
64 | * overly-agressive optimizing compilers from screwing things up. |
---|
65 | */ |
---|
66 | |
---|
67 | double _ntl_ldexp_zero = 0.0; |
---|
68 | |
---|
69 | double _ntl_ldexp(double x, long e) |
---|
70 | { |
---|
71 | if (e > NTL_MAX_INT) |
---|
72 | return x/_ntl_ldexp_zero; |
---|
73 | else if (e < NTL_MIN_INT) |
---|
74 | return x*_ntl_ldexp_zero; |
---|
75 | else |
---|
76 | return ldexp(x, ((int) e)); |
---|
77 | } |
---|
78 | |
---|
79 | |
---|
80 | void _ntl_abort(void) |
---|
81 | { |
---|
82 | _ntl_abort_cxx_callback(); |
---|
83 | abort(); |
---|
84 | } |
---|
85 | |
---|
86 | |
---|
87 | |
---|