#
source:
trunk/third/gmp/mpf/ui_sub.c
@
18191

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

Line | |
---|---|

1 | /* mpf_ui_sub -- Subtract a float from an unsigned long int. |

2 | |

3 | Copyright 1993, 1994, 1995, 1996, 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 | void |

26 | mpf_ui_sub (mpf_ptr r, unsigned long int u, mpf_srcptr v) |

27 | { |

28 | mp_srcptr up, vp; |

29 | mp_ptr rp, tp; |

30 | mp_size_t usize, vsize, rsize; |

31 | mp_size_t prec; |

32 | mp_exp_t uexp; |

33 | mp_size_t ediff; |

34 | int negate; |

35 | mp_limb_t ulimb; |

36 | TMP_DECL (marker); |

37 | |

38 | vsize = v->_mp_size; |

39 | |

40 | /* Handle special cases that don't work in generic code below. */ |

41 | if (u == 0) |

42 | { |

43 | mpf_neg (r, v); |

44 | return; |

45 | } |

46 | if (vsize == 0) |

47 | { |

48 | mpf_set_ui (r, u); |

49 | return; |

50 | } |

51 | |

52 | /* If signs of U and V are different, perform addition. */ |

53 | if (vsize < 0) |

54 | { |

55 | __mpf_struct v_negated; |

56 | v_negated._mp_size = -vsize; |

57 | v_negated._mp_exp = v->_mp_exp; |

58 | v_negated._mp_d = v->_mp_d; |

59 | mpf_add_ui (r, &v_negated, u); |

60 | return; |

61 | } |

62 | |

63 | TMP_MARK (marker); |

64 | |

65 | /* Signs are now known to be the same. */ |

66 | |

67 | ulimb = u; |

68 | /* Make U be the operand with the largest exponent. */ |

69 | if (1 < v->_mp_exp) |

70 | { |

71 | negate = 1; |

72 | usize = ABS (vsize); |

73 | vsize = 1; |

74 | up = v->_mp_d; |

75 | vp = &ulimb; |

76 | rp = r->_mp_d; |

77 | prec = r->_mp_prec + 1; |

78 | uexp = v->_mp_exp; |

79 | ediff = uexp - 1; |

80 | } |

81 | else |

82 | { |

83 | negate = 0; |

84 | usize = 1; |

85 | vsize = ABS (vsize); |

86 | up = &ulimb; |

87 | vp = v->_mp_d; |

88 | rp = r->_mp_d; |

89 | prec = r->_mp_prec; |

90 | uexp = 1; |

91 | ediff = 1 - v->_mp_exp; |

92 | } |

93 | |

94 | /* Ignore leading limbs in U and V that are equal. Doing |

95 | this helps increase the precision of the result. */ |

96 | if (ediff == 0) |

97 | { |

98 | /* This loop normally exits immediately. Optimize for that. */ |

99 | for (;;) |

100 | { |

101 | usize--; |

102 | vsize--; |

103 | if (up[usize] != vp[vsize]) |

104 | break; |

105 | uexp--; |

106 | if (usize == 0) |

107 | goto Lu0; |

108 | if (vsize == 0) |

109 | goto Lv0; |

110 | } |

111 | usize++; |

112 | vsize++; |

113 | /* Note that either operand (but not both operands) might now have |

114 | leading zero limbs. It matters only that U is unnormalized if |

115 | vsize is now zero, and vice versa. And it is only in that case |

116 | that we have to adjust uexp. */ |

117 | if (vsize == 0) |

118 | Lv0: |

119 | while (usize != 0 && up[usize - 1] == 0) |

120 | usize--, uexp--; |

121 | if (usize == 0) |

122 | Lu0: |

123 | while (vsize != 0 && vp[vsize - 1] == 0) |

124 | vsize--, uexp--; |

125 | } |

126 | |

127 | /* If U extends beyond PREC, ignore the part that does. */ |

128 | if (usize > prec) |

129 | { |

130 | up += usize - prec; |

131 | usize = prec; |

132 | } |

133 | |

134 | /* If V extends beyond PREC, ignore the part that does. |

135 | Note that this may make vsize negative. */ |

136 | if (vsize + ediff > prec) |

137 | { |

138 | vp += vsize + ediff - prec; |

139 | vsize = prec - ediff; |

140 | } |

141 | |

142 | /* Allocate temp space for the result. Allocate |

143 | just vsize + ediff later??? */ |

144 | tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB); |

145 | |

146 | if (ediff >= prec) |

147 | { |

148 | /* V completely cancelled. */ |

149 | if (tp != up) |

150 | MPN_COPY (rp, up, usize); |

151 | rsize = usize; |

152 | } |

153 | else |

154 | { |

155 | /* Locate the least significant non-zero limb in (the needed |

156 | parts of) U and V, to simplify the code below. */ |

157 | for (;;) |

158 | { |

159 | if (vsize == 0) |

160 | { |

161 | MPN_COPY (rp, up, usize); |

162 | rsize = usize; |

163 | goto done; |

164 | } |

165 | if (vp[0] != 0) |

166 | break; |

167 | vp++, vsize--; |

168 | } |

169 | for (;;) |

170 | { |

171 | if (usize == 0) |

172 | { |

173 | MPN_COPY (rp, vp, vsize); |

174 | rsize = vsize; |

175 | negate ^= 1; |

176 | goto done; |

177 | } |

178 | if (up[0] != 0) |

179 | break; |

180 | up++, usize--; |

181 | } |

182 | |

183 | /* uuuu | uuuu | uuuu | uuuu | uuuu */ |

184 | /* vvvvvvv | vv | vvvvv | v | vv */ |

185 | |

186 | if (usize > ediff) |

187 | { |

188 | /* U and V partially overlaps. */ |

189 | if (ediff == 0) |

190 | { |

191 | /* Have to compare the leading limbs of u and v |

192 | to determine whether to compute u - v or v - u. */ |

193 | if (usize > vsize) |

194 | { |

195 | /* uuuu */ |

196 | /* vv */ |

197 | int cmp; |

198 | cmp = mpn_cmp (up + usize - vsize, vp, vsize); |

199 | if (cmp >= 0) |

200 | { |

201 | mp_size_t size; |

202 | size = usize - vsize; |

203 | MPN_COPY (tp, up, size); |

204 | mpn_sub_n (tp + size, up + size, vp, vsize); |

205 | rsize = usize; |

206 | } |

207 | else |

208 | { |

209 | /* vv */ /* Swap U and V. */ |

210 | /* uuuu */ |

211 | mp_size_t size, i; |

212 | size = usize - vsize; |

213 | tp[0] = -up[0] & GMP_NUMB_MASK; |

214 | for (i = 1; i < size; i++) |

215 | tp[i] = ~up[i] & GMP_NUMB_MASK; |

216 | mpn_sub_n (tp + size, vp, up + size, vsize); |

217 | mpn_sub_1 (tp + size, tp + size, vsize, (mp_limb_t) 1); |

218 | negate ^= 1; |

219 | rsize = usize; |

220 | } |

221 | } |

222 | else if (usize < vsize) |

223 | { |

224 | /* uuuu */ |

225 | /* vvvvvvv */ |

226 | int cmp; |

227 | cmp = mpn_cmp (up, vp + vsize - usize, usize); |

228 | if (cmp > 0) |

229 | { |

230 | mp_size_t size, i; |

231 | size = vsize - usize; |

232 | tp[0] = -vp[0] & GMP_NUMB_MASK; |

233 | for (i = 1; i < size; i++) |

234 | tp[i] = ~vp[i] & GMP_NUMB_MASK; |

235 | mpn_sub_n (tp + size, up, vp + size, usize); |

236 | mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1); |

237 | rsize = vsize; |

238 | } |

239 | else |

240 | { |

241 | /* vvvvvvv */ /* Swap U and V. */ |

242 | /* uuuu */ |

243 | /* This is the only place we can get 0.0. */ |

244 | mp_size_t size; |

245 | size = vsize - usize; |

246 | MPN_COPY (tp, vp, size); |

247 | mpn_sub_n (tp + size, vp + size, up, usize); |

248 | negate ^= 1; |

249 | rsize = vsize; |

250 | } |

251 | } |

252 | else |

253 | { |

254 | /* uuuu */ |

255 | /* vvvv */ |

256 | int cmp; |

257 | cmp = mpn_cmp (up, vp + vsize - usize, usize); |

258 | if (cmp > 0) |

259 | { |

260 | mpn_sub_n (tp, up, vp, usize); |

261 | rsize = usize; |

262 | } |

263 | else |

264 | { |

265 | mpn_sub_n (tp, vp, up, usize); |

266 | negate ^= 1; |

267 | rsize = usize; |

268 | /* can give zero */ |

269 | } |

270 | } |

271 | } |

272 | else |

273 | { |

274 | if (vsize + ediff <= usize) |

275 | { |

276 | /* uuuu */ |

277 | /* v */ |

278 | mp_size_t size; |

279 | size = usize - ediff - vsize; |

280 | MPN_COPY (tp, up, size); |

281 | mpn_sub (tp + size, up + size, usize - size, vp, vsize); |

282 | rsize = usize; |

283 | } |

284 | else |

285 | { |

286 | /* uuuu */ |

287 | /* vvvvv */ |

288 | mp_size_t size, i; |

289 | size = vsize + ediff - usize; |

290 | tp[0] = -vp[0] & GMP_NUMB_MASK; |

291 | for (i = 1; i < size; i++) |

292 | tp[i] = ~vp[i] & GMP_NUMB_MASK; |

293 | mpn_sub (tp + size, up, usize, vp + size, usize - ediff); |

294 | mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1); |

295 | rsize = vsize + ediff; |

296 | } |

297 | } |

298 | } |

299 | else |

300 | { |

301 | /* uuuu */ |

302 | /* vv */ |

303 | mp_size_t size, i; |

304 | size = vsize + ediff - usize; |

305 | tp[0] = -vp[0] & GMP_NUMB_MASK; |

306 | for (i = 1; i < vsize; i++) |

307 | tp[i] = ~vp[i] & GMP_NUMB_MASK; |

308 | for (i = vsize; i < size; i++) |

309 | tp[i] = GMP_NUMB_MAX; |

310 | mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1); |

311 | rsize = size + usize; |

312 | } |

313 | |

314 | /* Full normalize. Optimize later. */ |

315 | while (rsize != 0 && tp[rsize - 1] == 0) |

316 | { |

317 | rsize--; |

318 | uexp--; |

319 | } |

320 | MPN_COPY (rp, tp, rsize); |

321 | } |

322 | |

323 | done: |

324 | r->_mp_size = negate ? -rsize : rsize; |

325 | r->_mp_exp = uexp; |

326 | TMP_FREE (marker); |

327 | } |

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