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
/*
Copyright (C) 2009, 2015 William Hart
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 "ulong_extras.h"
ulong
n_gcdinv(ulong * s, ulong x, ulong y)
{
slong v1, v2, t2;
ulong d, r, quot, rem;
FLINT_ASSERT(y > x);
v1 = 0;
v2 = 1;
r = x;
x = y;
/* y and x both have top bit set */
if ((slong) (x & r) < 0)
{
d = x - r;
t2 = v2;
x = r;
v2 = v1 - v2;
v1 = t2;
r = d;
}
/* second value has second msb set */
while ((slong) (r << 1) < 0)
{
d = x - r;
if (d < r) /* quot = 1 */
{
t2 = v2;
x = r;
v2 = v1 - v2;
v1 = t2;
r = d;
}
else if (d < (r << 1)) /* quot = 2 */
{
x = r;
t2 = v2;
v2 = v1 - ((ulong) v2 << 1);
v1 = t2;
r = d - x;
}
else /* quot = 3 */
{
x = r;
t2 = v2;
v2 = v1 - 3 * v2;
v1 = t2;
r = d - (x << 1);
}
}
while (r)
{
/* overflow not possible due to top 2 bits of r not being set */
if (x < (r << 2)) /* if quot < 4 */
{
d = x - r;
if (d < r) /* quot = 1 */
{
t2 = v2;
x = r;
v2 = v1 - v2;
v1 = t2;
r = d;
}
else if (d < (r << 1)) /* quot = 2 */
{
x = r;
t2 = v2;
v2 = v1 - ((ulong) v2 << 1);
v1 = t2;
r = d - x;
}
else /* quot = 3 */
{
x = r;
t2 = v2;
v2 = v1 - 3 * (ulong) v2;
v1 = t2;
r = d - (x << 1);
}
}
else
{
quot = x / r;
rem = x - r * quot;
x = r;
t2 = v2;
v2 = v1 - quot * v2;
v1 = t2;
r = rem;
}
}
if (v1 < WORD(0))
v1 += y;
(*s) = v1;
return x;
}