chiark / gitweb /
math/mpx-mul4-amd64-sse2.S: SSE2 multipliers for AMD64.
[catacomb] / math / mp-sqrt.c
CommitLineData
65397be7 1/* -*-c-*-
65397be7 2 *
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
48mp *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 *
74 * We use the Newton-Raphson recurrence relation
75 *
76 * %$x_{i+1} = x_i - \frac{x_i^2 - a}{2 x_i}$%
77 *
78 * We inspect the term %$q = x^2 - a$% to see when to stop. Increasing
79 * %$x$% is pointless when %$-q < 2 x + 1$%.
80 */
81
82 for (;;) {
83 q = mp_sqr(q, d);
84 q = mp_sub(q, q, a);
0c9d150d 85 if (MP_ZEROP(q))
65397be7 86 break;
a69a3efd 87 if (MP_NEGP(q)) {
65397be7 88 r = mp_lsl(r, d, 1);
89 r->f |= MP_NEG;
0c9d150d 90 if (MP_CMP(q, >=, r))
65397be7 91 break;
92 }
93 mp_div(&r, &q, q, d);
94 r = mp_lsr(r, r, 1);
95 if (r->v == r->vl)
96 d = mp_sub(d, d, MP_ONE);
97 else
98 d = mp_sub(d, d, r);
99 }
100
101 /* --- Finished, at last --- */
102
432c4e18 103 mp_drop(a);
65397be7 104 mp_drop(q);
f1140c41 105 mp_drop(r);
65397be7 106 return (d);
107}
108
109/*----- Test rig ----------------------------------------------------------*/
110
111#ifdef TEST_RIG
112
113#include <mLib/testrig.h>
114
115static int verify(dstr *v)
116{
117 mp *a = *(mp **)v[0].buf;
118 mp *qq = *(mp **)v[1].buf;
119 mp *q = mp_sqrt(MP_NEW, a);
120 int ok = 1;
121
4b536f42 122 if (!MP_EQ(q, qq)) {
65397be7 123 ok = 0;
124 fputs("\n*** sqrt failed", stderr);
45c0fd36 125 fputs("\n*** a = ", stderr); mp_writefile(a, stderr, 10);
65397be7 126 fputs("\n*** result = ", stderr); mp_writefile(q, stderr, 10);
127 fputs("\n*** expect = ", stderr); mp_writefile(qq, stderr, 10);
128 fputc('\n', stderr);
129 }
130
131 mp_drop(a);
132 mp_drop(q);
133 mp_drop(qq);
134 assert(mparena_count(MPARENA_GLOBAL) == 0);
135
136 return (ok);
137}
138
139static test_chunk tests[] = {
140 { "sqrt", verify, { &type_mp, &type_mp, 0 } },
141 { 0, 0, { 0 } },
142};
143
144int main(int argc, char *argv[])
145{
146 sub_init();
0f00dc4c 147 test_run(argc, argv, tests, SRCDIR "/t/mp");
65397be7 148 return (0);
149}
150
151#endif
152
153/*----- That's all, folks -------------------------------------------------*/