#
source:
trunk/third/gmp/extract-dbl.c
@
18191

Revision 18191, 6.4 KB checked in by ghudson, 22 years ago (diff) |
---|

Line | |
---|---|

1 | /* __gmp_extract_double -- convert from double to array of mp_limb_t. |

2 | |

3 | Copyright 1996, 1999, 2000, 2001, 2002 Free Software Foundation, Inc. |

4 | |

5 | This file is part of the GNU MP Library. |

6 | |

7 | The GNU MP Library is free software; you can redistribute it and/or modify |

8 | it under the terms of the GNU Lesser General Public License as published by |

9 | the Free Software Foundation; either version 2.1 of the License, or (at your |

10 | option) any later version. |

11 | |

12 | The GNU MP Library is distributed in the hope that it will be useful, but |

13 | WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |

14 | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |

15 | License for more details. |

16 | |

17 | You should have received a copy of the GNU Lesser General Public License |

18 | along with the GNU MP Library; see the file COPYING.LIB. If not, write to |

19 | the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |

20 | MA 02111-1307, USA. */ |

21 | |

22 | #include "gmp.h" |

23 | #include "gmp-impl.h" |

24 | |

25 | #ifdef XDEBUG |

26 | #undef _GMP_IEEE_FLOATS |

27 | #endif |

28 | |

29 | #ifndef _GMP_IEEE_FLOATS |

30 | #define _GMP_IEEE_FLOATS 0 |

31 | #endif |

32 | |

33 | #define BITS_IN_MANTISSA 53 |

34 | |

35 | /* Extract a non-negative double in d. */ |

36 | |

37 | int |

38 | __gmp_extract_double (mp_ptr rp, double d) |

39 | { |

40 | long exp; |

41 | unsigned sc; |

42 | #ifdef _LONG_LONG_LIMB |

43 | #define BITS_PER_PART 64 /* somewhat bogus */ |

44 | unsigned long long int manl; |

45 | #else |

46 | #define BITS_PER_PART GMP_LIMB_BITS |

47 | unsigned long int manh, manl; |

48 | #endif |

49 | |

50 | /* BUGS |

51 | |

52 | 1. Should handle Inf and NaN in IEEE specific code. |

53 | 2. Handle Inf and NaN also in default code, to avoid hangs. |

54 | 3. Generalize to handle all GMP_LIMB_BITS >= 32. |

55 | 4. This lits is incomplete and misspelled. |

56 | */ |

57 | |

58 | ASSERT (d >= 0.0); |

59 | |

60 | if (d == 0.0) |

61 | { |

62 | MPN_ZERO (rp, LIMBS_PER_DOUBLE); |

63 | return 0; |

64 | } |

65 | |

66 | #if _GMP_IEEE_FLOATS |

67 | { |

68 | #if defined (__alpha) && __GNUC__ == 2 && __GNUC_MINOR__ == 8 |

69 | /* Work around alpha-specific bug in GCC 2.8.x. */ |

70 | volatile |

71 | #endif |

72 | union ieee_double_extract x; |

73 | x.d = d; |

74 | exp = x.s.exp; |

75 | #if BITS_PER_PART == 64 /* generalize this to BITS_PER_PART > BITS_IN_MANTISSA */ |

76 | manl = (((mp_limb_t) 1 << 63) |

77 | | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11)); |

78 | if (exp == 0) |

79 | { |

80 | /* Denormalized number. Don't try to be clever about this, |

81 | since it is not an important case to make fast. */ |

82 | exp = 1; |

83 | do |

84 | { |

85 | manl = manl << 1; |

86 | exp--; |

87 | } |

88 | while ((mp_limb_signed_t) manl >= 0); |

89 | } |

90 | #endif |

91 | #if BITS_PER_PART == 32 |

92 | manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21); |

93 | manl = x.s.manl << 11; |

94 | if (exp == 0) |

95 | { |

96 | /* Denormalized number. Don't try to be clever about this, |

97 | since it is not an important case to make fast. */ |

98 | exp = 1; |

99 | do |

100 | { |

101 | manh = (manh << 1) | (manl >> 31); |

102 | manl = manl << 1; |

103 | exp--; |

104 | } |

105 | while ((mp_limb_signed_t) manh >= 0); |

106 | } |

107 | #endif |

108 | #if BITS_PER_PART != 32 && BITS_PER_PART != 64 |

109 | You need to generalize the code above to handle this. |

110 | #endif |

111 | exp -= 1022; /* Remove IEEE bias. */ |

112 | } |

113 | #else |

114 | { |

115 | /* Unknown (or known to be non-IEEE) double format. */ |

116 | exp = 0; |

117 | if (d >= 1.0) |

118 | { |

119 | if (d * 0.5 == d) |

120 | abort (); |

121 | |

122 | while (d >= 32768.0) |

123 | { |

124 | d *= (1.0 / 65536.0); |

125 | exp += 16; |

126 | } |

127 | while (d >= 1.0) |

128 | { |

129 | d *= 0.5; |

130 | exp += 1; |

131 | } |

132 | } |

133 | else if (d < 0.5) |

134 | { |

135 | while (d < (1.0 / 65536.0)) |

136 | { |

137 | d *= 65536.0; |

138 | exp -= 16; |

139 | } |

140 | while (d < 0.5) |

141 | { |

142 | d *= 2.0; |

143 | exp -= 1; |

144 | } |

145 | } |

146 | |

147 | d *= (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2))); |

148 | #if BITS_PER_PART == 64 |

149 | manl = d; |

150 | #else |

151 | manh = d; |

152 | manl = (d - manh) * (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2))); |

153 | #endif |

154 | } |

155 | #endif /* IEEE */ |

156 | |

157 | /* Up until here, we have ignored the actual limb size. Remains |

158 | to split manh,,manl into an array of LIMBS_PER_DOUBLE limbs. |

159 | */ |

160 | |

161 | sc = (unsigned) (exp + 64 * GMP_NUMB_BITS) % GMP_NUMB_BITS; |

162 | |

163 | /* We add something here to get rounding right. */ |

164 | exp = (exp + 64 * GMP_NUMB_BITS) / GMP_NUMB_BITS - 64 * GMP_NUMB_BITS / GMP_NUMB_BITS + 1; |

165 | |

166 | #if LIMBS_PER_DOUBLE == 2 |

167 | #if GMP_NAIL_BITS == 0 |

168 | if (sc != 0) |

169 | { |

170 | rp[1] = manl >> (GMP_LIMB_BITS - sc); |

171 | rp[0] = manl << sc; |

172 | } |

173 | else |

174 | { |

175 | rp[1] = manl; |

176 | rp[0] = 0; |

177 | exp--; |

178 | } |

179 | #else |

180 | if (sc > GMP_NAIL_BITS) |

181 | { |

182 | rp[1] = manl >> (GMP_LIMB_BITS - sc); |

183 | rp[0] = (manl << (sc - GMP_NAIL_BITS)) & GMP_NUMB_MASK; |

184 | } |

185 | else |

186 | { |

187 | if (sc == 0) |

188 | { |

189 | rp[1] = manl >> GMP_NAIL_BITS; |

190 | rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS) & GMP_NUMB_MASK; |

191 | exp--; |

192 | } |

193 | else |

194 | { |

195 | rp[1] = manl >> (GMP_LIMB_BITS - sc); |

196 | rp[0] = (manl >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK; |

197 | } |

198 | } |

199 | #endif |

200 | #endif |

201 | |

202 | #if LIMBS_PER_DOUBLE == 3 |

203 | #if GMP_NAIL_BITS == 0 |

204 | if (sc != 0) |

205 | { |

206 | rp[2] = manh >> (GMP_LIMB_BITS - sc); |

207 | rp[1] = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc)); |

208 | rp[0] = manl << sc; |

209 | } |

210 | else |

211 | { |

212 | rp[2] = manh; |

213 | rp[1] = manl; |

214 | rp[0] = 0; |

215 | exp--; |

216 | } |

217 | #else |

218 | if (sc > GMP_NAIL_BITS) |

219 | { |

220 | rp[2] = (manh >> (GMP_LIMB_BITS - sc)); |

221 | rp[1] = ((manh << (sc - GMP_NAIL_BITS)) | |

222 | (manl >> (GMP_LIMB_BITS - sc + GMP_NAIL_BITS))) & GMP_NUMB_MASK; |

223 | if (sc >= 2 * GMP_NAIL_BITS) |

224 | rp[0] = (manl << sc - 2 * GMP_NAIL_BITS) & GMP_NUMB_MASK; |

225 | else |

226 | rp[0] = manl >> (2 * GMP_NAIL_BITS - sc) & GMP_NUMB_MASK; |

227 | } |

228 | else |

229 | { |

230 | if (sc == 0) |

231 | { |

232 | rp[2] = manh >> GMP_NAIL_BITS; |

233 | rp[1] = ((manh << GMP_NUMB_BITS - GMP_NAIL_BITS) | (manl >> 2 * GMP_NAIL_BITS)) & GMP_NUMB_MASK; |

234 | rp[0] = 0; |

235 | exp--; |

236 | } |

237 | else |

238 | { |

239 | rp[2] = (manh >> (GMP_LIMB_BITS - sc)); |

240 | rp[1] = (manh >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK; |

241 | rp[0] = ((manh << (GMP_NUMB_BITS - GMP_NAIL_BITS + sc)) |

242 | | (manl >> (GMP_LIMB_BITS - (GMP_NUMB_BITS - GMP_NAIL_BITS + sc)))) & GMP_NUMB_MASK; |

243 | } |

244 | } |

245 | #endif |

246 | #endif |

247 | |

248 | #if LIMBS_PER_DOUBLE > 3 |

249 | /* Insert code for splitting manh,,manl into LIMBS_PER_DOUBLE |

250 | mp_limb_t's at rp. */ |

251 | if (sc != 0) |

252 | { |

253 | /* This is not perfect, and would fail for GMP_LIMB_BITS == 16. |

254 | The ASSERT_ALWAYS should catch the problematic cases. */ |

255 | ASSERT_ALWAYS ((manl << sc) == 0); |

256 | manl = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc)); |

257 | manh = manh >> (GMP_LIMB_BITS - sc); |

258 | } |

259 | { |

260 | int i; |

261 | for (i = LIMBS_PER_DOUBLE - 1; i >= 0; i--) |

262 | { |

263 | rp[i] = manh >> (BITS_PER_LONGINT - GMP_LIMB_BITS); |

264 | manh = ((manh << GMP_LIMB_BITS) |

265 | | (manl >> (BITS_PER_LONGINT - GMP_LIMB_BITS))); |

266 | manl = manl << GMP_LIMB_BITS; |

267 | } |

268 | } |

269 | #endif |

270 | |

271 | return exp; |

272 | } |

**Note:**See TracBrowser for help on using the repository browser.