Commit | Line | Data |
---|---|---|
65397be7 | 1 | /* -*-c-*- |
2 | * | |
65397be7 | 3 | * Compute integer square roots |
4 | * | |
5 | * (c) 2000 Straylight/Edgeware | |
6 | */ | |
7 | ||
45c0fd36 | 8 | /*----- Licensing notice --------------------------------------------------* |
65397be7 | 9 | * |
10 | * This file is part of Catacomb. | |
11 | * | |
12 | * Catacomb is free software; you can redistribute it and/or modify | |
13 | * it under the terms of the GNU Library General Public License as | |
14 | * published by the Free Software Foundation; either version 2 of the | |
15 | * License, or (at your option) any later version. | |
45c0fd36 | 16 | * |
65397be7 | 17 | * Catacomb is distributed in the hope that it will be useful, |
18 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
19 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
20 | * GNU Library General Public License for more details. | |
45c0fd36 | 21 | * |
65397be7 | 22 | * You should have received a copy of the GNU Library General Public |
23 | * License along with Catacomb; if not, write to the Free | |
24 | * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, | |
25 | * MA 02111-1307, USA. | |
26 | */ | |
27 | ||
65397be7 | 28 | /*----- Header files ------------------------------------------------------*/ |
29 | ||
30 | #include "mp.h" | |
31 | ||
32 | /*----- Main code ---------------------------------------------------------*/ | |
33 | ||
34 | /* --- @mp_sqrt@ --- * | |
35 | * | |
36 | * Arguments: @mp *d@ = pointer to destination integer | |
37 | * @mp *a@ = (nonnegative) integer to take square root of | |
38 | * | |
39 | * Returns: The largest integer %$x$% such that %$x^2 \le a$%. | |
40 | * | |
41 | * Use: Computes integer square roots. | |
42 | * | |
43 | * The current implementation isn't very good: it uses the | |
44 | * Newton-Raphson method to find an approximation to %$a$%. If | |
45 | * there's any demand for a better version, I'll write one. | |
46 | */ | |
47 | ||
48 | mp *mp_sqrt(mp *d, mp *a) | |
49 | { | |
50 | unsigned long z; | |
51 | mp *q = MP_NEW, *r = MP_NEW; | |
52 | ||
53 | /* --- Sanity preservation --- */ | |
54 | ||
a69a3efd | 55 | assert(!MP_NEGP(a)); |
65397be7 | 56 | |
57 | /* --- Deal with trivial cases --- */ | |
58 | ||
59 | MP_SHRINK(a); | |
0c9d150d | 60 | if (MP_ZEROP(a)) { |
f1140c41 | 61 | mp_drop(d); |
65397be7 | 62 | return (MP_ZERO); |
63 | } | |
64 | ||
65 | /* --- Find an initial guess of about the right size --- */ | |
66 | ||
67 | z = mp_bits(a); | |
68 | z >>= 1; | |
69 | mp_copy(a); | |
70 | d = mp_lsr(d, a, z); | |
65397be7 | 71 | |
72 | /* --- Main approximation --- * | |
73 | * | |
d9383a70 MW |
74 | * The Newton--Raphson method finds approximate zeroes of a function by |
75 | * starting with a guess and repeatedly refining the guess by approximating | |
76 | * the function near the guess by its tangent at the guess | |
77 | * %$x$%-coordinate, using where the tangent cuts the %$x$%-axis as the new | |
78 | * guess. | |
65397be7 | 79 | * |
d9383a70 MW |
80 | * Given a function %$f(x)$% and a guess %$x_i$%, the tangent has the |
81 | * equation %$y = f(x_i) + f'(x_i) (x - x_i)$%, which is zero when | |
65397be7 | 82 | * |
d9383a70 MW |
83 | * %$\displaystyle x = x_i - \frac{f(x_i)}{f'(x_i)} |
84 | * | |
85 | * We set %$f(x) = x^2 - a$%, so our recurrence will be | |
86 | * | |
87 | * %$\displaystyle x_{i+1} = x_i - \frac{x_i^2 - a}{2 x_i}$% | |
88 | * | |
89 | * It's possible to simplify this, but it's useful to see %$q = x_i^2 - a$% | |
90 | * so that we know when to stop. We want the largest integer not larger | |
91 | * than the true square root. If %$q > 0$% then %$x_i$% is definitely too | |
92 | * large, and we should decrease it by at least one even if the adjustment | |
93 | * term %$(x_i^2 - a)/2 x$% is less than one. | |
94 | * | |
95 | * Suppose, then, that %$q \le 0$%. Then %$(x_i + 1)^2 - a = {}$% | |
96 | * $%x_i^2 + 2 x_i + 1 - a = q + 2 x_i + 1$%. Hence, if %$q \ge -2 x_i$% | |
97 | * then %$x_i + 1$% is an overestimate and we should settle where we are. | |
98 | * Otherwise, %$x_i + 1$% is an underestimate -- but, in this case the | |
99 | * adjustment will always be at least one. | |
65397be7 | 100 | */ |
101 | ||
102 | for (;;) { | |
103 | q = mp_sqr(q, d); | |
104 | q = mp_sub(q, q, a); | |
0c9d150d | 105 | if (MP_ZEROP(q)) |
65397be7 | 106 | break; |
a69a3efd | 107 | if (MP_NEGP(q)) { |
65397be7 | 108 | r = mp_lsl(r, d, 1); |
109 | r->f |= MP_NEG; | |
0c9d150d | 110 | if (MP_CMP(q, >=, r)) |
65397be7 | 111 | break; |
112 | } | |
113 | mp_div(&r, &q, q, d); | |
114 | r = mp_lsr(r, r, 1); | |
115 | if (r->v == r->vl) | |
116 | d = mp_sub(d, d, MP_ONE); | |
117 | else | |
118 | d = mp_sub(d, d, r); | |
119 | } | |
120 | ||
121 | /* --- Finished, at last --- */ | |
122 | ||
432c4e18 | 123 | mp_drop(a); |
65397be7 | 124 | mp_drop(q); |
f1140c41 | 125 | mp_drop(r); |
65397be7 | 126 | return (d); |
127 | } | |
128 | ||
129 | /*----- Test rig ----------------------------------------------------------*/ | |
130 | ||
131 | #ifdef TEST_RIG | |
132 | ||
133 | #include <mLib/testrig.h> | |
134 | ||
135 | static int verify(dstr *v) | |
136 | { | |
137 | mp *a = *(mp **)v[0].buf; | |
138 | mp *qq = *(mp **)v[1].buf; | |
139 | mp *q = mp_sqrt(MP_NEW, a); | |
140 | int ok = 1; | |
141 | ||
4b536f42 | 142 | if (!MP_EQ(q, qq)) { |
65397be7 | 143 | ok = 0; |
144 | fputs("\n*** sqrt failed", stderr); | |
45c0fd36 | 145 | fputs("\n*** a = ", stderr); mp_writefile(a, stderr, 10); |
65397be7 | 146 | fputs("\n*** result = ", stderr); mp_writefile(q, stderr, 10); |
147 | fputs("\n*** expect = ", stderr); mp_writefile(qq, stderr, 10); | |
148 | fputc('\n', stderr); | |
149 | } | |
150 | ||
151 | mp_drop(a); | |
152 | mp_drop(q); | |
153 | mp_drop(qq); | |
154 | assert(mparena_count(MPARENA_GLOBAL) == 0); | |
155 | ||
156 | return (ok); | |
157 | } | |
158 | ||
159 | static test_chunk tests[] = { | |
160 | { "sqrt", verify, { &type_mp, &type_mp, 0 } }, | |
161 | { 0, 0, { 0 } }, | |
162 | }; | |
163 | ||
164 | int main(int argc, char *argv[]) | |
165 | { | |
166 | sub_init(); | |
0f00dc4c | 167 | test_run(argc, argv, tests, SRCDIR "/t/mp"); |
65397be7 | 168 | return (0); |
169 | } | |
170 | ||
171 | #endif | |
172 | ||
173 | /*----- That's all, folks -------------------------------------------------*/ |