[BACK]Return to sqrtrem.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / gmp-2.0.2-ssh-2 / mpn / generic

Annotation of OpenXM/src/kan96xx/gmp-2.0.2-ssh-2/mpn/generic/sqrtrem.c, Revision 1.1.1.1

1.1       takayama    1: /* mpn_sqrtrem (root_ptr, rem_ptr, op_ptr, op_size)
                      2:
                      3:    Write the square root of {OP_PTR, OP_SIZE} at ROOT_PTR.
                      4:    Write the remainder at REM_PTR, if REM_PTR != NULL.
                      5:    Return the size of the remainder.
                      6:    (The size of the root is always half of the size of the operand.)
                      7:
                      8:    OP_PTR and ROOT_PTR may not point to the same object.
                      9:    OP_PTR and REM_PTR may point to the same object.
                     10:
                     11:    If REM_PTR is NULL, only the root is computed and the return value of
                     12:    the function is 0 if OP is a perfect square, and *any* non-zero number
                     13:    otherwise.
                     14:
                     15: Copyright (C) 1993, 1994, 1996 Free Software Foundation, Inc.
                     16:
                     17: This file is part of the GNU MP Library.
                     18:
                     19: The GNU MP Library is free software; you can redistribute it and/or modify
                     20: it under the terms of the GNU Library General Public License as published by
                     21: the Free Software Foundation; either version 2 of the License, or (at your
                     22: option) any later version.
                     23:
                     24: The GNU MP Library is distributed in the hope that it will be useful, but
                     25: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
                     26: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
                     27: License for more details.
                     28:
                     29: You should have received a copy of the GNU Library General Public License
                     30: along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
                     31: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
                     32: MA 02111-1307, USA. */
                     33:
                     34: /* This code is just correct if "unsigned char" has at least 8 bits.  It
                     35:    doesn't help to use CHAR_BIT from limits.h, as the real problem is
                     36:    the static arrays.  */
                     37:
                     38: #include "gmp.h"
                     39: #include "gmp-impl.h"
                     40: #include "longlong.h"
                     41:
                     42: /* Square root algorithm:
                     43:
                     44:    1. Shift OP (the input) to the left an even number of bits s.t. there
                     45:       are an even number of words and either (or both) of the most
                     46:       significant bits are set.  This way, sqrt(OP) has exactly half as
                     47:       many words as OP, and has its most significant bit set.
                     48:
                     49:    2. Get a 9-bit approximation to sqrt(OP) using the pre-computed tables.
                     50:       This approximation is used for the first single-precision
                     51:       iterations of Newton's method, yielding a full-word approximation
                     52:       to sqrt(OP).
                     53:
                     54:    3. Perform multiple-precision Newton iteration until we have the
                     55:       exact result.  Only about half of the input operand is used in
                     56:       this calculation, as the square root is perfectly determinable
                     57:       from just the higher half of a number.  */
                     58: 
                     59: /* Define this macro for IEEE P854 machines with a fast sqrt instruction.  */
                     60: #if defined __GNUC__ && ! defined __SOFT_FLOAT
                     61:
                     62: #if defined __sparc__
                     63: #define SQRT(a) \
                     64:   ({                                                                   \
                     65:     double __sqrt_res;                                                 \
                     66:     asm ("fsqrtd %1,%0" : "=f" (__sqrt_res) : "f" (a));                        \
                     67:     __sqrt_res;                                                                \
                     68:   })
                     69: #endif
                     70:
                     71: #if defined __HAVE_68881__
                     72: #define SQRT(a) \
                     73:   ({                                                                   \
                     74:     double __sqrt_res;                                                 \
                     75:     asm ("fsqrtx %1,%0" : "=f" (__sqrt_res) : "f" (a));                        \
                     76:     __sqrt_res;                                                                \
                     77:   })
                     78: #endif
                     79:
                     80: #if defined __hppa
                     81: #define SQRT(a) \
                     82:   ({                                                                   \
                     83:     double __sqrt_res;                                                 \
                     84:     asm ("fsqrt,dbl %1,%0" : "=fx" (__sqrt_res) : "fx" (a));           \
                     85:     __sqrt_res;                                                                \
                     86:   })
                     87: #endif
                     88:
                     89: #if defined _ARCH_PWR2
                     90: #define SQRT(a) \
                     91:   ({                                                                   \
                     92:     double __sqrt_res;                                                 \
                     93:     asm ("fsqrt %0,%1" : "=f" (__sqrt_res) : "f" (a));                 \
                     94:     __sqrt_res;                                                                \
                     95:   })
                     96: #endif
                     97:
                     98: #endif
                     99:
                    100: #ifndef SQRT
                    101:
                    102: /* Tables for initial approximation of the square root.  These are
                    103:    indexed with bits 1-8 of the operand for which the square root is
                    104:    calculated, where bit 0 is the most significant non-zero bit.  I.e.
                    105:    the most significant one-bit is not used, since that per definition
                    106:    is one.  Likewise, the tables don't return the highest bit of the
                    107:    result.  That bit must be inserted by or:ing the returned value with
                    108:    0x100.  This way, we get a 9-bit approximation from 8-bit tables!  */
                    109:
                    110: /* Table to be used for operands with an even total number of bits.
                    111:    (Exactly as in the decimal system there are similarities between the
                    112:    square root of numbers with the same initial digits and an even
                    113:    difference in the total number of digits.  Consider the square root
                    114:    of 1, 10, 100, 1000, ...)  */
                    115: static unsigned char even_approx_tab[256] =
                    116: {
                    117:   0x6a, 0x6a, 0x6b, 0x6c, 0x6c, 0x6d, 0x6e, 0x6e,
                    118:   0x6f, 0x70, 0x71, 0x71, 0x72, 0x73, 0x73, 0x74,
                    119:   0x75, 0x75, 0x76, 0x77, 0x77, 0x78, 0x79, 0x79,
                    120:   0x7a, 0x7b, 0x7b, 0x7c, 0x7d, 0x7d, 0x7e, 0x7f,
                    121:   0x80, 0x80, 0x81, 0x81, 0x82, 0x83, 0x83, 0x84,
                    122:   0x85, 0x85, 0x86, 0x87, 0x87, 0x88, 0x89, 0x89,
                    123:   0x8a, 0x8b, 0x8b, 0x8c, 0x8d, 0x8d, 0x8e, 0x8f,
                    124:   0x8f, 0x90, 0x90, 0x91, 0x92, 0x92, 0x93, 0x94,
                    125:   0x94, 0x95, 0x96, 0x96, 0x97, 0x97, 0x98, 0x99,
                    126:   0x99, 0x9a, 0x9b, 0x9b, 0x9c, 0x9c, 0x9d, 0x9e,
                    127:   0x9e, 0x9f, 0xa0, 0xa0, 0xa1, 0xa1, 0xa2, 0xa3,
                    128:   0xa3, 0xa4, 0xa4, 0xa5, 0xa6, 0xa6, 0xa7, 0xa7,
                    129:   0xa8, 0xa9, 0xa9, 0xaa, 0xaa, 0xab, 0xac, 0xac,
                    130:   0xad, 0xad, 0xae, 0xaf, 0xaf, 0xb0, 0xb0, 0xb1,
                    131:   0xb2, 0xb2, 0xb3, 0xb3, 0xb4, 0xb5, 0xb5, 0xb6,
                    132:   0xb6, 0xb7, 0xb7, 0xb8, 0xb9, 0xb9, 0xba, 0xba,
                    133:   0xbb, 0xbb, 0xbc, 0xbd, 0xbd, 0xbe, 0xbe, 0xbf,
                    134:   0xc0, 0xc0, 0xc1, 0xc1, 0xc2, 0xc2, 0xc3, 0xc3,
                    135:   0xc4, 0xc5, 0xc5, 0xc6, 0xc6, 0xc7, 0xc7, 0xc8,
                    136:   0xc9, 0xc9, 0xca, 0xca, 0xcb, 0xcb, 0xcc, 0xcc,
                    137:   0xcd, 0xce, 0xce, 0xcf, 0xcf, 0xd0, 0xd0, 0xd1,
                    138:   0xd1, 0xd2, 0xd3, 0xd3, 0xd4, 0xd4, 0xd5, 0xd5,
                    139:   0xd6, 0xd6, 0xd7, 0xd7, 0xd8, 0xd9, 0xd9, 0xda,
                    140:   0xda, 0xdb, 0xdb, 0xdc, 0xdc, 0xdd, 0xdd, 0xde,
                    141:   0xde, 0xdf, 0xe0, 0xe0, 0xe1, 0xe1, 0xe2, 0xe2,
                    142:   0xe3, 0xe3, 0xe4, 0xe4, 0xe5, 0xe5, 0xe6, 0xe6,
                    143:   0xe7, 0xe7, 0xe8, 0xe8, 0xe9, 0xea, 0xea, 0xeb,
                    144:   0xeb, 0xec, 0xec, 0xed, 0xed, 0xee, 0xee, 0xef,
                    145:   0xef, 0xf0, 0xf0, 0xf1, 0xf1, 0xf2, 0xf2, 0xf3,
                    146:   0xf3, 0xf4, 0xf4, 0xf5, 0xf5, 0xf6, 0xf6, 0xf7,
                    147:   0xf7, 0xf8, 0xf8, 0xf9, 0xf9, 0xfa, 0xfa, 0xfb,
                    148:   0xfb, 0xfc, 0xfc, 0xfd, 0xfd, 0xfe, 0xfe, 0xff,
                    149: };
                    150:
                    151: /* Table to be used for operands with an odd total number of bits.
                    152:    (Further comments before previous table.)  */
                    153: static unsigned char odd_approx_tab[256] =
                    154: {
                    155:   0x00, 0x00, 0x00, 0x01, 0x01, 0x02, 0x02, 0x03,
                    156:   0x03, 0x04, 0x04, 0x05, 0x05, 0x06, 0x06, 0x07,
                    157:   0x07, 0x08, 0x08, 0x09, 0x09, 0x0a, 0x0a, 0x0b,
                    158:   0x0b, 0x0c, 0x0c, 0x0d, 0x0d, 0x0e, 0x0e, 0x0f,
                    159:   0x0f, 0x10, 0x10, 0x10, 0x11, 0x11, 0x12, 0x12,
                    160:   0x13, 0x13, 0x14, 0x14, 0x15, 0x15, 0x16, 0x16,
                    161:   0x16, 0x17, 0x17, 0x18, 0x18, 0x19, 0x19, 0x1a,
                    162:   0x1a, 0x1b, 0x1b, 0x1b, 0x1c, 0x1c, 0x1d, 0x1d,
                    163:   0x1e, 0x1e, 0x1f, 0x1f, 0x20, 0x20, 0x20, 0x21,
                    164:   0x21, 0x22, 0x22, 0x23, 0x23, 0x23, 0x24, 0x24,
                    165:   0x25, 0x25, 0x26, 0x26, 0x27, 0x27, 0x27, 0x28,
                    166:   0x28, 0x29, 0x29, 0x2a, 0x2a, 0x2a, 0x2b, 0x2b,
                    167:   0x2c, 0x2c, 0x2d, 0x2d, 0x2d, 0x2e, 0x2e, 0x2f,
                    168:   0x2f, 0x30, 0x30, 0x30, 0x31, 0x31, 0x32, 0x32,
                    169:   0x32, 0x33, 0x33, 0x34, 0x34, 0x35, 0x35, 0x35,
                    170:   0x36, 0x36, 0x37, 0x37, 0x37, 0x38, 0x38, 0x39,
                    171:   0x39, 0x39, 0x3a, 0x3a, 0x3b, 0x3b, 0x3b, 0x3c,
                    172:   0x3c, 0x3d, 0x3d, 0x3d, 0x3e, 0x3e, 0x3f, 0x3f,
                    173:   0x40, 0x40, 0x40, 0x41, 0x41, 0x41, 0x42, 0x42,
                    174:   0x43, 0x43, 0x43, 0x44, 0x44, 0x45, 0x45, 0x45,
                    175:   0x46, 0x46, 0x47, 0x47, 0x47, 0x48, 0x48, 0x49,
                    176:   0x49, 0x49, 0x4a, 0x4a, 0x4b, 0x4b, 0x4b, 0x4c,
                    177:   0x4c, 0x4c, 0x4d, 0x4d, 0x4e, 0x4e, 0x4e, 0x4f,
                    178:   0x4f, 0x50, 0x50, 0x50, 0x51, 0x51, 0x51, 0x52,
                    179:   0x52, 0x53, 0x53, 0x53, 0x54, 0x54, 0x54, 0x55,
                    180:   0x55, 0x56, 0x56, 0x56, 0x57, 0x57, 0x57, 0x58,
                    181:   0x58, 0x59, 0x59, 0x59, 0x5a, 0x5a, 0x5a, 0x5b,
                    182:   0x5b, 0x5b, 0x5c, 0x5c, 0x5d, 0x5d, 0x5d, 0x5e,
                    183:   0x5e, 0x5e, 0x5f, 0x5f, 0x60, 0x60, 0x60, 0x61,
                    184:   0x61, 0x61, 0x62, 0x62, 0x62, 0x63, 0x63, 0x63,
                    185:   0x64, 0x64, 0x65, 0x65, 0x65, 0x66, 0x66, 0x66,
                    186:   0x67, 0x67, 0x67, 0x68, 0x68, 0x68, 0x69, 0x69,
                    187: };
                    188: #endif
                    189:
                    190: 
                    191: mp_size_t
                    192: #if __STDC__
                    193: mpn_sqrtrem (mp_ptr root_ptr, mp_ptr rem_ptr, mp_srcptr op_ptr, mp_size_t op_size)
                    194: #else
                    195: mpn_sqrtrem (root_ptr, rem_ptr, op_ptr, op_size)
                    196:      mp_ptr root_ptr;
                    197:      mp_ptr rem_ptr;
                    198:      mp_srcptr op_ptr;
                    199:      mp_size_t op_size;
                    200: #endif
                    201: {
                    202:   /* R (root result) */
                    203:   mp_ptr rp;                   /* Pointer to least significant word */
                    204:   mp_size_t rsize;             /* The size in words */
                    205:
                    206:   /* T (OP shifted to the left a.k.a. normalized) */
                    207:   mp_ptr tp;                   /* Pointer to least significant word */
                    208:   mp_size_t tsize;             /* The size in words */
                    209:   mp_ptr t_end_ptr;            /* Pointer right beyond most sign. word */
                    210:   mp_limb_t t_high0, t_high1;  /* The two most significant words */
                    211:
                    212:   /* TT (temporary for numerator/remainder) */
                    213:   mp_ptr ttp;                  /* Pointer to least significant word */
                    214:
                    215:   /* X (temporary for quotient in main loop) */
                    216:   mp_ptr xp;                   /* Pointer to least significant word */
                    217:   mp_size_t xsize;             /* The size in words */
                    218:
                    219:   unsigned cnt;
                    220:   mp_limb_t initial_approx;    /* Initially made approximation */
                    221:   mp_size_t tsizes[BITS_PER_MP_LIMB];  /* Successive calculation precisions */
                    222:   mp_size_t tmp;
                    223:   mp_size_t i;
                    224:
                    225:   mp_limb_t cy_limb;
                    226:   TMP_DECL (marker);
                    227:
                    228:   /* If OP is zero, both results are zero.  */
                    229:   if (op_size == 0)
                    230:     return 0;
                    231:
                    232:   count_leading_zeros (cnt, op_ptr[op_size - 1]);
                    233:   tsize = op_size;
                    234:   if ((tsize & 1) != 0)
                    235:     {
                    236:       cnt += BITS_PER_MP_LIMB;
                    237:       tsize++;
                    238:     }
                    239:
                    240:   rsize = tsize / 2;
                    241:   rp = root_ptr;
                    242:
                    243:   TMP_MARK (marker);
                    244:
                    245:   /* Shift OP an even number of bits into T, such that either the most or
                    246:      the second most significant bit is set, and such that the number of
                    247:      words in T becomes even.  This way, the number of words in R=sqrt(OP)
                    248:      is exactly half as many as in OP, and the most significant bit of R
                    249:      is set.
                    250:
                    251:      Also, the initial approximation is simplified by this up-shifted OP.
                    252:
                    253:      Finally, the Newtonian iteration which is the main part of this
                    254:      program performs division by R.  The fast division routine expects
                    255:      the divisor to be "normalized" in exactly the sense of having the
                    256:      most significant bit set.  */
                    257:
                    258:   tp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB);
                    259:
                    260:   if ((cnt & ~1) % BITS_PER_MP_LIMB != 0)
                    261:     t_high0 = mpn_lshift (tp + cnt / BITS_PER_MP_LIMB, op_ptr, op_size,
                    262:                          (cnt & ~1) % BITS_PER_MP_LIMB);
                    263:   else
                    264:     MPN_COPY (tp + cnt / BITS_PER_MP_LIMB, op_ptr, op_size);
                    265:
                    266:   if (cnt >= BITS_PER_MP_LIMB)
                    267:     tp[0] = 0;
                    268:
                    269:   t_high0 = tp[tsize - 1];
                    270:   t_high1 = tp[tsize - 2];     /* Never stray.  TSIZE is >= 2.  */
                    271:
                    272: /* Is there a fast sqrt instruction defined for this machine?  */
                    273: #ifdef SQRT
                    274:   {
                    275:     initial_approx = SQRT (t_high0 * 2.0
                    276:                           * ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1))
                    277:                           + t_high1);
                    278:     /* If t_high0,,t_high1 is big, the result in INITIAL_APPROX might have
                    279:        become incorrect due to overflow in the conversion from double to
                    280:        mp_limb_t above.  It will typically be zero in that case, but might be
                    281:        a small number on some machines.  The most significant bit of
                    282:        INITIAL_APPROX should be set, so that bit is a good overflow
                    283:        indication.  */
                    284:     if ((mp_limb_signed_t) initial_approx >= 0)
                    285:       initial_approx = ~(mp_limb_t)0;
                    286:   }
                    287: #else
                    288:   /* Get a 9 bit approximation from the tables.  The tables expect to
                    289:      be indexed with the 8 high bits right below the highest bit.
                    290:      Also, the highest result bit is not returned by the tables, and
                    291:      must be or:ed into the result.  The scheme gives 9 bits of start
                    292:      approximation with just 256-entry 8 bit tables.  */
                    293:
                    294:   if ((cnt & 1) == 0)
                    295:     {
                    296:       /* The most sign bit of t_high0 is set.  */
                    297:       initial_approx = t_high0 >> (BITS_PER_MP_LIMB - 8 - 1);
                    298:       initial_approx &= 0xff;
                    299:       initial_approx = even_approx_tab[initial_approx];
                    300:     }
                    301:   else
                    302:     {
                    303:       /* The most significant bit of T_HIGH0 is unset,
                    304:         the second most significant is set.  */
                    305:       initial_approx = t_high0 >> (BITS_PER_MP_LIMB - 8 - 2);
                    306:       initial_approx &= 0xff;
                    307:       initial_approx = odd_approx_tab[initial_approx];
                    308:     }
                    309:   initial_approx |= 0x100;
                    310:   initial_approx <<= BITS_PER_MP_LIMB - 8 - 1;
                    311:
                    312:   /* Perform small precision Newtonian iterations to get a full word
                    313:      approximation.  For small operands, these iteration will make the
                    314:      entire job.  */
                    315:   if (t_high0 == ~(mp_limb_t)0)
                    316:     initial_approx = t_high0;
                    317:   else
                    318:     {
                    319:       mp_limb_t quot;
                    320:
                    321:       if (t_high0 >= initial_approx)
                    322:        initial_approx = t_high0 + 1;
                    323:
                    324:       /* First get about 18 bits with pure C arithmetics.  */
                    325:       quot = t_high0 / (initial_approx >> BITS_PER_MP_LIMB/2) << BITS_PER_MP_LIMB/2;
                    326:       initial_approx = (initial_approx + quot) / 2;
                    327:       initial_approx |= (mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1);
                    328:
                    329:       /* Now get a full word by one (or for > 36 bit machines) several
                    330:         iterations.  */
                    331:       for (i = 16; i < BITS_PER_MP_LIMB; i <<= 1)
                    332:        {
                    333:          mp_limb_t ignored_remainder;
                    334:
                    335:          udiv_qrnnd (quot, ignored_remainder,
                    336:                      t_high0, t_high1, initial_approx);
                    337:          initial_approx = (initial_approx + quot) / 2;
                    338:          initial_approx |= (mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1);
                    339:        }
                    340:     }
                    341: #endif
                    342:
                    343:   rp[0] = initial_approx;
                    344:   rsize = 1;
                    345:
                    346: #ifdef DEBUG
                    347:          printf ("\n\nT = ");
                    348:          mpn_dump (tp, tsize);
                    349: #endif
                    350:
                    351:   if (tsize > 2)
                    352:     {
                    353:       /* Determine the successive precisions to use in the iteration.  We
                    354:         minimize the precisions, beginning with the highest (i.e. last
                    355:         iteration) to the lowest (i.e. first iteration).  */
                    356:
                    357:       xp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB);
                    358:       ttp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB);
                    359:
                    360:       t_end_ptr = tp + tsize;
                    361:
                    362:       tmp = tsize / 2;
                    363:       for (i = 0;; i++)
                    364:        {
                    365:          tsize = (tmp + 1) / 2;
                    366:          if (tmp == tsize)
                    367:            break;
                    368:          tsizes[i] = tsize + tmp;
                    369:          tmp = tsize;
                    370:        }
                    371:
                    372:       /* Main Newton iteration loop.  For big arguments, most of the
                    373:         time is spent here.  */
                    374:
                    375:       /* It is possible to do a great optimization here.  The successive
                    376:         divisors in the mpn_divmod call below has more and more leading
                    377:         words equal to its predecessor.  Therefore the beginning of
                    378:         each division will repeat the same work as did the last
                    379:         division.  If we could guarantee that the leading words of two
                    380:         consecutive divisors are the same (i.e. in this case, a later
                    381:         divisor has just more digits at the end) it would be a simple
                    382:         matter of just using the old remainder of the last division in
                    383:         a subsequent division, to take care of this optimization.  This
                    384:         idea would surely make a difference even for small arguments.  */
                    385:
                    386:       /* Loop invariants:
                    387:
                    388:         R <= shiftdown_to_same_size(floor(sqrt(OP))) < R + 1.
                    389:         X - 1 < shiftdown_to_same_size(floor(sqrt(OP))) <= X.
                    390:         R <= shiftdown_to_same_size(X).  */
                    391:
                    392:       while (--i >= 0)
                    393:        {
                    394:          mp_limb_t cy;
                    395: #ifdef DEBUG
                    396:          mp_limb_t old_least_sign_r = rp[0];
                    397:          mp_size_t old_rsize = rsize;
                    398:
                    399:          printf ("R = ");
                    400:          mpn_dump (rp, rsize);
                    401: #endif
                    402:          tsize = tsizes[i];
                    403:
                    404:          /* Need to copy the numerator into temporary space, as
                    405:             mpn_divmod overwrites its numerator argument with the
                    406:             remainder (which we currently ignore).  */
                    407:          MPN_COPY (ttp, t_end_ptr - tsize, tsize);
                    408:          cy = mpn_divmod (xp, ttp, tsize, rp, rsize);
                    409:          xsize = tsize - rsize;
                    410:
                    411: #ifdef DEBUG
                    412:          printf ("X =%d ", cy);
                    413:          mpn_dump (xp, xsize);
                    414: #endif
                    415:
                    416:          /* Add X and R with the most significant limbs aligned,
                    417:             temporarily ignoring at least one limb at the low end of X.  */
                    418:          tmp = xsize - rsize;
                    419:          cy += mpn_add_n (xp + tmp, rp, xp + tmp, rsize);
                    420:
                    421:          /* If T begins with more than 2 x BITS_PER_MP_LIMB of ones, we get
                    422:             intermediate roots that'd need an extra bit.  We don't want to
                    423:             handle that since it would make the subsequent divisor
                    424:             non-normalized, so round such roots down to be only ones in the
                    425:             current precision.  */
                    426:          if (cy == 2)
                    427:            {
                    428:              mp_size_t j;
                    429:              for (j = xsize; j >= 0; j--)
                    430:                xp[j] = ~(mp_limb_t)0;
                    431:            }
                    432:
                    433:          /* Divide X by 2 and put the result in R.  This is the new
                    434:             approximation.  Shift in the carry from the addition.  */
                    435:          mpn_rshift (rp, xp, xsize, 1);
                    436:          rp[xsize - 1] |= ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1));
                    437:          rsize = xsize;
                    438: #ifdef DEBUG
                    439:          if (old_least_sign_r != rp[rsize - old_rsize])
                    440:            printf (">>>>>>>> %d: %0*lX, %0*lX <<<<<<<<\n",
                    441:                    i, 2 * BYTES_PER_MP_LIMB, old_least_sign_r,
                    442:                    2 * BYTES_PER_MP_LIMB, rp[rsize - old_rsize]);
                    443: #endif
                    444:        }
                    445:     }
                    446:
                    447: #ifdef DEBUG
                    448:   printf ("(final) R = ");
                    449:   mpn_dump (rp, rsize);
                    450: #endif
                    451:
                    452:   /* We computed the square root of OP * 2**(2*floor(cnt/2)).
                    453:      This has resulted in R being 2**floor(cnt/2) to large.
                    454:      Shift it down here to fix that.  */
                    455:   if (cnt / 2 != 0)
                    456:     {
                    457:       mpn_rshift (rp, rp, rsize, cnt/2);
                    458:       rsize -= rp[rsize - 1] == 0;
                    459:     }
                    460:
                    461:   /* Calculate the remainder.  */
                    462:   mpn_mul_n (tp, rp, rp, rsize);
                    463:   tsize = rsize + rsize;
                    464:   tsize -= tp[tsize - 1] == 0;
                    465:   if (op_size < tsize
                    466:       || (op_size == tsize && mpn_cmp (op_ptr, tp, op_size) < 0))
                    467:     {
                    468:       /* R is too large.  Decrement it.  */
                    469:
                    470:       /* These operations can't overflow.  */
                    471:       cy_limb  = mpn_sub_n (tp, tp, rp, rsize);
                    472:       cy_limb += mpn_sub_n (tp, tp, rp, rsize);
                    473:       mpn_sub_1 (tp + rsize, tp + rsize, tsize - rsize, cy_limb);
                    474:       mpn_add_1 (tp, tp, tsize, (mp_limb_t) 1);
                    475:
                    476:       mpn_sub_1 (rp, rp, rsize, (mp_limb_t) 1);
                    477:
                    478: #ifdef DEBUG
                    479:       printf ("(adjusted) R = ");
                    480:       mpn_dump (rp, rsize);
                    481: #endif
                    482:     }
                    483:
                    484:   if (rem_ptr != NULL)
                    485:     {
                    486:       cy_limb = mpn_sub (rem_ptr, op_ptr, op_size, tp, tsize);
                    487:       MPN_NORMALIZE (rem_ptr, op_size);
                    488:       TMP_FREE (marker);
                    489:       return op_size;
                    490:     }
                    491:   else
                    492:     {
                    493:       int res;
                    494:       res = op_size != tsize || mpn_cmp (op_ptr, tp, op_size);
                    495:       TMP_FREE (marker);
                    496:       return res;
                    497:     }
                    498: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>