summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/rec_sqrt.c
blob: ef686ae4e59fd5c994123470b6513f7de182fd56 (plain)
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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
/* mpfr_rec_sqrt -- inverse square root

Copyright 2008-2018 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.

This file is part of the GNU MPFR Library.

The GNU MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

The GNU MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */

#define MPFR_NEED_LONGLONG_H /* for umul_ppmm */
#include "mpfr-impl.h"

#define LIMB_SIZE(x) ((((x)-1)>>MPFR_LOG2_GMP_NUMB_BITS) + 1)

#define MPFR_COM_N(x,y,n)                               \
  do                                                    \
    {                                                   \
      mp_size_t i;                                      \
      for (i = 0; i < n; i++)                           \
        *((x)+i) = ~*((y)+i);                           \
    }                                                   \
  while (0)

/* Put in X a p-bit approximation of 1/sqrt(A),
   where X = {x, n}/B^n, n = ceil(p/GMP_NUMB_BITS),
   A = 2^(1+as)*{a, an}/B^an, as is 0 or 1, an = ceil(ap/GMP_NUMB_BITS),
   where B = 2^GMP_NUMB_BITS.

   We have 1 <= A < 4 and 1/2 <= X < 1.

   The error in the approximate result with respect to the true
   value 1/sqrt(A) is bounded by 1 ulp(X), i.e., 2^{-p} since 1/2 <= X < 1.

   Note: x and a are left-aligned, i.e., the most significant bit of
   a[an-1] is set, and so is the most significant bit of the output x[n-1].

   If p is not a multiple of GMP_NUMB_BITS, the extra low bits of the input
   A are taken into account to compute the approximation of 1/sqrt(A), but
   whether or not they are zero, the error between X and 1/sqrt(A) is bounded
   by 1 ulp(X) [in precision p].
   The extra low bits of the output X (if p is not a multiple of GMP_NUMB_BITS)
   are set to 0.

   Assumptions:
   (1) A should be normalized, i.e., the most significant bit of a[an-1]
       should be 1. If as=0, we have 1 <= A < 2; if as=1, we have 2 <= A < 4.
   (2) p >= 12
   (3) {a, an} and {x, n} should not overlap
   (4) GMP_NUMB_BITS >= 12 and is even

   Note: this routine is much more efficient when ap is small compared to p,
   including the case where ap <= GMP_NUMB_BITS, thus it can be used to
   implement an efficient mpfr_rec_sqrt_ui function.

   References:
   [1] Modern Computer Algebra, Richard Brent and Paul Zimmermann,
   https://members.loria.fr/PZimmermann/mca/pub226.html
*/
static void
mpfr_mpn_rec_sqrt (mpfr_limb_ptr x, mpfr_prec_t p,
                   mpfr_limb_srcptr a, mpfr_prec_t ap, int as)

{
  /* the following T1 and T2 are bipartite tables giving initial
     approximation for the inverse square root, with 13-bit input split in
     5+4+4, and 11-bit output. More precisely, if 2048 <= i < 8192,
     with i = a*2^8 + b*2^4 + c, we use for approximation of
     2048/sqrt(i/2048) the value x = T1[16*(a-8)+b] + T2[16*(a-8)+c].
     The largest error is obtained for i = 2054, where x = 2044,
     and 2048/sqrt(i/2048) = 2045.006576...
  */
  static short int T1[384] = {
2040, 2033, 2025, 2017, 2009, 2002, 1994, 1987, 1980, 1972, 1965, 1958, 1951,
1944, 1938, 1931, /* a=8 */
1925, 1918, 1912, 1905, 1899, 1892, 1886, 1880, 1874, 1867, 1861, 1855, 1849,
1844, 1838, 1832, /* a=9 */
1827, 1821, 1815, 1810, 1804, 1799, 1793, 1788, 1783, 1777, 1772, 1767, 1762,
1757, 1752, 1747, /* a=10 */
1742, 1737, 1733, 1728, 1723, 1718, 1713, 1709, 1704, 1699, 1695, 1690, 1686,
1681, 1677, 1673, /* a=11 */
1669, 1664, 1660, 1656, 1652, 1647, 1643, 1639, 1635, 1631, 1627, 1623, 1619,
1615, 1611, 1607, /* a=12 */
1603, 1600, 1596, 1592, 1588, 1585, 1581, 1577, 1574, 1570, 1566, 1563, 1559,
1556, 1552, 1549, /* a=13 */
1545, 1542, 1538, 1535, 1532, 1528, 1525, 1522, 1518, 1515, 1512, 1509, 1505,
1502, 1499, 1496, /* a=14 */
1493, 1490, 1487, 1484, 1481, 1478, 1475, 1472, 1469, 1466, 1463, 1460, 1457,
1454, 1451, 1449, /* a=15 */
1446, 1443, 1440, 1438, 1435, 1432, 1429, 1427, 1424, 1421, 1419, 1416, 1413,
1411, 1408, 1405, /* a=16 */
1403, 1400, 1398, 1395, 1393, 1390, 1388, 1385, 1383, 1380, 1378, 1375, 1373,
1371, 1368, 1366, /* a=17 */
1363, 1360, 1358, 1356, 1353, 1351, 1349, 1346, 1344, 1342, 1340, 1337, 1335,
1333, 1331, 1329, /* a=18 */
1327, 1325, 1323, 1321, 1319, 1316, 1314, 1312, 1310, 1308, 1306, 1304, 1302,
1300, 1298, 1296, /* a=19 */
1294, 1292, 1290, 1288, 1286, 1284, 1282, 1280, 1278, 1276, 1274, 1272, 1270,
1268, 1266, 1265, /* a=20 */
1263, 1261, 1259, 1257, 1255, 1253, 1251, 1250, 1248, 1246, 1244, 1242, 1241,
1239, 1237, 1235, /* a=21 */
1234, 1232, 1230, 1229, 1227, 1225, 1223, 1222, 1220, 1218, 1217, 1215, 1213,
1212, 1210, 1208, /* a=22 */
1206, 1204, 1203, 1201, 1199, 1198, 1196, 1195, 1193, 1191, 1190, 1188, 1187,
1185, 1184, 1182, /* a=23 */
1181, 1180, 1178, 1177, 1175, 1174, 1172, 1171, 1169, 1168, 1166, 1165, 1163,
1162, 1160, 1159, /* a=24 */
1157, 1156, 1154, 1153, 1151, 1150, 1149, 1147, 1146, 1144, 1143, 1142, 1140,
1139, 1137, 1136, /* a=25 */
1135, 1133, 1132, 1131, 1129, 1128, 1127, 1125, 1124, 1123, 1121, 1120, 1119,
1117, 1116, 1115, /* a=26 */
1114, 1113, 1111, 1110, 1109, 1108, 1106, 1105, 1104, 1103, 1101, 1100, 1099,
1098, 1096, 1095, /* a=27 */
1093, 1092, 1091, 1090, 1089, 1087, 1086, 1085, 1084, 1083, 1081, 1080, 1079,
1078, 1077, 1076, /* a=28 */
1075, 1073, 1072, 1071, 1070, 1069, 1068, 1067, 1065, 1064, 1063, 1062, 1061,
1060, 1059, 1058, /* a=29 */
1057, 1056, 1055, 1054, 1052, 1051, 1050, 1049, 1048, 1047, 1046, 1045, 1044,
1043, 1042, 1041, /* a=30 */
1040, 1039, 1038, 1037, 1036, 1035, 1034, 1033, 1032, 1031, 1030, 1029, 1028,
1027, 1026, 1025 /* a=31 */
};
  static unsigned char T2[384] = {
    7, 7, 6, 6, 5, 5, 4, 4, 4, 3, 3, 2, 2, 1, 1, 0, /* a=8 */
    6, 5, 5, 5, 4, 4, 3, 3, 3, 2, 2, 2, 1, 1, 0, 0, /* a=9 */
    5, 5, 4, 4, 4, 3, 3, 3, 2, 2, 2, 1, 1, 1, 0, 0, /* a=10 */
    4, 4, 3, 3, 3, 3, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, /* a=11 */
    3, 3, 3, 3, 2, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, 0, /* a=12 */
    3, 3, 3, 2, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, 0, 0, /* a=13 */
    3, 3, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 0, 0, 0, 0, /* a=14 */
    2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, /* a=15 */
    2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, /* a=16 */
    2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, /* a=17 */
    3, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 0, /* a=18 */
    2, 2, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, /* a=19 */
    1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, /* a=20 */
    1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, /* a=21 */
    1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, /* a=22 */
    2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, /* a=23 */
    1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, /* a=24 */
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, /* a=25 */
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, /* a=26 */
    1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, /* a=27 */
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, /* a=28 */
    1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, /* a=29 */
    1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, /* a=30 */
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0  /* a=31 */
};
  mp_size_t n = LIMB_SIZE(p);   /* number of limbs of X */
  mp_size_t an = LIMB_SIZE(ap); /* number of limbs of A */

  /* A should be normalized */
  MPFR_ASSERTD((a[an - 1] & MPFR_LIMB_HIGHBIT) != 0);
  /* We should have enough bits in one limb and GMP_NUMB_BITS should be even.
     Since that does not depend on MPFR, we always check this. */
  MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS >= 12 && (GMP_NUMB_BITS & 1) == 0);
  /* {a, an} and {x, n} should not overlap */
  MPFR_ASSERTD((a + an <= x) || (x + n <= a));
  MPFR_ASSERTD(p >= 11);

  if (MPFR_UNLIKELY(an > n)) /* we can cut the input to n limbs */
    {
      a += an - n;
      an = n;
    }

  if (p == 11) /* should happen only from recursive calls */
    {
      unsigned long i, ab, ac;
      mp_limb_t t;

      /* take the 12+as most significant bits of A */
      i = a[an - 1] >> (GMP_NUMB_BITS - (12 + as));
      /* if one wants faithful rounding for p=11, replace #if 0 by #if 1 */
      ab = i >> 4;
      ac = (ab & 0x3F0) | (i & 0x0F);
      t = (mp_limb_t) T1[ab - 0x80] + (mp_limb_t) T2[ac - 0x80];
      x[0] = t << (GMP_NUMB_BITS - p);
    }
  else /* p >= 12 */
    {
      mpfr_prec_t h, pl;
      mpfr_limb_ptr r, s, t, u;
      mp_size_t xn, rn, th, ln, tn, sn, ahn, un;
      mp_limb_t neg, cy, cu;
      MPFR_TMP_DECL(marker);

      /* compared to Algorithm 3.9 of [1], we have {a, an} = A/2 if as=0,
         and A/4 if as=1. */

      /* h = max(11, ceil((p+3)/2)) is the bitsize of the recursive call */
      h = (p < 18) ? 11 : (p >> 1) + 2;

      xn = LIMB_SIZE(h);       /* limb size of the recursive Xh */
      rn = LIMB_SIZE(2 * h);   /* a priori limb size of Xh^2 */
      ln = n - xn;             /* remaining limbs to be computed */

      /* Since |Xh - A^{-1/2}| <= 2^{-h}, then by multiplying by Xh + A^{-1/2}
         we get |Xh^2 - 1/A| <= 2^{-h+1}, thus |A*Xh^2 - 1| <= 2^{-h+3},
         thus the h-3 most significant bits of t should be zero,
         which is in fact h+1+as-3 because of the normalization of A.
         This corresponds to th=floor((h+1+as-3)/GMP_NUMB_BITS) limbs.

         More precisely we have |Xh^2 - 1/A| <= 2^{-h} * (Xh + A^{-1/2})
         <= 2^{-h} * (2 A^{-1/2} + 2^{-h}) <= 2.001 * 2^{-h} * A^{-1/2}
         since A < 4 and h >= 11, thus
         |A*Xh^2 - 1| <= 2.001 * 2^{-h} * A^{1/2} <= 1.001 * 2^{2-h}.
         This is sufficient to prove that the upper limb of {t,tn} below is
         less that 0.501 * 2^GMP_NUMB_BITS, thus cu = 0 below.
      */
      th = (h + 1 + as - 3) >> MPFR_LOG2_GMP_NUMB_BITS;
      tn = LIMB_SIZE(2 * h + 1 + as);

      /* we need h+1+as bits of a */
      ahn = LIMB_SIZE(h + 1 + as); /* number of high limbs of A
                                      needed for the recursive call*/
      if (MPFR_UNLIKELY(ahn > an))
        ahn = an;
      mpfr_mpn_rec_sqrt (x + ln, h, a + an - ahn, ahn * GMP_NUMB_BITS, as);
      /* the most h significant bits of X are set, X has ceil(h/GMP_NUMB_BITS)
         limbs, the low (-h) % GMP_NUMB_BITS bits are zero */

      /* compared to Algorithm 3.9 of [1], we have {x+ln,xn} = X_h */

      MPFR_TMP_MARK (marker);
      /* first step: square X in r, result is exact */
      un = xn + (tn - th);
      /* We use the same temporary buffer to store r and u: r needs 2*xn
         limbs where u needs xn+(tn-th) limbs. Since tn can store at least
         2h bits, and th at most h bits, then tn-th can store at least h bits,
         thus tn - th >= xn, and reserving the space for u is enough. */
      MPFR_ASSERTD(2 * xn <= un);
      u = r = MPFR_TMP_LIMBS_ALLOC (un);
      if (2 * h <= GMP_NUMB_BITS) /* xn=rn=1, and since p <= 2h-3, n=1,
                                     thus ln = 0 */
        {
          MPFR_ASSERTD(ln == 0);
          cy = x[0] >> (GMP_NUMB_BITS >> 1);
          r ++;
          r[0] = cy * cy;
        }
      else if (xn == 1) /* xn=1, rn=2 */
        umul_ppmm(r[1], r[0], x[ln], x[ln]);
      else
        {
          mpn_mul_n (r, x + ln, x + ln, xn);
          /* we have {r, 2*xn} = X_h^2 */
          if (rn < 2 * xn)
            r ++;
        }
      /* now the 2h most significant bits of {r, rn} contains X^2, r has rn
         limbs, and the low (-2h) % GMP_NUMB_BITS bits are zero */

      /* Second step: s <- A * (r^2), and truncate the low ap bits,
         i.e., at weight 2^{-2h} (s is aligned to the low significant bits)
       */
      sn = an + rn;
      s = MPFR_TMP_LIMBS_ALLOC (sn);
      if (rn == 1) /* rn=1 implies n=1, since rn*GMP_NUMB_BITS >= 2h,
                           and 2h >= p+3 */
        {
          /* necessarily p <= GMP_NUMB_BITS-3: we can ignore the two low
             bits from A */
          /* since n=1, and we ensured an <= n, we also have an=1 */
          MPFR_ASSERTD(an == 1);
          umul_ppmm (s[1], s[0], r[0], a[0]);
        }
      else
        {
          /* we have p <= n * GMP_NUMB_BITS
             2h <= rn * GMP_NUMB_BITS with p+3 <= 2h <= p+4
             thus n <= rn <= n + 1 */
          MPFR_ASSERTD(rn <= n + 1);
          /* since we ensured an <= n, we have an <= rn */
          MPFR_ASSERTD(an <= rn);
          mpn_mul (s, r, rn, a, an);
          /* s should be near B^sn/2^(1+as), thus s[sn-1] is either
             100000... or 011111... if as=0, or
             010000... or 001111... if as=1.
             We ignore the bits of s after the first 2h+1+as ones.
             We have {s, rn+an} = A*X_h^2/2 if as=0, A*X_h^2/4 if as=1. */
        }

      /* We ignore the bits of s after the first 2h+1+as ones: s has rn + an
         limbs, where rn = LIMBS(2h), an=LIMBS(a), and tn = LIMBS(2h+1+as). */
      t = s + sn - tn; /* pointer to low limb of the high part of t */
      /* the upper h-3 bits of 1-t should be zero,
         where 1 corresponds to the most significant bit of t[tn-1] if as=0,
         and to the 2nd most significant bit of t[tn-1] if as=1 */

      /* compute t <- 1 - t, which is B^tn - {t, tn+1},
         with rounding toward -Inf, i.e., rounding the input t toward +Inf.
         We could only modify the low tn - th limbs from t, but it gives only
         a small speedup, and would make the code more complex.
      */
      neg = t[tn - 1] & (MPFR_LIMB_HIGHBIT >> as);
      if (neg == 0) /* Ax^2 < 1: we have t = th + eps, where 0 <= eps < ulp(th)
                       is the part truncated above, thus 1 - t rounded to -Inf
                       is 1 - th - ulp(th) */
        {
          /* since the 1+as most significant bits of t are zero, set them
             to 1 before the one-complement */
          t[tn - 1] |= MPFR_LIMB_HIGHBIT | (MPFR_LIMB_HIGHBIT >> as);
          MPFR_COM_N (t, t, tn);
          /* we should add 1 here to get 1-th complement, and subtract 1 for
             -ulp(th), thus we do nothing */
        }
      else /* negative case: we want 1 - t rounded toward -Inf, i.e.,
              th + eps rounded toward +Inf, which is th + ulp(th):
              we discard the bit corresponding to 1,
              and we add 1 to the least significant bit of t */
        {
          t[tn - 1] ^= neg;
          mpn_add_1 (t, t, tn, 1);
        }
      tn -= th; /* we know at least th = floor((h+1+as-3)/GMP_NUMB_LIMBS) of
                   the high limbs of {t, tn} are zero */

      /* tn = rn - th, where rn * GMP_NUMB_BITS >= 2*h and
         th * GMP_NUMB_BITS <= h+1+as-3, thus tn > 0 */
      MPFR_ASSERTD(tn > 0);

      /* u <- x * t, where {t, tn} contains at least h+3 bits,
         and {x, xn} contains h bits, thus tn >= xn */
      MPFR_ASSERTD(tn >= xn);
      if (tn == 1) /* necessarily xn=1 */
        umul_ppmm (u[1], u[0], t[0], x[ln]);
      else
        mpn_mul (u, t, tn, x + ln, xn);

      /* we have {u, tn+xn} = T_l X_h/2 if as=0, T_l X_h/4 if as=1 */

      /* we have already discarded the upper th high limbs of t, thus we only
         have to consider the upper n - th limbs of u */
      un = n - th; /* un cannot be zero, since p <= n*GMP_NUMB_BITS,
                      h = ceil((p+3)/2) <= (p+4)/2,
                      th*GMP_NUMB_BITS <= h-1 <= p/2+1,
                      thus (n-th)*GMP_NUMB_BITS >= p/2-1.
                   */
      MPFR_ASSERTD(un > 0);
      u += (tn + xn) - un; /* xn + tn - un = xn + (original_tn - th) - (n - th)
                                           = xn + original_tn - n
                              = LIMBS(h) + LIMBS(2h+1+as) - LIMBS(p) > 0
                              since 2h >= p+3 */
      MPFR_ASSERTD(tn + xn > un); /* will allow to access u[-1] below */

      /* In case as=0, u contains |x*(1-Ax^2)/2|, which is exactly what we
         need to add or subtract.
         In case as=1, u contains |x*(1-Ax^2)/4|, thus we need to multiply
         u by 2. */

      if (as == 1)
        /* shift on un+1 limbs to get most significant bit of u[-1] into
           least significant bit of u[0] */
        mpn_lshift (u - 1, u - 1, un + 1, 1);

      /* now {u,un} represents U / 2 from Algorithm 3.9 */

      pl = n * GMP_NUMB_BITS - p;       /* low bits from x */
      /* We want that the low pl bits are zero after rounding to nearest,
         thus we round u to nearest at bit pl-1 of u[0] */
      if (pl > 0)
        {
          cu = mpn_add_1 (u, u, un, u[0] & (MPFR_LIMB_ONE << (pl - 1)));
          /* mask bits 0..pl-1 of u[0] */
          u[0] &= ~MPFR_LIMB_MASK(pl);
        }
      else /* round bit is in u[-1] */
        cu = mpn_add_1 (u, u, un, u[-1] >> (GMP_NUMB_BITS - 1));
      MPFR_ASSERTN(cu == 0);

      /* We already have filled {x + ln, xn = n - ln}, and we want to add or
         subtract {u, un} at position x.
         un = n - th, where th contains <= h+1+as-3<=h-1 bits
         ln = n - xn, where xn contains >= h bits
         thus un > ln.
         Warning: ln might be zero.
      */
      MPFR_ASSERTD(un > ln);
      /* we can have un = ln + 2, for example with GMP_NUMB_BITS=32 and
         p=62, as=0, then h=33, n=2, th=0, xn=2, thus un=2 and ln=0. */
      MPFR_ASSERTD(un == ln + 1 || un == ln + 2);
      /* the high un-ln limbs of u will overlap the low part of {x+ln,xn},
         we need to add or subtract the overlapping part {u + ln, un - ln} */
      /* Warning! th may be 0, in which case the mpn_add_1 and mpn_sub_1
         below (with size = th) mustn't be used. */
      if (neg == 0)
        {
          if (ln > 0)
            MPN_COPY (x, u, ln);
          cy = mpn_add (x + ln, x + ln, xn, u + ln, un - ln);
          /* cy is the carry at x + (ln + xn) = x + n */
        }
      else /* negative case */
        {
          /* subtract {u+ln, un-ln} from {x+ln,un} */
          cy = mpn_sub (x + ln, x + ln, xn, u + ln, un - ln);
          /* cy is the borrow at x + (ln + xn) = x + n */

          /* cy cannot be non-zero, since the most significant bit of Xh is 1,
             and the correction is bounded by 2^{-h+3} */
          MPFR_ASSERTD(cy == 0);
          if (ln > 0)
            {
              MPFR_COM_N (x, u, ln);
              /* we must add one for the 2-complement ... */
              cy = mpn_add_1 (x, x, n, MPFR_LIMB_ONE);
              /* ... and subtract 1 at x[ln], where n = ln + xn */
              cy -= mpn_sub_1 (x + ln, x + ln, xn, MPFR_LIMB_ONE);
            }
        }

      /* cy can be 1 when A=1, i.e., {a, n} = B^n. In that case we should
         have X = B^n, and setting X to 1-2^{-p} satisfies the error bound
         of 1 ulp. */
      if (MPFR_UNLIKELY(cy != 0))
        {
          cy -= mpn_sub_1 (x, x, n, MPFR_LIMB_ONE << pl);
          MPFR_ASSERTD(cy == 0);
        }

      MPFR_TMP_FREE (marker);
    }
}

int
mpfr_rec_sqrt (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
{
  mpfr_prec_t rp, up, wp;
  mp_size_t rn, wn;
  int s, cy, inex;
  mpfr_limb_ptr x;
  MPFR_TMP_DECL(marker);

  MPFR_LOG_FUNC
    (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (u), mpfr_log_prec, u, rnd_mode),
     ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (r), mpfr_log_prec, r, inex));

  /* special values */
  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(u)))
    {
      if (MPFR_IS_NAN(u))
        {
          MPFR_SET_NAN(r);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_ZERO(u)) /* 1/sqrt(+0) = 1/sqrt(-0) = +Inf */
        {
          /* +0 or -0 */
          MPFR_SET_INF(r);
          MPFR_SET_POS(r);
          MPFR_SET_DIVBY0 ();
          MPFR_RET(0); /* Inf is exact */
        }
      else
        {
          MPFR_ASSERTD(MPFR_IS_INF(u));
          /* 1/sqrt(-Inf) = NAN */
          if (MPFR_IS_NEG(u))
            {
              MPFR_SET_NAN(r);
              MPFR_RET_NAN;
            }
          /* 1/sqrt(+Inf) = +0 */
          MPFR_SET_POS(r);
          MPFR_SET_ZERO(r);
          MPFR_RET(0);
        }
    }

  /* if u < 0, 1/sqrt(u) is NaN */
  if (MPFR_UNLIKELY(MPFR_IS_NEG(u)))
    {
      MPFR_SET_NAN(r);
      MPFR_RET_NAN;
    }

  MPFR_SET_POS(r);

  rp = MPFR_PREC(r); /* output precision */
  up = MPFR_PREC(u); /* input precision */
  wp = rp + 11;      /* initial working precision */

  /* Let u = U*2^e, where e = EXP(u), and 1/2 <= U < 1.
     If e is even, we compute an approximation of X of (4U)^{-1/2},
     and the result is X*2^(-(e-2)/2) [case s=1].
     If e is odd, we compute an approximation of X of (2U)^{-1/2},
     and the result is X*2^(-(e-1)/2) [case s=0]. */

  /* parity of the exponent of u */
  s = 1 - ((mpfr_uexp_t) MPFR_GET_EXP (u) & 1);

  rn = LIMB_SIZE(rp);

  /* for the first iteration, if rp + 11 fits into rn limbs, we round up
     up to a full limb to maximize the chance of rounding, while avoiding
     to allocate extra space */
  wp = rp + 11;
  if (wp < rn * GMP_NUMB_BITS)
    wp = rn * GMP_NUMB_BITS;
  for (;;)
    {
      MPFR_TMP_MARK (marker);
      wn = LIMB_SIZE(wp);
      if (r == u || wn > rn) /* out of place, i.e., we cannot write to r */
        x = MPFR_TMP_LIMBS_ALLOC (wn);
      else
        x = MPFR_MANT(r);
      mpfr_mpn_rec_sqrt (x, wp, MPFR_MANT(u), up, s);
      /* If the input was not truncated, the error is at most one ulp;
         if the input was truncated, the error is at most two ulps
         (see algorithms.tex). */
      if (MPFR_LIKELY (mpfr_round_p (x, wn, wp - (wp < up),
                                     rp + (rnd_mode == MPFR_RNDN))))
        break;

      /* We detect only now the exact case where u=2^(2e), to avoid
         slowing down the average case. This can happen only when the
         mantissa is exactly 1/2 and the exponent is odd. */
      if (s == 0 && mpfr_cmp_ui_2exp (u, 1, MPFR_EXP(u) - 1) == 0)
        {
          mpfr_prec_t pl = wn * GMP_NUMB_BITS - wp;

          /* we should have x=111...111 */
          mpn_add_1 (x, x, wn, MPFR_LIMB_ONE << pl);
          x[wn - 1] = MPFR_LIMB_HIGHBIT;
          s += 2;
          break; /* go through */
        }
      MPFR_TMP_FREE(marker);

      wp += GMP_NUMB_BITS;
    }
  cy = mpfr_round_raw (MPFR_MANT(r), x, wp, 0, rp, rnd_mode, &inex);
  MPFR_EXP(r) = - (MPFR_EXP(u) - 1 - s) / 2;
  if (MPFR_UNLIKELY(cy != 0))
    {
      MPFR_EXP(r) ++;
      MPFR_MANT(r)[rn - 1] = MPFR_LIMB_HIGHBIT;
    }
  MPFR_TMP_FREE(marker);
  return mpfr_check_range (r, inex, rnd_mode);
}