1/* Implementations of operations between mpfr and mpz/mpq data 2 3Copyright 2001, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#define MPFR_NEED_LONGLONG_H 24#include "mpfr-impl.h" 25 26/* Init and set a mpfr_t with enough precision to store a mpz */ 27static void 28init_set_z (mpfr_ptr t, mpz_srcptr z) 29{ 30 mpfr_prec_t p; 31 int i; 32 33 if (mpz_size (z) <= 1) 34 p = GMP_NUMB_BITS; 35 else 36 MPFR_MPZ_SIZEINBASE2 (p, z); 37 mpfr_init2 (t, p); 38 i = mpfr_set_z (t, z, MPFR_RNDN); 39 MPFR_ASSERTD (i == 0); 40} 41 42/* Init, set a mpfr_t with enough precision to store a mpz_t without round, 43 call the function, and clear the allocated mpfr_t */ 44static int 45foo (mpfr_ptr x, mpfr_srcptr y, mpz_srcptr z, mpfr_rnd_t r, 46 int (*f)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t)) 47{ 48 mpfr_t t; 49 int i; 50 init_set_z (t, z); 51 i = (*f) (x, y, t, r); 52 mpfr_clear (t); 53 return i; 54} 55 56int 57mpfr_mul_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r) 58{ 59 return foo (y, x, z, r, mpfr_mul); 60} 61 62int 63mpfr_div_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r) 64{ 65 return foo (y, x, z, r, mpfr_div); 66} 67 68int 69mpfr_add_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r) 70{ 71 /* Mpz 0 is unsigned */ 72 if (MPFR_UNLIKELY (mpz_sgn (z) == 0)) 73 return mpfr_set (y, x, r); 74 else 75 return foo (y, x, z, r, mpfr_add); 76} 77 78int 79mpfr_sub_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r) 80{ 81 /* Mpz 0 is unsigned */ 82 if (MPFR_UNLIKELY (mpz_sgn (z) == 0)) 83 return mpfr_set (y, x, r); 84 else 85 return foo (y, x, z, r, mpfr_sub); 86} 87 88int 89mpfr_cmp_z (mpfr_srcptr x, mpz_srcptr z) 90{ 91 mpfr_t t; 92 int res; 93 init_set_z (t, z); 94 res = mpfr_cmp (x, t); 95 mpfr_clear (t); 96 return res; 97} 98 99int 100mpfr_mul_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode) 101{ 102 mpfr_t tmp; 103 int res; 104 mpfr_prec_t p; 105 106 if (MPFR_UNLIKELY (mpq_sgn (z) == 0)) 107 return mpfr_mul_ui (y, x, 0, rnd_mode); 108 else 109 { 110 MPFR_MPZ_SIZEINBASE2 (p, mpq_numref (z)); 111 mpfr_init2 (tmp, MPFR_PREC (x) + p); 112 res = mpfr_mul_z (tmp, x, mpq_numref(z), MPFR_RNDN ); 113 MPFR_ASSERTD (res == 0); 114 res = mpfr_div_z (y, tmp, mpq_denref(z), rnd_mode); 115 mpfr_clear (tmp); 116 return res; 117 } 118} 119 120int 121mpfr_div_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode) 122{ 123 mpfr_t tmp; 124 int res; 125 mpfr_prec_t p; 126 127 if (MPFR_UNLIKELY (mpq_sgn (z) == 0)) 128 return mpfr_div_ui (y, x, 0, rnd_mode); 129 else if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (z)) == 0)) 130 p = 0; 131 else 132 MPFR_MPZ_SIZEINBASE2 (p, mpq_denref (z)); 133 mpfr_init2 (tmp, MPFR_PREC(x) + p); 134 res = mpfr_mul_z (tmp, x, mpq_denref(z), MPFR_RNDN ); 135 MPFR_ASSERTD( res == 0 ); 136 res = mpfr_div_z (y, tmp, mpq_numref(z), rnd_mode); 137 mpfr_clear (tmp); 138 return res; 139} 140 141int 142mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode) 143{ 144 mpfr_t t,q; 145 mpfr_prec_t p; 146 mpfr_exp_t err; 147 int res; 148 MPFR_ZIV_DECL (loop); 149 150 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 151 { 152 if (MPFR_IS_NAN (x)) 153 { 154 MPFR_SET_NAN (y); 155 MPFR_RET_NAN; 156 } 157 else if (MPFR_IS_INF (x)) 158 { 159 MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0); 160 MPFR_SET_INF (y); 161 MPFR_SET_SAME_SIGN (y, x); 162 MPFR_RET (0); 163 } 164 else 165 { 166 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 167 if (MPFR_UNLIKELY (mpq_sgn (z) == 0)) 168 return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */ 169 else 170 return mpfr_set_q (y, z, rnd_mode); 171 } 172 } 173 174 p = MPFR_PREC (y) + 10; 175 mpfr_init2 (t, p); 176 mpfr_init2 (q, p); 177 178 MPFR_ZIV_INIT (loop, p); 179 for (;;) 180 { 181 res = mpfr_set_q (q, z, MPFR_RNDN); /* Error <= 1/2 ulp(q) */ 182 /* If z if @INF@ (1/0), res = 0, so it quits immediately */ 183 if (MPFR_UNLIKELY (res == 0)) 184 /* Result is exact so we can add it directly! */ 185 { 186 res = mpfr_add (y, x, q, rnd_mode); 187 break; 188 } 189 mpfr_add (t, x, q, MPFR_RNDN); /* Error <= 1/2 ulp(t) */ 190 /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t)) 191 If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t))) 192 <= 2^(EXP(q)-EXP(t)) 193 If EXP(q)-EXP(t)<0, <= 2^0 */ 194 /* We can get 0, but we can't round since q is inexact */ 195 if (MPFR_LIKELY (!MPFR_IS_ZERO (t))) 196 { 197 err = (mpfr_exp_t) p - 1 - MAX (MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0); 198 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, MPFR_PREC (y), rnd_mode))) 199 { 200 res = mpfr_set (y, t, rnd_mode); 201 break; 202 } 203 } 204 MPFR_ZIV_NEXT (loop, p); 205 mpfr_set_prec (t, p); 206 mpfr_set_prec (q, p); 207 } 208 MPFR_ZIV_FREE (loop); 209 mpfr_clear (t); 210 mpfr_clear (q); 211 return res; 212} 213 214int 215mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mpfr_rnd_t rnd_mode) 216{ 217 mpfr_t t,q; 218 mpfr_prec_t p; 219 int res; 220 mpfr_exp_t err; 221 MPFR_ZIV_DECL (loop); 222 223 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 224 { 225 if (MPFR_IS_NAN (x)) 226 { 227 MPFR_SET_NAN (y); 228 MPFR_RET_NAN; 229 } 230 else if (MPFR_IS_INF (x)) 231 { 232 MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0); 233 MPFR_SET_INF (y); 234 MPFR_SET_SAME_SIGN (y, x); 235 MPFR_RET (0); 236 } 237 else 238 { 239 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 240 241 if (MPFR_UNLIKELY (mpq_sgn (z) == 0)) 242 return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */ 243 else 244 { 245 res = mpfr_set_q (y, z, MPFR_INVERT_RND (rnd_mode)); 246 MPFR_CHANGE_SIGN (y); 247 return -res; 248 } 249 } 250 } 251 252 p = MPFR_PREC (y) + 10; 253 mpfr_init2 (t, p); 254 mpfr_init2 (q, p); 255 256 MPFR_ZIV_INIT (loop, p); 257 for(;;) 258 { 259 res = mpfr_set_q(q, z, MPFR_RNDN); /* Error <= 1/2 ulp(q) */ 260 /* If z if @INF@ (1/0), res = 0, so it quits immediately */ 261 if (MPFR_UNLIKELY (res == 0)) 262 /* Result is exact so we can add it directly!*/ 263 { 264 res = mpfr_sub (y, x, q, rnd_mode); 265 break; 266 } 267 mpfr_sub (t, x, q, MPFR_RNDN); /* Error <= 1/2 ulp(t) */ 268 /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t)) 269 If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t))) 270 <= 2^(EXP(q)-EXP(t)) 271 If EXP(q)-EXP(t)<0, <= 2^0 */ 272 /* We can get 0, but we can't round since q is inexact */ 273 if (MPFR_LIKELY (!MPFR_IS_ZERO (t))) 274 { 275 err = (mpfr_exp_t) p - 1 - MAX (MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0); 276 res = MPFR_CAN_ROUND (t, err, MPFR_PREC (y), rnd_mode); 277 if (MPFR_LIKELY (res != 0)) /* We can round! */ 278 { 279 res = mpfr_set (y, t, rnd_mode); 280 break; 281 } 282 } 283 MPFR_ZIV_NEXT (loop, p); 284 mpfr_set_prec (t, p); 285 mpfr_set_prec (q, p); 286 } 287 MPFR_ZIV_FREE (loop); 288 mpfr_clear (t); 289 mpfr_clear (q); 290 return res; 291} 292 293int 294mpfr_cmp_q (mpfr_srcptr x, mpq_srcptr z) 295{ 296 mpfr_t t; 297 int res; 298 mpfr_prec_t p; 299 /* x < a/b ? <=> x*b < a */ 300 MPFR_ASSERTD (mpz_sgn (mpq_denref (z)) != 0); 301 MPFR_MPZ_SIZEINBASE2 (p, mpq_denref (z)); 302 mpfr_init2 (t, MPFR_PREC(x) + p); 303 res = mpfr_mul_z (t, x, mpq_denref (z), MPFR_RNDN ); 304 MPFR_ASSERTD (res == 0); 305 res = mpfr_cmp_z (t, mpq_numref (z) ); 306 mpfr_clear (t); 307 return res; 308} 309 310int 311mpfr_cmp_f (mpfr_srcptr x, mpf_srcptr z) 312{ 313 mpfr_t t; 314 int res; 315 316 mpfr_init2 (t, MPFR_PREC_MIN + ABS(SIZ(z)) * GMP_NUMB_BITS ); 317 res = mpfr_set_f (t, z, MPFR_RNDN); 318 MPFR_ASSERTD (res == 0); 319 res = mpfr_cmp (x, t); 320 mpfr_clear (t); 321 return res; 322} 323