1*86d7f5d3SJohn Marino /* mpf_div -- Divide two floats.
2*86d7f5d3SJohn Marino
3*86d7f5d3SJohn Marino Copyright 1993, 1994, 1996, 2000, 2001, 2002, 2004, 2005, 2010 Free Software
4*86d7f5d3SJohn Marino Foundation, Inc.
5*86d7f5d3SJohn Marino
6*86d7f5d3SJohn Marino This file is part of the GNU MP Library.
7*86d7f5d3SJohn Marino
8*86d7f5d3SJohn Marino The GNU MP Library is free software; you can redistribute it and/or modify
9*86d7f5d3SJohn Marino it under the terms of the GNU Lesser General Public License as published by
10*86d7f5d3SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
11*86d7f5d3SJohn Marino option) any later version.
12*86d7f5d3SJohn Marino
13*86d7f5d3SJohn Marino The GNU MP Library is distributed in the hope that it will be useful, but
14*86d7f5d3SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15*86d7f5d3SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16*86d7f5d3SJohn Marino License for more details.
17*86d7f5d3SJohn Marino
18*86d7f5d3SJohn Marino You should have received a copy of the GNU Lesser General Public License
19*86d7f5d3SJohn Marino along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
20*86d7f5d3SJohn Marino
21*86d7f5d3SJohn Marino #include "gmp.h"
22*86d7f5d3SJohn Marino #include "gmp-impl.h"
23*86d7f5d3SJohn Marino
24*86d7f5d3SJohn Marino
25*86d7f5d3SJohn Marino /* Not done:
26*86d7f5d3SJohn Marino
27*86d7f5d3SJohn Marino No attempt is made to identify an overlap u==v. The result will be
28*86d7f5d3SJohn Marino correct (1.0), but a full actual division is done whereas of course
29*86d7f5d3SJohn Marino x/x==1 needs no work. Such a call is not a sensible thing to make, and
30*86d7f5d3SJohn Marino it's left to an application to notice and optimize if it might arise
31*86d7f5d3SJohn Marino somehow through pointer aliasing or whatever.
32*86d7f5d3SJohn Marino
33*86d7f5d3SJohn Marino Enhancements:
34*86d7f5d3SJohn Marino
35*86d7f5d3SJohn Marino The high quotient limb is non-zero when high{up,vsize} >= {vp,vsize}. We
36*86d7f5d3SJohn Marino could make that comparison and use qsize==prec instead of qsize==prec+1,
37*86d7f5d3SJohn Marino to save one limb in the division.
38*86d7f5d3SJohn Marino
39*86d7f5d3SJohn Marino If r==u but the size is enough bigger than prec that there won't be an
40*86d7f5d3SJohn Marino overlap between quotient and dividend in mpn_tdiv_qr, then we can avoid
41*86d7f5d3SJohn Marino copying up,usize. This would only arise from a prec reduced with
42*86d7f5d3SJohn Marino mpf_set_prec_raw and will be pretty unusual, but might be worthwhile if
43*86d7f5d3SJohn Marino it could be worked into the copy_u decision cleanly. */
44*86d7f5d3SJohn Marino
45*86d7f5d3SJohn Marino void
mpf_div(mpf_ptr r,mpf_srcptr u,mpf_srcptr v)46*86d7f5d3SJohn Marino mpf_div (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
47*86d7f5d3SJohn Marino {
48*86d7f5d3SJohn Marino mp_srcptr up, vp;
49*86d7f5d3SJohn Marino mp_ptr rp, tp, new_vp;
50*86d7f5d3SJohn Marino mp_size_t usize, vsize, rsize, prospective_rsize, tsize, zeros;
51*86d7f5d3SJohn Marino mp_size_t sign_quotient, prec, high_zero, chop;
52*86d7f5d3SJohn Marino mp_exp_t rexp;
53*86d7f5d3SJohn Marino int copy_u;
54*86d7f5d3SJohn Marino TMP_DECL;
55*86d7f5d3SJohn Marino
56*86d7f5d3SJohn Marino usize = SIZ(u);
57*86d7f5d3SJohn Marino vsize = SIZ(v);
58*86d7f5d3SJohn Marino sign_quotient = usize ^ vsize;
59*86d7f5d3SJohn Marino usize = ABS (usize);
60*86d7f5d3SJohn Marino vsize = ABS (vsize);
61*86d7f5d3SJohn Marino prec = PREC(r);
62*86d7f5d3SJohn Marino
63*86d7f5d3SJohn Marino if (vsize == 0)
64*86d7f5d3SJohn Marino DIVIDE_BY_ZERO;
65*86d7f5d3SJohn Marino
66*86d7f5d3SJohn Marino if (usize == 0)
67*86d7f5d3SJohn Marino {
68*86d7f5d3SJohn Marino SIZ(r) = 0;
69*86d7f5d3SJohn Marino EXP(r) = 0;
70*86d7f5d3SJohn Marino return;
71*86d7f5d3SJohn Marino }
72*86d7f5d3SJohn Marino
73*86d7f5d3SJohn Marino TMP_MARK;
74*86d7f5d3SJohn Marino rexp = EXP(u) - EXP(v) + 1;
75*86d7f5d3SJohn Marino
76*86d7f5d3SJohn Marino rp = PTR(r);
77*86d7f5d3SJohn Marino up = PTR(u);
78*86d7f5d3SJohn Marino vp = PTR(v);
79*86d7f5d3SJohn Marino
80*86d7f5d3SJohn Marino prospective_rsize = usize - vsize + 1; /* quot from using given u,v sizes */
81*86d7f5d3SJohn Marino rsize = prec + 1; /* desired quot */
82*86d7f5d3SJohn Marino
83*86d7f5d3SJohn Marino zeros = rsize - prospective_rsize; /* padding u to give rsize */
84*86d7f5d3SJohn Marino copy_u = (zeros > 0 || rp == up); /* copy u if overlap or padding */
85*86d7f5d3SJohn Marino
86*86d7f5d3SJohn Marino chop = MAX (-zeros, 0); /* negative zeros means shorten u */
87*86d7f5d3SJohn Marino up += chop;
88*86d7f5d3SJohn Marino usize -= chop;
89*86d7f5d3SJohn Marino zeros += chop; /* now zeros >= 0 */
90*86d7f5d3SJohn Marino
91*86d7f5d3SJohn Marino tsize = usize + zeros; /* size for possible copy of u */
92*86d7f5d3SJohn Marino
93*86d7f5d3SJohn Marino /* copy and possibly extend u if necessary */
94*86d7f5d3SJohn Marino if (copy_u)
95*86d7f5d3SJohn Marino {
96*86d7f5d3SJohn Marino tp = TMP_ALLOC_LIMBS (tsize + 1); /* +1 for mpn_div_q's scratch needs */
97*86d7f5d3SJohn Marino MPN_ZERO (tp, zeros);
98*86d7f5d3SJohn Marino MPN_COPY (tp+zeros, up, usize);
99*86d7f5d3SJohn Marino up = tp;
100*86d7f5d3SJohn Marino usize = tsize;
101*86d7f5d3SJohn Marino }
102*86d7f5d3SJohn Marino else
103*86d7f5d3SJohn Marino {
104*86d7f5d3SJohn Marino tp = TMP_ALLOC_LIMBS (usize + 1);
105*86d7f5d3SJohn Marino }
106*86d7f5d3SJohn Marino
107*86d7f5d3SJohn Marino /* ensure divisor doesn't overlap quotient */
108*86d7f5d3SJohn Marino if (rp == vp)
109*86d7f5d3SJohn Marino {
110*86d7f5d3SJohn Marino new_vp = TMP_ALLOC_LIMBS (vsize);
111*86d7f5d3SJohn Marino MPN_COPY (new_vp, vp, vsize);
112*86d7f5d3SJohn Marino vp = new_vp;
113*86d7f5d3SJohn Marino }
114*86d7f5d3SJohn Marino
115*86d7f5d3SJohn Marino ASSERT (usize-vsize+1 == rsize);
116*86d7f5d3SJohn Marino mpn_div_q (rp, up, usize, vp, vsize, tp);
117*86d7f5d3SJohn Marino
118*86d7f5d3SJohn Marino /* strip possible zero high limb */
119*86d7f5d3SJohn Marino high_zero = (rp[rsize-1] == 0);
120*86d7f5d3SJohn Marino rsize -= high_zero;
121*86d7f5d3SJohn Marino rexp -= high_zero;
122*86d7f5d3SJohn Marino
123*86d7f5d3SJohn Marino SIZ(r) = sign_quotient >= 0 ? rsize : -rsize;
124*86d7f5d3SJohn Marino EXP(r) = rexp;
125*86d7f5d3SJohn Marino TMP_FREE;
126*86d7f5d3SJohn Marino }
127