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
|
This is work in progress. Before continuing, please read
http://websympa.loria.fr/wwsympa/arc/mpfr/2011-01/msg00039.html.
/* mpfr_out_raw -- output a floating-point number to binary portable format
Copyright 2011-2017 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
The GNU MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.
The GNU MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
#include "mpfr-impl.h"
/* format for out_raw:
A mpfr_t is represented by up to 3 fields, each one is represented by a
sequence of 32-bit words. 32-bit words are stored as 4 bytes in little
endian format:
(a) a field for sign, precision and other bit-fields:
- the sign (1 bit)
- 2 bits for NaN, Inf, 0, normal numbers:
00: 0
01: normal
10: Inf
11: NaN
- 1 bit for the exponent encoding (exp_enc) for normal numbers,
otherwise 0 is stored here (reserved for future extensions)
- 1 bit for the precision encoding (prec_enc)
If prec_enc=0, the remaining 27 bits encode the precision (< 2^27)
If prec_enc=1, the precision is stored in the following 27 bits
(high part) and then 32 bits (low part). Thus the maximal precision
is 59 bits.
(b) (optional) a field for the exponent:
- if the number is NaN, Inf, 0, this field is empty
- if exp_enc=0, this field contains one 32-bit (signed) word encoding
the exponent
- if exp_enc=1, a first 32-bit word encodes a positive integer m,
and the following m 32-bit words encode the exponent (in 2-complement
representation, with least significant words first)
(c) (optional) a field for the significand:
- if the number is NaN, Inf, 0, this field is empty
- otherwise, let p = ceil(prec/32), the significand is represented
by p consecutive 32-bit words (least significant words first).
Thus on a little-endian machine the significand can be directly
copied using memcopy.
Examples:
- a normal binary32 IEEE-754 number uses 96 bits: 32 for (a), 32 for (b),
and 32 for (c);
- a normal binary64 IEEE-754 number uses 128 bits: 32 for (a), 32 for (b),
and 64 for (c) (idem for a significand of 64 bits, as in Intel x86
double-extended format);
- a normal binary128 IEEE-754 number uses 192 bits: 32 for (a), 32 for (b),
and 128 for (c).
*/
size_t
mpfr_out_raw (FILE *stream, mpfr_srcptr x)
{
size_t n; /* number of bytes of the output */
mpfr_prec_t prec = MPFR_PREC(x);
int prec_enc, exp_enc = 0;
unsigned char *s, *t;
MPFR_ASSERTN (CHAR_BIT == 8);
prec_enc = (prec >= 134217728);
n = 4 + 4 * prec_enc;
if (MPFR_LIKELY (! MPFR_IS_SINGULAR(x)))
{
mpfr_exp_t e = MPFR_EXP(x);
mpfr_prec_t p = (prec - 1) / 32 + 1; /* ceil(prec/32) */
exp_enc = (e < -2147483648 || 2147483647 < e);
n += 4 + 4 * exp_enc + 4 * p;
}
t = s = (unsigned char *) malloc (n * sizeof (char));
t[0] = (mpfr_signbit (x) != 0) << 7;
if (MPFR_IS_NAN(x))
t[0] += 3 << 5;
else if (MPFR_IS_INF(x))
t[0] += 2 << 5;
else if (MPFR_IS_ZERO(x))
t[0] += 1 << 5;
t[0] += exp_enc << 4;
t[0] += prec_enc << 3;
if (prec_enc)
{
/* FIXME: the shift count may be too large for the type size. */
/* FIXME: check that the precision is not too large. */
t[0] += prec >> 56; /* reduction mod 8 is implicit */
t[1] = prec >> 48;
t[2] = prec >> 40;
t[3] = prec >> 32;
t += 4;
t[0] = 0;
}
t[0] += prec >> 24; /* reduction mod 256 is implicit */
t[1] = prec >> 16;
t[2] = prec >> 8;
t[3] = prec;
free (s);
return n;
}
|
