Coverage Report

Created: 2024-11-21 06:47

/src/libgmp/mpn/hgcd2.c
Line
Count
Source (jump to first uncovered line)
1
/* hgcd2.c
2
3
   THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
4
   SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
5
   GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
6
7
Copyright 1996, 1998, 2000-2004, 2008, 2012, 2019 Free Software Foundation,
8
Inc.
9
10
This file is part of the GNU MP Library.
11
12
The GNU MP Library is free software; you can redistribute it and/or modify
13
it under the terms of either:
14
15
  * the GNU Lesser General Public License as published by the Free
16
    Software Foundation; either version 3 of the License, or (at your
17
    option) any later version.
18
19
or
20
21
  * the GNU General Public License as published by the Free Software
22
    Foundation; either version 2 of the License, or (at your option) any
23
    later version.
24
25
or both in parallel, as here.
26
27
The GNU MP Library is distributed in the hope that it will be useful, but
28
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
29
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
30
for more details.
31
32
You should have received copies of the GNU General Public License and the
33
GNU Lesser General Public License along with the GNU MP Library.  If not,
34
see https://www.gnu.org/licenses/.  */
35
36
#include "gmp-impl.h"
37
#include "longlong.h"
38
39
#include "mpn/generic/hgcd2-div.h"
40
41
#if GMP_NAIL_BITS != 0
42
#error Nails not implemented
43
#endif
44
45
/* Reduces a,b until |a-b| (almost) fits in one limb + 1 bit. Constructs
46
   matrix M. Returns 1 if we make progress, i.e. can perform at least
47
   one subtraction. Otherwise returns zero. */
48
49
/* FIXME: Possible optimizations:
50
51
   The div2 function starts with checking the most significant bit of
52
   the numerator. We can maintained normalized operands here, call
53
   hgcd with normalized operands only, which should make the code
54
   simpler and possibly faster.
55
56
   Experiment with table lookups on the most significant bits.
57
58
   This function is also a candidate for assembler implementation.
59
*/
60
int
61
mpn_hgcd2 (mp_limb_t ah, mp_limb_t al, mp_limb_t bh, mp_limb_t bl,
62
     struct hgcd_matrix1 *M)
63
87.9k
{
64
87.9k
  mp_limb_t u00, u01, u10, u11;
65
66
87.9k
  if (ah < 2 || bh < 2)
67
313
    return 0;
68
69
87.6k
  if (ah > bh || (ah == bh && al > bl))
70
42.6k
    {
71
42.6k
      sub_ddmmss (ah, al, ah, al, bh, bl);
72
42.6k
      if (ah < 2)
73
368
  return 0;
74
75
42.3k
      u00 = u01 = u11 = 1;
76
42.3k
      u10 = 0;
77
42.3k
    }
78
44.9k
  else
79
44.9k
    {
80
44.9k
      sub_ddmmss (bh, bl, bh, bl, ah, al);
81
44.9k
      if (bh < 2)
82
1.60k
  return 0;
83
84
43.3k
      u00 = u10 = u11 = 1;
85
43.3k
      u01 = 0;
86
43.3k
    }
87
88
85.6k
  if (ah < bh)
89
42.7k
    goto subtract_a;
90
91
42.8k
  for (;;)
92
811k
    {
93
811k
      ASSERT (ah >= bh);
94
811k
      if (ah == bh)
95
407
  goto done;
96
97
811k
      if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2)))
98
42.2k
  {
99
42.2k
    ah = (ah << (GMP_LIMB_BITS / 2) ) + (al >> (GMP_LIMB_BITS / 2));
100
42.2k
    bh = (bh << (GMP_LIMB_BITS / 2) ) + (bl >> (GMP_LIMB_BITS / 2));
101
102
42.2k
    break;
103
42.2k
  }
104
105
      /* Subtract a -= q b, and multiply M from the right by (1 q ; 0
106
   1), affecting the second column of M. */
107
769k
      ASSERT (ah > bh);
108
769k
      sub_ddmmss (ah, al, ah, al, bh, bl);
109
110
769k
      if (ah < 2)
111
44
  goto done;
112
113
769k
      if (ah <= bh)
114
306k
  {
115
    /* Use q = 1 */
116
306k
    u01 += u00;
117
306k
    u11 += u10;
118
306k
  }
119
462k
      else
120
462k
  {
121
462k
    mp_limb_t r[2];
122
462k
    mp_limb_t q = div2 (r, ah, al, bh, bl);
123
462k
    al = r[0]; ah = r[1];
124
462k
    if (ah < 2)
125
290
      {
126
        /* A is too small, but q is correct. */
127
290
        u01 += q * u00;
128
290
        u11 += q * u10;
129
290
        goto done;
130
290
      }
131
462k
    q++;
132
462k
    u01 += q * u00;
133
462k
    u11 += q * u10;
134
462k
  }
135
811k
    subtract_a:
136
811k
      ASSERT (bh >= ah);
137
811k
      if (ah == bh)
138
389
  goto done;
139
140
811k
      if (bh < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2)))
141
42.0k
  {
142
42.0k
    ah = (ah << (GMP_LIMB_BITS / 2) ) + (al >> (GMP_LIMB_BITS / 2));
143
42.0k
    bh = (bh << (GMP_LIMB_BITS / 2) ) + (bl >> (GMP_LIMB_BITS / 2));
144
145
42.0k
    goto subtract_a1;
146
42.0k
  }
147
148
      /* Subtract b -= q a, and multiply M from the right by (1 0 ; q
149
   1), affecting the first column of M. */
150
769k
      sub_ddmmss (bh, bl, bh, bl, ah, al);
151
152
769k
      if (bh < 2)
153
34
  goto done;
154
155
769k
      if (bh <= ah)
156
306k
  {
157
    /* Use q = 1 */
158
306k
    u00 += u01;
159
306k
    u10 += u11;
160
306k
  }
161
463k
      else
162
463k
  {
163
463k
    mp_limb_t r[2];
164
463k
    mp_limb_t q = div2 (r, bh, bl, ah, al);
165
463k
    bl = r[0]; bh = r[1];
166
463k
    if (bh < 2)
167
290
      {
168
        /* B is too small, but q is correct. */
169
290
        u00 += q * u01;
170
290
        u10 += q * u11;
171
290
        goto done;
172
290
      }
173
463k
    q++;
174
463k
    u00 += q * u01;
175
463k
    u10 += q * u11;
176
463k
  }
177
769k
    }
178
179
  /* NOTE: Since we discard the least significant half limb, we don't get a
180
     truly maximal M (corresponding to |a - b| < 2^{GMP_LIMB_BITS +1}). */
181
  /* Single precision loop */
182
42.2k
  for (;;)
183
712k
    {
184
712k
      ASSERT (ah >= bh);
185
186
712k
      ah -= bh;
187
712k
      if (ah < (CNST_LIMB (1) << (GMP_LIMB_BITS / 2 + 1)))
188
19.1k
  break;
189
190
693k
      if (ah <= bh)
191
289k
  {
192
    /* Use q = 1 */
193
289k
    u01 += u00;
194
289k
    u11 += u10;
195
289k
  }
196
404k
      else
197
404k
  {
198
404k
    mp_double_limb_t rq = div1 (ah, bh);
199
404k
    mp_limb_t q = rq.d1;
200
404k
    ah = rq.d0;
201
202
404k
    if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1)))
203
23.1k
      {
204
        /* A is too small, but q is correct. */
205
23.1k
        u01 += q * u00;
206
23.1k
        u11 += q * u10;
207
23.1k
        break;
208
23.1k
      }
209
381k
    q++;
210
381k
    u01 += q * u00;
211
381k
    u11 += q * u10;
212
381k
  }
213
712k
    subtract_a1:
214
712k
      ASSERT (bh >= ah);
215
216
712k
      bh -= ah;
217
712k
      if (bh < (CNST_LIMB (1) << (GMP_LIMB_BITS / 2 + 1)))
218
18.9k
  break;
219
220
693k
      if (bh <= ah)
221
289k
  {
222
    /* Use q = 1 */
223
289k
    u00 += u01;
224
289k
    u10 += u11;
225
289k
  }
226
404k
      else
227
404k
  {
228
404k
    mp_double_limb_t rq = div1 (bh, ah);
229
404k
    mp_limb_t q = rq.d1;
230
404k
    bh = rq.d0;
231
232
404k
    if (bh < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1)))
233
22.9k
      {
234
        /* B is too small, but q is correct. */
235
22.9k
        u00 += q * u01;
236
22.9k
        u10 += q * u11;
237
22.9k
        break;
238
22.9k
      }
239
381k
    q++;
240
381k
    u00 += q * u01;
241
381k
    u10 += q * u11;
242
381k
  }
243
693k
    }
244
245
85.6k
 done:
246
85.6k
  M->u[0][0] = u00; M->u[0][1] = u01;
247
85.6k
  M->u[1][0] = u10; M->u[1][1] = u11;
248
249
85.6k
  return 1;
250
42.2k
}
251
252
/* Sets (r;b) = (a;b) M, with M = (u00, u01; u10, u11). Vector must
253
 * have space for n + 1 limbs. Uses three buffers to avoid a copy*/
254
mp_size_t
255
mpn_hgcd_mul_matrix1_vector (const struct hgcd_matrix1 *M,
256
           mp_ptr rp, mp_srcptr ap, mp_ptr bp, mp_size_t n)
257
0
{
258
0
  mp_limb_t ah, bh;
259
260
  /* Compute (r,b) <-- (u00 a + u10 b, u01 a + u11 b) as
261
262
     r  = u00 * a
263
     r += u10 * b
264
     b *= u11
265
     b += u01 * a
266
  */
267
268
0
#if HAVE_NATIVE_mpn_addaddmul_1msb0
269
0
  ah = mpn_addaddmul_1msb0 (rp, ap, bp, n, M->u[0][0], M->u[1][0]);
270
0
  bh = mpn_addaddmul_1msb0 (bp, bp, ap, n, M->u[1][1], M->u[0][1]);
271
#else
272
  ah =     mpn_mul_1 (rp, ap, n, M->u[0][0]);
273
  ah += mpn_addmul_1 (rp, bp, n, M->u[1][0]);
274
275
  bh =     mpn_mul_1 (bp, bp, n, M->u[1][1]);
276
  bh += mpn_addmul_1 (bp, ap, n, M->u[0][1]);
277
#endif
278
0
  rp[n] = ah;
279
0
  bp[n] = bh;
280
281
0
  n += (ah | bh) > 0;
282
0
  return n;
283
0
}