← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/ecl/intrepid

« back to all changes in this revision

Viewing changes to src/gmp/tests/refmpn.c

  • Committer: Bazaar Package Importer
  • Author(s): Peter Van Eynde
  • Date: 2007-04-09 11:51:51 UTC
  • mfrom: (1.1.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20070409115151-ql8cr0kalzx1jmla
Tags: 0.9i-20070324-2
Upload to unstable. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/gmp/tests/refmpn.c
1
1
/* Reference mpn functions, designed to be simple, portable and independent
2
2
   of the normal gmp code.  Speed isn't a consideration.
3
3
 
4
 
Copyright 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2004 Free Software
5
 
Foundation, Inc.
 
4
Copyright 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 Free
 
5
Software Foundation, Inc.
6
6
 
7
7
This file is part of the GNU MP Library.
8
8
 
18
18
 
19
19
You should have received a copy of the GNU Lesser General Public License
20
20
along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
21
 
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
22
 
MA 02111-1307, USA. */
 
21
the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 
22
MA 02110-1301, USA. */
23
23
 
24
24
 
25
25
/* Most routines have assertions representing what the mpn routines are
111
111
  ASSERT (size >= 0);
112
112
  if (size == 0)
113
113
    size = 1;
114
 
  p = (mp_ptr) malloc (size * BYTES_PER_MP_LIMB);
 
114
  p = (mp_ptr) malloc ((size_t) (size * BYTES_PER_MP_LIMB));
115
115
  ASSERT (p != NULL);
116
116
  return p;
117
117
}
118
118
 
 
119
/* Free limbs allocated by refmpn_malloc_limbs. NOTE: Can't free
 
120
 * memory allocated by refmpn_malloc_limbs_aligned. */
 
121
void
 
122
refmpn_free_limbs (mp_ptr p)
 
123
{
 
124
  free (p);
 
125
}
 
126
 
119
127
mp_ptr
120
128
refmpn_memdup_limbs (mp_srcptr ptr, mp_size_t size)
121
129
{
127
135
 
128
136
/* malloc n limbs on a multiple of m bytes boundary */
129
137
mp_ptr
130
 
refmpn_malloc_limbs_aligned (size_t n, size_t m)
 
138
refmpn_malloc_limbs_aligned (mp_size_t n, size_t m)
131
139
{
132
140
  return (mp_ptr) align_pointer (refmpn_malloc_limbs (n + m-1), m);
133
141
}
148
156
  refmpn_fill (ptr, size, CNST_LIMB(0));
149
157
}
150
158
 
 
159
void
 
160
refmpn_zero_extend (mp_ptr ptr, mp_size_t oldsize, mp_size_t newsize)
 
161
{
 
162
  ASSERT (newsize >= oldsize);
 
163
  refmpn_zero (ptr+oldsize, newsize-oldsize);
 
164
}
 
165
 
151
166
int
152
167
refmpn_zero_p (mp_srcptr ptr, mp_size_t size)
153
168
{
158
173
  return 1;
159
174
}
160
175
 
 
176
mp_size_t
 
177
refmpn_normalize (mp_srcptr ptr, mp_size_t size)
 
178
{
 
179
  ASSERT (size >= 0);
 
180
  while (size > 0 && ptr[size-1] == 0)
 
181
    size--;
 
182
  return size;
 
183
}
 
184
 
161
185
/* the highest one bit in x */
162
186
mp_limb_t
163
187
refmpn_msbone (mp_limb_t x)
183
207
  return (refmpn_msbone (x) << 1) - 1;
184
208
}
185
209
 
 
210
/* How many digits in the given base will fit in a limb.
 
211
   Notice that the product b is allowed to be equal to the limit
 
212
   2^GMP_NUMB_BITS, this ensures the result for base==2 will be
 
213
   GMP_NUMB_BITS (and similarly other powers of 2).  */
 
214
int
 
215
refmpn_chars_per_limb (int base)
 
216
{
 
217
  mp_limb_t  limit[2], b[2];
 
218
  int        chars_per_limb;
 
219
 
 
220
  ASSERT (base >= 2);
 
221
 
 
222
  limit[0] = 0;  /* limit = 2^GMP_NUMB_BITS */
 
223
  limit[1] = 1;
 
224
  b[0] = 1;      /* b = 1 */
 
225
  b[1] = 0;
 
226
 
 
227
  chars_per_limb = 0;
 
228
  for (;;)
 
229
    {
 
230
      if (refmpn_mul_1 (b, b, (mp_size_t) 2, (mp_limb_t) base))
 
231
        break;
 
232
      if (refmpn_cmp (b, limit, (mp_size_t) 2) > 0)
 
233
        break;
 
234
      chars_per_limb++;
 
235
    }
 
236
  return chars_per_limb;
 
237
}
 
238
 
 
239
/* The biggest value base**n which fits in GMP_NUMB_BITS. */
 
240
mp_limb_t
 
241
refmpn_big_base (int base)
 
242
{
 
243
  int        chars_per_limb = refmpn_chars_per_limb (base);
 
244
  int        i;
 
245
  mp_limb_t  bb;
 
246
 
 
247
  ASSERT (base >= 2);
 
248
  bb = 1;
 
249
  for (i = 0; i < chars_per_limb; i++)
 
250
    bb *= base;
 
251
  return bb;
 
252
}
 
253
 
186
254
 
187
255
void
188
256
refmpn_setbit (mp_ptr ptr, unsigned long bit)
252
320
    rp[i] = sp[i];
253
321
}
254
322
 
 
323
/* Copy {xp,xsize} to {wp,wsize}.  If x is shorter, then pad w with low
 
324
   zeros to wsize.  If x is longer, then copy just the high wsize limbs.  */
 
325
void
 
326
refmpn_copy_extend (mp_ptr wp, mp_size_t wsize, mp_srcptr xp, mp_size_t xsize)
 
327
{
 
328
  ASSERT (wsize >= 0);
 
329
  ASSERT (xsize >= 0);
 
330
 
 
331
  /* high part of x if x bigger than w */
 
332
  if (xsize > wsize)
 
333
    {
 
334
      xp += xsize - wsize;
 
335
      xsize = wsize;
 
336
    }
 
337
 
 
338
  refmpn_copy (wp + wsize-xsize, xp, xsize);
 
339
  refmpn_zero (wp, wsize-xsize);
 
340
}
 
341
 
255
342
void
256
343
refmpn_com_n (mp_ptr rp, mp_srcptr sp, mp_size_t size)
257
344
{
380
467
  LOGOPS ((s1p[i] ^ s2p[i]) ^ GMP_NUMB_MASK);
381
468
}
382
469
 
 
470
 
 
471
/* set *dh,*dl to mh:ml - sh:sl, in full limbs */
 
472
void
 
473
refmpn_sub_ddmmss (mp_limb_t *dh, mp_limb_t *dl,
 
474
                   mp_limb_t mh, mp_limb_t ml, mp_limb_t sh, mp_limb_t sl)
 
475
{
 
476
  *dl = ml - sl;
 
477
  *dh = mh - sh - (ml < sl);
 
478
}
 
479
 
 
480
 
383
481
/* set *w to x+y, return 0 or 1 carry */
384
482
mp_limb_t
385
 
add (mp_limb_t *w, mp_limb_t x, mp_limb_t y)
 
483
ref_addc_limb (mp_limb_t *w, mp_limb_t x, mp_limb_t y)
386
484
{
387
485
  mp_limb_t  sum, cy;
388
486
 
402
500
 
403
501
/* set *w to x-y, return 0 or 1 borrow */
404
502
mp_limb_t
405
 
sub (mp_limb_t *w, mp_limb_t x, mp_limb_t y)
 
503
ref_subc_limb (mp_limb_t *w, mp_limb_t x, mp_limb_t y)
406
504
{
407
505
  mp_limb_t  diff, cy;
408
506
 
430
528
  ASSERT_LIMB (y);
431
529
  ASSERT (c == 0 || c == 1);
432
530
 
433
 
  r = add (w, x, y);
434
 
  return r + add (w, *w, c);
 
531
  r = ref_addc_limb (w, x, y);
 
532
  return r + ref_addc_limb (w, *w, c);
435
533
}
436
534
 
437
535
/* set *w to x-y-c (where c is 0 or 1), return 0 or 1 borrow */
444
542
  ASSERT_LIMB (y);
445
543
  ASSERT (c == 0 || c == 1);
446
544
 
447
 
  r = sub (w, x, y);
448
 
  return r + sub (w, *w, c);
 
545
  r = ref_subc_limb (w, x, y);
 
546
  return r + ref_subc_limb (w, *w, c);
449
547
}
450
548
 
451
549
 
466
564
mp_limb_t
467
565
refmpn_add_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t n)
468
566
{
469
 
  AORS_1 (add);
 
567
  AORS_1 (ref_addc_limb);
470
568
}
471
569
mp_limb_t
472
570
refmpn_sub_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t n)
473
571
{
474
 
  AORS_1 (sub);
 
572
  AORS_1 (ref_subc_limb);
475
573
}
476
574
 
477
575
#define AORS_NC(operation)                                              \
514
612
  return refmpn_sub_nc (rp, s1p, s2p, size, CNST_LIMB(0));
515
613
}
516
614
 
 
615
mp_limb_t
 
616
refmpn_addlsh1_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
 
617
{
 
618
  mp_limb_t cy;
 
619
  mp_ptr tp;
 
620
 
 
621
  ASSERT (refmpn_overlap_fullonly_two_p (rp, up, vp, n));
 
622
  ASSERT (n >= 1);
 
623
  ASSERT_MPN (up, n);
 
624
  ASSERT_MPN (vp, n);
 
625
 
 
626
  tp = refmpn_malloc_limbs (n);
 
627
  cy  = refmpn_lshift (tp, vp, n, 1);
 
628
  cy += refmpn_add_n (rp, up, tp, n);
 
629
  free (tp);
 
630
  return cy;
 
631
}
 
632
mp_limb_t
 
633
refmpn_sublsh1_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
 
634
{
 
635
  mp_limb_t cy;
 
636
  mp_ptr tp;
 
637
 
 
638
  ASSERT (refmpn_overlap_fullonly_two_p (rp, up, vp, n));
 
639
  ASSERT (n >= 1);
 
640
  ASSERT_MPN (up, n);
 
641
  ASSERT_MPN (vp, n);
 
642
 
 
643
  tp = refmpn_malloc_limbs (n);
 
644
  cy  = mpn_lshift (tp, vp, n, 1);
 
645
  cy += mpn_sub_n (rp, up, tp, n);
 
646
  free (tp);
 
647
  return cy;
 
648
}
 
649
mp_limb_t
 
650
refmpn_rsh1add_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
 
651
{
 
652
  mp_limb_t cya, cys;
 
653
 
 
654
  ASSERT (refmpn_overlap_fullonly_two_p (rp, up, vp, n));
 
655
  ASSERT (n >= 1);
 
656
  ASSERT_MPN (up, n);
 
657
  ASSERT_MPN (vp, n);
 
658
 
 
659
  cya = mpn_add_n (rp, up, vp, n);
 
660
  cys = mpn_rshift (rp, rp, n, 1) >> (GMP_NUMB_BITS - 1);
 
661
  rp[n - 1] |= cya << (GMP_NUMB_BITS - 1);
 
662
  return cys;
 
663
}
 
664
mp_limb_t
 
665
refmpn_rsh1sub_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
 
666
{
 
667
  mp_limb_t cya, cys;
 
668
 
 
669
  ASSERT (refmpn_overlap_fullonly_two_p (rp, up, vp, n));
 
670
  ASSERT (n >= 1);
 
671
  ASSERT_MPN (up, n);
 
672
  ASSERT_MPN (vp, n);
 
673
 
 
674
  cya = mpn_sub_n (rp, up, vp, n);
 
675
  cys = mpn_rshift (rp, rp, n, 1) >> (GMP_NUMB_BITS - 1);
 
676
  rp[n - 1] |= cya << (GMP_NUMB_BITS - 1);
 
677
  return cys;
 
678
}
517
679
 
518
680
/* Twos complement, return borrow. */
519
681
mp_limb_t
567
729
#define LOWPART(x)   ((x) & LOWMASK)
568
730
#define HIGHPART(x)  SHIFTLOW((x) & HIGHMASK)
569
731
 
570
 
/* Set *hi,*lo to x*y, using full limbs not nails. */
 
732
/* Set return:*lo to x*y, using full limbs not nails. */
571
733
mp_limb_t
572
734
refmpn_umul_ppmm (mp_limb_t *lo, mp_limb_t x, mp_limb_t y)
573
735
{
615
777
    {
616
778
      hi = refmpn_umul_ppmm (&lo, sp[i], multiplier);
617
779
      lo >>= GMP_NAIL_BITS;
618
 
      ASSERT_NOCARRY (add (&hi, hi, add (&lo, lo, carry)));
 
780
      ASSERT_NOCARRY (ref_addc_limb (&hi, hi, ref_addc_limb (&lo, lo, carry)));
619
781
      rp[i] = lo;
620
782
      carry = hi;
621
783
    }
636
798
  mp_limb_t  c;
637
799
 
638
800
  ASSERT (refmpn_overlap_fullonly_p (dst, src, size));
639
 
  ASSERT (! refmpn_overlap_p (dst, size+1, mult, 2));
 
801
  ASSERT (! refmpn_overlap_p (dst, size+1, mult, (mp_size_t) 2));
640
802
  ASSERT (size >= 1);
641
803
  ASSERT_MPN (mult, 2);
642
804
 
694
856
 
695
857
 
696
858
mp_limb_t
 
859
refmpn_addmul_N (mp_ptr dst, mp_srcptr src, mp_size_t size,
 
860
                 mp_srcptr mult, mp_size_t msize)
 
861
{
 
862
  mp_ptr     src_copy;
 
863
  mp_limb_t  ret;
 
864
  mp_size_t  i;
 
865
 
 
866
  ASSERT (dst == src || ! refmpn_overlap_p (dst, size+msize-1, src, size));
 
867
  ASSERT (! refmpn_overlap_p (dst, size+msize-1, mult, msize));
 
868
  ASSERT (size >= msize);
 
869
  ASSERT_MPN (mult, msize);
 
870
 
 
871
  /* in case dst==src */
 
872
  src_copy = refmpn_malloc_limbs (size);
 
873
  refmpn_copyi (src_copy, src, size);
 
874
  src = src_copy;
 
875
 
 
876
  for (i = 0; i < msize-1; i++)
 
877
    dst[size+i] = refmpn_addmul_1 (dst+i, src, size, mult[i]);
 
878
  ret = refmpn_addmul_1 (dst+i, src, size, mult[i]);
 
879
 
 
880
  free (src_copy);
 
881
  return ret;
 
882
}
 
883
 
 
884
mp_limb_t
 
885
refmpn_addmul_2 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
 
886
{
 
887
  return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 2);
 
888
}
 
889
mp_limb_t
 
890
refmpn_addmul_3 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
 
891
{
 
892
  return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 3);
 
893
}
 
894
mp_limb_t
 
895
refmpn_addmul_4 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
 
896
{
 
897
  return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 4);
 
898
}
 
899
mp_limb_t
 
900
refmpn_addmul_5 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
 
901
{
 
902
  return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 5);
 
903
}
 
904
mp_limb_t
 
905
refmpn_addmul_6 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
 
906
{
 
907
  return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 6);
 
908
}
 
909
mp_limb_t
 
910
refmpn_addmul_7 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
 
911
{
 
912
  return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 7);
 
913
}
 
914
mp_limb_t
 
915
refmpn_addmul_8 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
 
916
{
 
917
  return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 8);
 
918
}
 
919
 
 
920
mp_limb_t
697
921
refmpn_addsub_nc (mp_ptr r1p, mp_ptr r2p,
698
922
                  mp_srcptr s1p, mp_srcptr s2p, mp_size_t size,
699
923
                  mp_limb_t carry)
792
1016
}
793
1017
 
794
1018
 
 
1019
/* accepting shift==0 and doing a plain copyi or copyd in that case */
795
1020
mp_limb_t
796
1021
refmpn_rshift_or_copy (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
797
1022
{
805
1030
      return refmpn_rshift (rp, sp, size, shift);
806
1031
    }
807
1032
}
808
 
 
809
1033
mp_limb_t
810
1034
refmpn_lshift_or_copy (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
811
1035
{
820
1044
    }
821
1045
}
822
1046
 
 
1047
/* accepting size==0 too */
 
1048
mp_limb_t
 
1049
refmpn_rshift_or_copy_any (mp_ptr rp, mp_srcptr sp, mp_size_t size,
 
1050
                           unsigned shift)
 
1051
{
 
1052
  return (size == 0 ? 0 : refmpn_rshift_or_copy (rp, sp, size, shift));
 
1053
}
 
1054
mp_limb_t
 
1055
refmpn_lshift_or_copy_any (mp_ptr rp, mp_srcptr sp, mp_size_t size,
 
1056
                           unsigned shift)
 
1057
{
 
1058
  return (size == 0 ? 0 : refmpn_lshift_or_copy (rp, sp, size, shift));
 
1059
}
823
1060
 
824
1061
/* Divide h,l by d, return quotient, store remainder to *rp.
825
1062
   Operates on full limbs, not nails.
868
1105
                             mp_limb_t divisor, mp_limb_t carry)
869
1106
{
870
1107
  mp_size_t  i;
 
1108
  mp_limb_t rem[1];
871
1109
  for (i = size-1; i >= 0; i--)
872
1110
    {
873
 
      rp[i] = refmpn_udiv_qrnnd (&carry, carry,
 
1111
      rp[i] = refmpn_udiv_qrnnd (rem, carry,
874
1112
                                 sp[i] << GMP_NAIL_BITS,
875
1113
                                 divisor << GMP_NAIL_BITS);
876
 
      carry >>= GMP_NAIL_BITS;
 
1114
      carry = *rem >> GMP_NAIL_BITS;
877
1115
    }
878
1116
  return carry;
879
1117
}
1007
1245
{
1008
1246
  mp_limb_t r;
1009
1247
  ASSERT (d & GMP_LIMB_HIGHBIT);
1010
 
  return refmpn_udiv_qrnnd (&r, -d-1, -1, d);
 
1248
  return refmpn_udiv_qrnnd (&r, -d-1, MP_LIMB_T_MAX, d);
1011
1249
}
1012
1250
 
1013
1251
 
1046
1284
  cs = refmpn_sub_1 (spcopy, sp, size, carry);
1047
1285
 
1048
1286
  for (c = 0; c <= 2; c++)
1049
 
    if (refmpn_divmod_1c (rp, spcopy, size, 3, c) == 0)
 
1287
    if (refmpn_divmod_1c (rp, spcopy, size, CNST_LIMB(3), c) == 0)
1050
1288
      goto done;
1051
1289
  ASSERT_FAIL (no value of c satisfies);
1052
1290
 
1268
1506
  return xsize;
1269
1507
}
1270
1508
 
 
1509
unsigned long
 
1510
ref_popc_limb (mp_limb_t src)
 
1511
{
 
1512
  unsigned long  count;
 
1513
  int  i;
 
1514
 
 
1515
  count = 0;
 
1516
  for (i = 0; i < GMP_LIMB_BITS; i++)
 
1517
    {
 
1518
      count += (src & 1);
 
1519
      src >>= 1;
 
1520
    }
 
1521
  return count;
 
1522
}
1271
1523
 
1272
1524
unsigned long
1273
1525
refmpn_popcount (mp_srcptr sp, mp_size_t size)
1274
1526
{
1275
1527
  unsigned long  count = 0;
1276
1528
  mp_size_t  i;
1277
 
  int        j;
1278
 
  mp_limb_t  l;
1279
1529
 
1280
1530
  ASSERT (size >= 0);
1281
1531
  ASSERT_MPN (sp, size);
1282
1532
 
1283
1533
  for (i = 0; i < size; i++)
1284
 
    {
1285
 
      l = sp[i];
1286
 
      for (j = 0; j < GMP_NUMB_BITS; j++)
1287
 
        {
1288
 
          count += (l & 1);
1289
 
          l >>= 1;
1290
 
        }
1291
 
    }
 
1534
    count += ref_popc_limb (sp[i]);
1292
1535
  return count;
1293
1536
}
1294
1537
 
1320
1563
  mp_limb_t  D[2];
1321
1564
  int        n;
1322
1565
 
1323
 
  ASSERT (! refmpn_overlap_p (r, 2, d, 2));
 
1566
  ASSERT (! refmpn_overlap_p (r, (mp_size_t) 2, d, (mp_size_t) 2));
1324
1567
  ASSERT_MPN (a, 2);
1325
1568
  ASSERT_MPN (d, 2);
1326
1569
 
1332
1575
    {
1333
1576
      if (D[1] & GMP_NUMB_HIGHBIT)
1334
1577
        break;
1335
 
      if (refmpn_cmp (r, D, 2) <= 0)
 
1578
      if (refmpn_cmp (r, D, (mp_size_t) 2) <= 0)
1336
1579
        break;
1337
 
      refmpn_lshift (D, D, 2, 1);
 
1580
      refmpn_lshift (D, D, (mp_size_t) 2, 1);
1338
1581
      n++;
1339
1582
      ASSERT (n <= GMP_NUMB_BITS);
1340
1583
    }
1341
1584
 
1342
1585
  while (n >= 0)
1343
1586
    {
1344
 
      if (refmpn_cmp (r, D, 2) >= 0)
1345
 
        ASSERT_NOCARRY (refmpn_sub_n (r, r, D, 2));
1346
 
      refmpn_rshift (D, D, 2, 1);
 
1587
      if (refmpn_cmp (r, D, (mp_size_t) 2) >= 0)
 
1588
        ASSERT_NOCARRY (refmpn_sub_n (r, r, D, (mp_size_t) 2));
 
1589
      refmpn_rshift (D, D, (mp_size_t) 2, 1);
1347
1590
      n--;
1348
1591
    }
1349
1592
 
1350
 
  ASSERT (refmpn_cmp (r, d, 2) < 0);
 
1593
  ASSERT (refmpn_cmp (r, d, (mp_size_t) 2) < 0);
1351
1594
}
1352
1595
 
1353
1596
 
1354
 
/* Find n where 0<n<2^GMP_NUMB_BITS and n==c mod m */
 
1597
/* Return n with 0 < n < 2^GMP_NUMB_BITS such that there exists 0 < |d| <
 
1598
   2^GMP_NUMB_BITS, and n == d * c mod 2^(2*GMP_NUMB_BITS).  */
1355
1599
mp_limb_t
1356
1600
refmpn_gcd_finda (const mp_limb_t c[2])
1357
1601
{
1394
1638
  ASSERT (nsize >= dsize);
1395
1639
  /* ASSERT (dsize > 2); */
1396
1640
  ASSERT (dsize >= 2);
1397
 
  ASSERT (dp[dsize-1] & GMP_LIMB_HIGHBIT);
 
1641
  ASSERT (dp[dsize-1] & GMP_NUMB_HIGHBIT);
1398
1642
  ASSERT (! refmpn_overlap_p (qp, nsize-dsize, np, nsize) || qp+dsize >= np);
1399
1643
  ASSERT_MPN (np, nsize);
1400
1644
  ASSERT_MPN (dp, dsize);
1414
1658
 
1415
1659
      ASSERT (n0 <= d1);
1416
1660
      if (n0 == d1)
1417
 
        q = MP_LIMB_T_MAX;
 
1661
        q = GMP_NUMB_MAX;
1418
1662
      else
1419
 
        q = refmpn_udiv_qrnnd (&dummy_r, n0, n1, d1);
 
1663
        q = refmpn_udiv_qrnnd (&dummy_r, n0, n1 << GMP_NAIL_BITS,
 
1664
                               d1 << GMP_NAIL_BITS);
1420
1665
 
1421
1666
      n0 -= refmpn_submul_1 (np+i, dp, dsize, q);
1422
1667
      ASSERT (n0 == 0 || n0 == MP_LIMB_T_MAX);
1464
1709
  ASSERT (qxn == 0);
1465
1710
  ASSERT_MPN (np, nsize);
1466
1711
  ASSERT_MPN (dp, dsize);
 
1712
  ASSERT (dsize > 0);
 
1713
  ASSERT (dp[dsize-1] != 0);
1467
1714
 
1468
1715
  if (dsize == 1)
1469
1716
    {
1474
1721
    {
1475
1722
      mp_ptr  n2p = refmpn_malloc_limbs (nsize+1);
1476
1723
      mp_ptr  d2p = refmpn_malloc_limbs (dsize);
1477
 
      int     norm = refmpn_count_leading_zeros (dp[dsize-1]);
 
1724
      int     norm = refmpn_count_leading_zeros (dp[dsize-1]) - GMP_NAIL_BITS;
1478
1725
 
1479
1726
      n2p[nsize] = refmpn_lshift_or_copy (n2p, np, nsize, norm);
1480
1727
      ASSERT_NOCARRY (refmpn_lshift_or_copy (d2p, dp, dsize, norm));
1488
1735
    }
1489
1736
}
1490
1737
 
 
1738
 
1491
1739
size_t
1492
1740
refmpn_get_str (unsigned char *dst, int base, mp_ptr src, mp_size_t size)
1493
1741
{
1502
1750
 
1503
1751
  MPN_SIZEINBASE (dsize, src, size, base);
1504
1752
  ASSERT (dsize >= 1);
1505
 
  ASSERT (! byte_overlap_p (dst, dsize, src, size * BYTES_PER_MP_LIMB));
 
1753
  ASSERT (! byte_overlap_p (dst, (mp_size_t) dsize, src, size * BYTES_PER_MP_LIMB));
1506
1754
 
1507
1755
  if (size == 0)
1508
1756
    {
1519
1767
    {
1520
1768
      d--;
1521
1769
      ASSERT (d >= dst);
1522
 
      *d = refmpn_divrem_1 (src, 0, src, size, (mp_limb_t) base);
 
1770
      *d = refmpn_divrem_1 (src, (mp_size_t) 0, src, size, (mp_limb_t) base);
1523
1771
      size -= (src[size-1] == 0);
1524
1772
    }
1525
1773
  while (size != 0);
1526
1774
 
1527
 
  /* at most one leading zero */
1528
 
  ASSERT (d == dst || d == dst+1);
 
1775
  /* Move result back and decrement dsize if we didn't generate
 
1776
     the maximum possible digits.  */
1529
1777
  if (d != dst)
1530
 
    *dst = 0;
 
1778
    {
 
1779
      size_t i;
 
1780
      dsize -= d - dst;
 
1781
      for (i = 0; i < dsize; i++)
 
1782
        dst[i] = d[i];
 
1783
    }
1531
1784
 
1532
1785
  if (POW2_P (base))
1533
1786
    free (src);
1537
1790
 
1538
1791
 
1539
1792
mp_limb_t
1540
 
refmpn_bswap_limb (mp_limb_t src)
 
1793
ref_bswap_limb (mp_limb_t src)
1541
1794
{
1542
1795
  mp_limb_t  dst;
1543
1796
  int        i;
1635
1888
      bit = GMP_NUMB_HIGHBIT;
1636
1889
    }
1637
1890
}
 
1891
 
 
1892
/* This is a simple bitwise algorithm working high to low across "s" and
 
1893
   testing each time whether setting the bit would make s^2 exceed n.  */
 
1894
mp_size_t
 
1895
refmpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr np, mp_size_t nsize)
 
1896
{
 
1897
  mp_ptr     tp, dp;
 
1898
  mp_size_t  ssize, talloc, tsize, dsize, ret, ilimbs;
 
1899
  unsigned   ibit;
 
1900
  long       i;
 
1901
  mp_limb_t  c;
 
1902
 
 
1903
  ASSERT (nsize >= 0);
 
1904
 
 
1905
  /* If n==0, then s=0 and r=0.  */
 
1906
  if (nsize == 0)
 
1907
    return 0;
 
1908
 
 
1909
  ASSERT (np[nsize - 1] != 0);
 
1910
  ASSERT (rp == NULL || MPN_SAME_OR_SEPARATE_P (np, rp, nsize));
 
1911
  ASSERT (rp == NULL || ! MPN_OVERLAP_P (sp, (nsize + 1) / 2, rp, nsize));
 
1912
  ASSERT (! MPN_OVERLAP_P (sp, (nsize + 1) / 2, np, nsize));
 
1913
 
 
1914
  /* root */
 
1915
  ssize = (nsize+1)/2;
 
1916
  refmpn_zero (sp, ssize);
 
1917
 
 
1918
  /* the remainder so far */
 
1919
  dp = refmpn_memdup_limbs (np, nsize);
 
1920
  dsize = nsize;
 
1921
 
 
1922
  /* temporary */
 
1923
  talloc = 2*ssize + 1;
 
1924
  tp = refmpn_malloc_limbs (talloc);
 
1925
 
 
1926
  for (i = GMP_NUMB_BITS * ssize - 1; i >= 0; i--)
 
1927
    {
 
1928
      /* t = 2*s*2^i + 2^(2*i), being the amount s^2 will increase by if 2^i
 
1929
         is added to it */
 
1930
 
 
1931
      ilimbs = (i+1) / GMP_NUMB_BITS;
 
1932
      ibit = (i+1) % GMP_NUMB_BITS;
 
1933
      refmpn_zero (tp, ilimbs);
 
1934
      c = refmpn_lshift_or_copy (tp+ilimbs, sp, ssize, ibit);
 
1935
      tsize = ilimbs + ssize;
 
1936
      tp[tsize] = c;
 
1937
      tsize += (c != 0);
 
1938
 
 
1939
      ilimbs = (2*i) / GMP_NUMB_BITS;
 
1940
      ibit = (2*i) % GMP_NUMB_BITS;
 
1941
      if (ilimbs + 1 > tsize)
 
1942
        {
 
1943
          refmpn_zero_extend (tp, tsize, ilimbs + 1);
 
1944
          tsize = ilimbs + 1;
 
1945
        }
 
1946
      c = refmpn_add_1 (tp+ilimbs, tp+ilimbs, tsize-ilimbs,
 
1947
                        CNST_LIMB(1) << ibit);
 
1948
      ASSERT (tsize < talloc);
 
1949
      tp[tsize] = c;
 
1950
      tsize += (c != 0);
 
1951
 
 
1952
      if (refmpn_cmp_twosizes (dp, dsize, tp, tsize) >= 0)
 
1953
        {
 
1954
          /* set this bit in s and subtract from the remainder */
 
1955
          refmpn_setbit (sp, i);
 
1956
 
 
1957
          ASSERT_NOCARRY (refmpn_sub_n (dp, dp, tp, dsize));
 
1958
          dsize = refmpn_normalize (dp, dsize);
 
1959
        }
 
1960
    }
 
1961
 
 
1962
  if (rp == NULL)
 
1963
    {
 
1964
      ret = ! refmpn_zero_p (dp, dsize);
 
1965
    }
 
1966
  else
 
1967
    {
 
1968
      ASSERT (dsize == 0 || dp[dsize-1] != 0);
 
1969
      refmpn_copy (rp, dp, dsize);
 
1970
      ret = dsize;
 
1971
    }
 
1972
 
 
1973
  free (dp);
 
1974
  free (tp);
 
1975
  return ret;
 
1976
}