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
/*
Copyright (C) 2014 Fredrik Johansson
This file is part of FLINT.
FLINT is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include "mpn_extras.h"
#include "arf.h"
void __gmpn_div_q(nn_ptr, nn_srcptr, slong, nn_srcptr, slong, nn_ptr);
static void
arf_div_special(arf_t z, const arf_t x, const arf_t y)
{
if ((arf_is_zero(x) && !arf_is_zero(y) && !arf_is_nan(y)) ||
(arf_is_inf(y) && !arf_is_special(x)))
{
arf_zero(z);
}
else if (arf_is_zero(y) || (arf_is_special(x) && arf_is_special(y)) ||
arf_is_nan(x) || arf_is_nan(y))
{
arf_nan(z);
}
else if (arf_sgn(x) == arf_sgn(y))
arf_pos_inf(z);
else
arf_neg_inf(z);
}
int
arf_div(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd)
{
slong xn, yn, zn, sn, tn, alloc;
nn_srcptr xptr, yptr;
nn_ptr tmp;
nn_ptr tptr, zptr;
int inexact;
slong fix, fix2;
ARF_MUL_TMP_DECL
if (arf_is_special(x) || arf_is_special(y))
{
arf_div_special(z, x, y);
return 0;
}
ARF_GET_MPN_READONLY(xptr, xn, x);
ARF_GET_MPN_READONLY(yptr, yn, y);
/* Division by a power of two */
if (yn == 1 && yptr[0] == LIMB_TOP)
{
fmpz_t t;
fmpz_init_set(t, ARF_EXPREF(y));
if (ARF_SGNBIT(y))
inexact = arf_neg_round(z, x, prec, rnd);
else
inexact = arf_set_round(z, x, prec, rnd);
_fmpz_sub2_fast(ARF_EXPREF(z), ARF_EXPREF(z), t, 1);
fmpz_clear(t);
return inexact;
}
sn = FLINT_MAX(prec - xn * FLINT_BITS + yn * FLINT_BITS, 0) + 32;
sn = (sn + FLINT_BITS - 1) / FLINT_BITS;
/* Need space for numerator (tn), quotient (zn), and one extra limb at the
top of tptr for the multiplication to check the remainder. */
tn = xn + sn;
zn = tn - yn + 1;
alloc = zn + (tn + 1);
/* need tn + 1 extra temporary limbs, which we store at the end of tptr */
alloc += tn + 1;
ARF_MUL_TMP_ALLOC(tmp, alloc)
zptr = tmp;
tptr = tmp + zn;
flint_mpn_zero(tptr, sn);
flint_mpn_copyi(tptr + sn, xptr, xn);
/* uses tn + 1 extra temporary limbs, tn limbs after tptr */
__gmpn_div_q(zptr, tptr, tn, yptr, yn, tptr + tn);
if (zptr[zn - 1] == 0)
{
zn--;
fix2 = 0;
}
else
{
fix2 = FLINT_BITS;
}
/* The remainder can only be zero if the last several bits of the
extended quotient are zero. */
if ((zptr[0] & ((LIMB_ONE << 24) - 1)) == 0)
{
/* The quotient is exact iff multiplying by y gives back the
input. Note: the multiplication may write sn + xn + 1 limbs, but
tptr[sn + xn] is guaranteed to be zero in that case since the
approximate quotient cannot be larger than the true quotient. */
if (zn >= yn)
flint_mpn_mul(tptr, zptr, zn, yptr, yn);
else
flint_mpn_mul(tptr, yptr, yn, zptr, zn);
/* The quotient is not exact. Perturbing the approximate quotient
and rounding gives the correct the result. */
if (!flint_mpn_zero_p(tptr, sn) ||
mpn_cmp(tptr + sn, xptr, xn) != 0)
{
zptr[0]++;
}
}
inexact = _arf_set_round_mpn(z, &fix, zptr, zn,
ARF_SGNBIT(x) ^ ARF_SGNBIT(y), prec, rnd);
_fmpz_sub2_fast(ARF_EXPREF(z), ARF_EXPREF(x), ARF_EXPREF(y), fix + fix2);
ARF_MUL_TMP_FREE(tmp, alloc)
return inexact;
}