Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpn/hgcd2.c
Line
Count
Source
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
50.7k
{
64
50.7k
  mp_limb_t u00, u01, u10, u11;
65
66
50.7k
  if (ah < 2 || bh < 2)
67
100
    return 0;
68
69
50.6k
  if (ah > bh || (ah == bh && al > bl))
70
24.4k
    {
71
24.4k
      sub_ddmmss (ah, al, ah, al, bh, bl);
72
24.4k
      if (ah < 2)
73
196
  return 0;
74
75
24.2k
      u00 = u01 = u11 = 1;
76
24.2k
      u10 = 0;
77
24.2k
    }
78
26.2k
  else
79
26.2k
    {
80
26.2k
      sub_ddmmss (bh, bl, bh, bl, ah, al);
81
26.2k
      if (bh < 2)
82
530
  return 0;
83
84
25.6k
      u00 = u10 = u11 = 1;
85
25.6k
      u01 = 0;
86
25.6k
    }
87
88
49.8k
  if (ah < bh)
89
24.7k
    goto subtract_a;
90
91
25.1k
  for (;;)
92
473k
    {
93
473k
      ASSERT (ah >= bh);
94
473k
      if (ah == bh)
95
174
  goto done;
96
97
473k
      if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2)))
98
24.7k
  {
99
24.7k
    ah = (ah << (GMP_LIMB_BITS / 2) ) + (al >> (GMP_LIMB_BITS / 2));
100
24.7k
    bh = (bh << (GMP_LIMB_BITS / 2) ) + (bl >> (GMP_LIMB_BITS / 2));
101
102
24.7k
    break;
103
24.7k
  }
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
448k
      ASSERT (ah > bh);
108
448k
      sub_ddmmss (ah, al, ah, al, bh, bl);
109
110
448k
      if (ah < 2)
111
12
  goto done;
112
113
448k
      if (ah <= bh)
114
177k
  {
115
    /* Use q = 1 */
116
177k
    u01 += u00;
117
177k
    u11 += u10;
118
177k
  }
119
270k
      else
120
270k
  {
121
270k
    mp_limb_t r[2];
122
270k
    mp_limb_t q = div2 (r, ah, al, bh, bl);
123
270k
    al = r[0]; ah = r[1];
124
270k
    if (ah < 2)
125
167
      {
126
        /* A is too small, but q is correct. */
127
167
        u01 += q * u00;
128
167
        u11 += q * u10;
129
167
        goto done;
130
167
      }
131
270k
    q++;
132
270k
    u01 += q * u00;
133
270k
    u11 += q * u10;
134
270k
  }
135
473k
    subtract_a:
136
473k
      ASSERT (bh >= ah);
137
473k
      if (ah == bh)
138
143
  goto done;
139
140
473k
      if (bh < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2)))
141
24.5k
  {
142
24.5k
    ah = (ah << (GMP_LIMB_BITS / 2) ) + (al >> (GMP_LIMB_BITS / 2));
143
24.5k
    bh = (bh << (GMP_LIMB_BITS / 2) ) + (bl >> (GMP_LIMB_BITS / 2));
144
145
24.5k
    goto subtract_a1;
146
24.5k
  }
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
448k
      sub_ddmmss (bh, bl, bh, bl, ah, al);
151
152
448k
      if (bh < 2)
153
6
  goto done;
154
155
448k
      if (bh <= ah)
156
179k
  {
157
    /* Use q = 1 */
158
179k
    u00 += u01;
159
179k
    u10 += u11;
160
179k
  }
161
269k
      else
162
269k
  {
163
269k
    mp_limb_t r[2];
164
269k
    mp_limb_t q = div2 (r, bh, bl, ah, al);
165
269k
    bl = r[0]; bh = r[1];
166
269k
    if (bh < 2)
167
104
      {
168
        /* B is too small, but q is correct. */
169
104
        u00 += q * u01;
170
104
        u10 += q * u11;
171
104
        goto done;
172
104
      }
173
269k
    q++;
174
269k
    u00 += q * u01;
175
269k
    u10 += q * u11;
176
269k
  }
177
448k
    }
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
24.7k
  for (;;)
183
418k
    {
184
418k
      ASSERT (ah >= bh);
185
186
418k
      ah -= bh;
187
418k
      if (ah < (CNST_LIMB (1) << (GMP_LIMB_BITS / 2 + 1)))
188
11.0k
  break;
189
190
407k
      if (ah <= bh)
191
170k
  {
192
    /* Use q = 1 */
193
170k
    u01 += u00;
194
170k
    u11 += u10;
195
170k
  }
196
236k
      else
197
236k
  {
198
236k
    mp_double_limb_t rq = div1 (ah, bh);
199
236k
    mp_limb_t q = rq.d1;
200
236k
    ah = rq.d0;
201
202
236k
    if (ah < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1)))
203
13.3k
      {
204
        /* A is too small, but q is correct. */
205
13.3k
        u01 += q * u00;
206
13.3k
        u11 += q * u10;
207
13.3k
        break;
208
13.3k
      }
209
223k
    q++;
210
223k
    u01 += q * u00;
211
223k
    u11 += q * u10;
212
223k
  }
213
418k
    subtract_a1:
214
418k
      ASSERT (bh >= ah);
215
216
418k
      bh -= ah;
217
418k
      if (bh < (CNST_LIMB (1) << (GMP_LIMB_BITS / 2 + 1)))
218
11.2k
  break;
219
220
407k
      if (bh <= ah)
221
170k
  {
222
    /* Use q = 1 */
223
170k
    u00 += u01;
224
170k
    u10 += u11;
225
170k
  }
226
236k
      else
227
236k
  {
228
236k
    mp_double_limb_t rq = div1 (bh, ah);
229
236k
    mp_limb_t q = rq.d1;
230
236k
    bh = rq.d0;
231
232
236k
    if (bh < (CNST_LIMB(1) << (GMP_LIMB_BITS / 2 + 1)))
233
13.5k
      {
234
        /* B is too small, but q is correct. */
235
13.5k
        u00 += q * u01;
236
13.5k
        u10 += q * u11;
237
13.5k
        break;
238
13.5k
      }
239
223k
    q++;
240
223k
    u00 += q * u01;
241
223k
    u10 += q * u11;
242
223k
  }
243
407k
    }
244
245
49.8k
 done:
246
49.8k
  M->u[0][0] = u00; M->u[0][1] = u01;
247
49.8k
  M->u[1][0] = u10; M->u[1][1] = u11;
248
249
49.8k
  return 1;
250
24.7k
}
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
41.7k
{
258
41.7k
  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
41.7k
#if HAVE_NATIVE_mpn_addaddmul_1msb0
269
41.7k
  ah = mpn_addaddmul_1msb0 (rp, ap, bp, n, M->u[0][0], M->u[1][0]);
270
41.7k
  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
41.7k
  rp[n] = ah;
279
41.7k
  bp[n] = bh;
280
281
41.7k
  n += (ah | bh) > 0;
282
41.7k
  return n;
283
41.7k
}