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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
|
Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation.
Contributed by the Spaces project, INRIA Lorraine.
This file is part of the MPFR Library.
The MPFR Library is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License (either version 2.1
of the License, or, at your option, any later version) and the GNU General
Public License as published by the Free Software Foundation (most of MPFR is
under the former, some under the latter).
The 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 MPFR Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
02111-1307, USA.
##############################################################################
Documentation:
- add a description of the algorithms used + proof of correctness
Installation:
- from Kevin Ryde <user42@zip.com.au>:
Determine the exp2/exp3 thresholds using tune/tuneup.c.
A start has been made on this, but there's a noticable step-effect
in the times making them cross back and forward between which is
faster. Hopefully this will go away with improvements to the exp
code.
Changes in existing functions:
- use (LDBL_MAX/2)+(LDBL_MAX/4*LDBL_EPSILON) instead of loop in tset_ld
to get largest power of 2 that fits in a long double
[and avoid overflow exceptions in tests, for example adding
#if MPFR_HAVE_TESTS_x86
x86_fldcw (x86_fstcw () & ~8);
#endif
in mpfr_test_init to enable overflow exceptions]
- allow the ISO C "P" exponent separator for base 16 in input functions
- invert 2nd and 3rd arguments of mpfr_round_prec (to be coherent with
mpfr_set_prec)
- merge mpfr_inp_str and mpfr_set_str (cf glibc sscanf/fscanf)
[mpfr_set_str -> with final '\0',
mpfr_strtofr -> without final '\0', returns the number of characters
read (or name mpfr_strtod suggested by Kevin, see below)]
Suggestion from Kevin:
- Function: char * mpf_strto (mpz_t ROP, const char *STR)
Read a floating point number from a string. If successful the
result is stored in ROP and the return value points to the
character after those parsed. If STR doesn't start with a valid
number then the return value is `NULL' and the value in ROP is
undefined.
Parsing follows the standard C `strtod' function (*note Parsing of
Integers: (libc)Parsing of Integers.). This means optional
leading whitespace, an optional `+' or `-', mantissa digits, and an
optional exponent consisting of an `e' or `E', an optional sign,
and digits. A hex mantissa can be given with a leading `0x' or
`0X', in which case `p' or `P' introduces the exponent (still in
decimal). In addition `inf' or `infinity' with an optional sign,
or `nan' or `nan(..chars..)', all non case significant, can be
given.
There must be at least one digit in the mantissa for the number to
be valid. If an exponent has no digits it's ignored and parsing
stops after the mantissa. If an `0x' or `0X' is not followed by
hexadecimal digits, parsing stops after the `0'.
Note that in the hex format the exponent represents a power of 2,
ie. the result is MANTISSA*2^EXPONENT. This is as per the `%a'
format in `printf' (*note Formatted Output Strings::).
`mpf_t' does not currently support NaNs or infinities. The value
stored to ROP for these is undefined.
The decimal point character, or string, expected is taken from the
current locale on systems providing `localeconv'. [end of suggestion]
- in mpfr_set_str, make string comparisons case insensitive and use @Inf@
and @NaN@ instead of the possibly ambiguous Inf and NaN in bases > 10.
Modify the other functions to make them consistent with this choice.
Possibly accept other strings, like those accepted by strtod.
New functions to implement:
- modf (to extract integer and fractional parts), suggested by
Dmitry Antipov <dmitry.antipov@mail.ru> Thu, 13 Jun 2002
- mpz_set_fr (to set a mpz from a mpfr, with a rounding mode)
-> cf code in leibniz:mpfr/mpz_set_fr.c
- those from LIA: missing secant, cosecant, cotangent (trigo/hyperbolic)
- atan2
- mpfr_set_float (MPF -> float) or mpfr_strtof (idem for mpfr_get_*)
- implement accurate summation algorithms from Demmel
(http://www.cs.berkeley.edu/~demmel/AccurateSummation.ps)
Rounding:
- mpfr_pow isn't completely specified (concerning signed zeros).
Efficiency:
- implement range reduction in sin/cos/tan for large arguments
(currently too slow for 2^1024)
- mpfr_asin/acos are too slow for small values (2^(-1021) for example)
- idem for mpfr_atanh (2^(-1021) for example)
- improve generic.c to work for number of terms <> 2^k
- rewrite mpfr_greater_p...
- combine tests for singular values [from Torbjorn Granlund]:
For example, mpfr_add could have this structure:
if (MPFR_SINGULAR (a))
...
else if (MPFR_SINGULAR (b)
...
else
plain operands
- mpf_t uses a scheme where the number of limbs actually present can
be less than the selected precision, thereby allowing low precision
values (for instance small integers) to be stored and manipulated in
an mpf_t efficiently.
Perhaps mpfr should get something similar, especially if looking to
replace mpf with mpfr, though it'd be a major change. Alternately
perhaps those mpfr routines like mpfr_mul where optimizations are
possible through stripping low zero bits or limbs could check for
that (this would be less efficient but easier).
Miscellaneous:
- rename mpf2mpfr.h to gmp-mpf2mpfr.h?
- from Kevin Ryde <user42@zip.com.au>:
Also for pi.c, a pre-calculated compiled-in pi to a few thousand
digits would be good value I think. After all, say 10000 bits using
1250 bytes would still be small compared to the code size!
Store pi in round to zero mode (to recover other modes).
- problem when reading a float followed by a character, for example 1.5*x
[from Fabrice.Rouillier@loria.fr, Mon, 04 Dec 2000]
- use AC_CHECK_FUNCS for __setfpucw rather than grepping a header file.
- rename (and rewrite) mpfr_isinteger to mpfr_integer_p.
(Kevin Ryde, 05 Mar 2002)
- use GMP rand/srand functions instead for SEED_RAND (time(NULL)). Kevin:
GMP tests/misc.c has a scheme to print the seed chosen and allow for a
re-run with a given seed. If you moved towards _gmp_rand and friends
at some stage you could probably use that scheme almost unchanged.
An advantage of _gmp_rand of course is that it's the same on all
systems, apart from limb size, so you're not at the mercy of random or
lrand48, etc.
- randomized tests: do like GMP (see misc.c), but using MPFR_CHECK_RANDOMIZE
- add tests for nextabove/nextbelow/nexttoward.
- documentation and tests for:
= < > unordered
mpfr_greater_p 0 0 1 0
mpfr_greaterequal_p 1 0 1 0
mpfr_less_p 0 1 0 0
mpfr_lessequal_p 1 1 0 0
mpfr_lessgreater_p 0 1 1 0
mpfr_equal_p 1 0 0 0
mpfr_unordered_p 0 0 0 1
Reentrancy / Thread-Safety:
- Temporary changes to emin/emax are not safe (all uses of
mpfr_save_emin_emax, eg. mpfr_set_q).
- Global variables for caching in mpfr_const_log2 and mpfr_const_pi
are not safe.
Portability:
- MINGW has neither random() nor rand48(), making the setups in
mpfr-test.h fail. It does have rand(), returning just 15 bits.
Switching to gmp_randstate_t etc would be better.
- _mpfr_ceil_log2 is IEEE specific. Many callers only seem to want
this operation on an integer, which could be done with
count_leading_zeros from longlong.h.
Possible future MPF / MPFR integration:
- mpf routines can become "extern inline"s calling mpfr equivalents,
probably just with GMP_RNDZ hard coded, since that's what mpf has
always done.
- Want to preserve the mpf_t size, for binary compatibility. Breaking
compatibility would cause lots of pain and potential subtle breakage
for users. If the fields in mpf_t are not enough then extra space
under _mp_d can be used.
- mpf_sgn has been a macro directly accessing the _mp_size field, so a
compatible representation would be required. At worst that field
could be maintained for mpf_sgn, but not otherwise used internally.
mpf_sgn should probably throw an exception if called with NaN, since
there's no useful value it can return, so it might want to become a
function. Inlined copies in existing binaries would hopefully never
see a NaN, if they only do old-style mpf things.
- mpfr routines replacing mpf routines must be reentrant and thread
safe, since of course that's what has been documented for mpf.
|
