Changeset 31390 in project


Ignore:
Timestamp:
09/11/14 21:35:58 (7 years ago)
Author:
sjamaan
Message:

numbers: Convert part of big-quotient-fix (ugly medium and small denominator devision algorithms) to core naming conventions. Still a lot remains to be done! (but at least the ugly #defines are gone from divide_digit_subtract, at the cost of some code repetition) Performance is a tiny bit better, too.

Location:
release/4/numbers/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • release/4/numbers/trunk/numbers-c.c

    r31376 r31390  
    10701070
    10711071static void
    1072 big_quotient_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
    1073 {
    1074   bignum_type bigx = big_of(x);
    1075   y = C_unfix(y);
    1076 
    1077   if (y == 1)
    1078     C_kontinue(k, x);
    1079   else if (y == -1)
    1080     C_return_bignum(k, bignum_new_sign(bigx, !(BIGNUM_NEGATIVE_P(bigx))));
    1081 
    1082   /* Too bad, we really need to do some work... */
    1083   {
    1084     int neg_p = (y < 0) ? !(BIGNUM_NEGATIVE_P(bigx)) : BIGNUM_NEGATIVE_P(bigx);
    1085     bignum_digit_type abs_y = (y < 0) ? -y : y;
    1086     bignum_type quotient;
    1087    
    1088     if (y == C_MOST_NEGATIVE_FIXNUM) {
    1089       if (!BIGNUM_NEGATIVE_P(bigx) && BIGNUM_LENGTH(bigx) == 1
    1090           && BIGNUM_REF(bigx, 1) == 1 && BIGNUM_REF(bigx, 0) == 0) {
    1091         /*
    1092          * Very very special case:
    1093          * quotient(MOST_NEGATIVE_FIXNUM, -(MOST_NEGATIVE_FIXNUM)) => -1
    1094          */
    1095         C_kontinue(k, C_fix(-1));
    1096       } else {
    1097         /* This is the only case we need to go allocate a bignum for */
    1098         bignum_type bigy =
    1099           bignum_allocate_from_fixnum(C_fix(C_MOST_NEGATIVE_FIXNUM));
    1100 
    1101         bignum_divide_unsigned_large_denominator
    1102           (bigx, bigy, (&quotient), ((bignum_type *) 0), neg_p, 0);
    1103         BIGNUM_DEALLOCATE(bigy);
    1104         C_return_bignum(k, quotient);
    1105       }
    1106     } else if (abs_y < BIGNUM_RADIX_ROOT) {
    1107       bignum_divide_unsigned_small_denominator
    1108         (bigx, abs_y, (&quotient), ((bignum_type *) 0), neg_p, 0);
    1109       C_return_bignum(k, quotient);
    1110     } else {
    1111       bignum_divide_unsigned_medium_denominator
    1112         (bigx, abs_y, (&quotient), ((bignum_type *) 0), neg_p, 0);
    1113       C_return_bignum(k, quotient);
    1114     }
    1115   }
    1116 }
    1117 
    1118 static void
    11191072big_quotient_big(C_word c, C_word self, C_word k, C_word x, C_word y)
    11201073{
     
    12471200static void bignum_posneg_bitwise_op(C_word c, C_word self, C_word result);
    12481201static void bignum_pospos_bitwise_op(C_word c, C_word self, C_word result);
     1202static void bignum_destructive_normalize(C_word target, C_word source, C_word shift_left);
     1203static void bignum_destructive_divide_unsigned_halfdigit(C_word c, C_word self, C_word quotient);
     1204static void bignum_destructive_divide_unsigned_digit(C_word c, C_word self, C_word quotient);
     1205static C_word bignum_divide_digit(C_word uh, C_word ul, C_word v, C_word *q);
     1206static C_word bignum_divide_digit_subtract(C_word v1, C_word v2, C_word guess, C_word *u);
    12491207
    12501208/* Eventually this will probably need to be integrated into C_2_plus. */
     
    16041562}
    16051563
     1564/* XXX TODO: Maybe pass true/false/bignum as initp, to allow copying data? */
    16061565void C_ccall
    16071566C_allocate_bignum(C_word c, C_word self, C_word k, C_word size, C_word negp, C_word initp)
     
    16951654/* Given (denominator > 1), it is fairly easy to show that
    16961655   quotient_high fits a bignum halfdigit, after which it is easy to
    1697    see that all digits fit a bignum full digit. */
     1656   see that all digits fit a bignum full digit.
     1657   
     1658   This works because denominator is known to fit a halfdigit, which
     1659   means that the remainder modulo denominator will also fit a halfdigit. */
    16981660static C_word
    16991661bignum_digits_destructive_scale_down(C_word *start, C_word *end, C_word denominator)
     
    18231785        ((R * R) + (R * (R - 2)) + (R - 1))
    18241786   Maximum value for carry: ((R * (R - 1)) + (R - 1))
    1825         where R == 2^HALF_DIGIT_SHIFT */
     1787        where R == 2^HALF_DIGIT_LENGTH */
    18261788static void
    18271789bignum_times_bignum_unsigned(C_word k, C_word x, C_word y, C_word negp)
     
    25602522  C_kontinue(k, C_bignum_normalize(result));
    25612523}
     2524
     2525void C_ccall
     2526C_u_bignum_quotient_fixnum(C_word c, C_word self, C_word k, C_word x, C_word y)
     2527{
     2528  y = C_unfix(y);
     2529 
     2530  if (y == 1) {
     2531    C_kontinue(k, x);
     2532  } else if (y == -1) {
     2533    C_u_bignum_negate(1, (C_word)NULL, k, x);
     2534  } else if (y == C_MOST_NEGATIVE_FIXNUM) {
     2535    /* This is the only case we need to go allocate a bignum for */
     2536    C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
     2537    bignum_type quotient;
     2538
     2539    y = C_a_u_i_fix_to_big(&a, C_fix(C_MOST_NEGATIVE_FIXNUM));
     2540
     2541    /* XXX TODO */
     2542    bignum_divide_unsigned_large_denominator
     2543      (big_of(x), big_of(y), (&quotient), (bignum_type)NULL, !C_bignum_negativep(x), 0);
     2544    C_return_bignum(k, quotient);
     2545  } else {
     2546    C_word negp = C_mk_bool((y < 0) ? !C_bignum_negativep(x) :
     2547                                      C_bignum_negativep(x)),
     2548           absy = (y < 0) ? -y : y,
     2549           kab[C_SIZEOF_CLOSURE(7)], *ka = kab, k2,
     2550           size, func;
     2551
     2552    if (C_fitsinbignumhalfdigitp(absy)) {
     2553      size = C_fix(C_bignum_size(x));
     2554      func = (C_word)bignum_destructive_divide_unsigned_halfdigit;
     2555    } else {
     2556      size = C_fix(C_bignum_size(x)) + 1; /* Due to normalization */
     2557      func = (C_word)bignum_destructive_divide_unsigned_digit;
     2558    }
     2559
     2560    k2 = C_closure(&ka, 5, func, k, x, C_fix(absy),
     2561                   /* Return quotient, not remainder */
     2562                   C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE);
     2563    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
     2564  }
     2565}
     2566
     2567static void
     2568bignum_destructive_divide_unsigned_halfdigit(C_word c, C_word self, C_word quotient)
     2569{
     2570  C_word k = C_block_item(self, 1),
     2571         numerator = C_block_item(self, 2),
     2572         denominator = C_unfix(C_block_item(self, 3)),
     2573         /* return_quotient = C_block_item(self, 4), */
     2574         return_remainder = C_block_item(self, 5),
     2575         remainder_negp = C_block_item(self, 6),
     2576         *start = C_bignum_digits(quotient),
     2577         *end = start + C_bignum_size(quotient) - 1,
     2578         remainder;
     2579
     2580  C_memcpy(start, C_bignum_digits(numerator),
     2581           /* TODO: This is currently in bytes.  If we change the
     2582            * representation that needs to change!
     2583            * We subtract the size of the header, too.
     2584            */
     2585           C_header_size(C_internal_bignum(numerator))-C_wordstobytes(1));
     2586
     2587  remainder = bignum_digits_destructive_scale_down(start, end, denominator);
     2588  assert(C_fitsinbignumhalfdigitp(remainder));
     2589
     2590  C_bignum_destructive_trim(quotient);
     2591  quotient = C_bignum_normalize(quotient);
     2592
     2593  if (C_truep(return_remainder)) {
     2594    remainder = C_truep(remainder_negp) ? -remainder : remainder;
     2595    C_values(4, C_SCHEME_UNDEFINED, k, quotient, C_fix(remainder));
     2596  } else {
     2597    C_kontinue(k, quotient);
     2598  }
     2599}
     2600
     2601/* XXX TODO: This seems unnecessary */
     2602static void
     2603bignum_destructive_divide_unsigned_digit(C_word c, C_word self, C_word quotient)
     2604{
     2605  C_word k = C_block_item(self, 1),
     2606         numerator = C_block_item(self, 2),
     2607         denominator = C_unfix(C_block_item(self, 3)),
     2608         return_quotient = C_block_item(self, 4),
     2609         return_remainder = C_block_item(self, 5),
     2610         remainder_negp = C_block_item(self, 6),
     2611         *start = C_bignum_digits(quotient),
     2612         *scan = start + C_bignum_size(quotient),
     2613         qj, remainder = 0, shift = 0;
     2614
     2615  /* Because `bignum_divide_digit' requires a normalized denominator. */
     2616  while (denominator < ((C_word)1 << (C_BIGNUM_DIGIT_LENGTH-1))) {
     2617    denominator <<= 1;
     2618    shift++;
     2619  }
     2620
     2621  if (shift) {
     2622    bignum_destructive_normalize(quotient, numerator, shift);
     2623  } else {
     2624    /* Could alo use destructive_normalize here, but this is faster */
     2625    C_memcpy(start, C_bignum_digits(numerator),
     2626             /* TODO: This is currently in bytes.  If we change the
     2627              * representation that needs to change!
     2628              * We subtract the size of the header, too.
     2629              */
     2630             C_header_size(C_internal_bignum(numerator))-C_wordstobytes(1));
     2631
     2632    *(scan-1) = 0; /* Most significant digit, not copied */
     2633  }
     2634
     2635  if (C_truep(return_quotient)) {
     2636    while (start < scan) {
     2637      remainder = bignum_divide_digit(remainder, (*--scan), denominator, &qj);
     2638      (*scan) = qj;
     2639    }
     2640    C_bignum_destructive_trim(quotient);
     2641    quotient = C_bignum_normalize(quotient);
     2642  } else {
     2643    while (start < scan)
     2644      remainder = bignum_divide_digit(remainder, (*--scan), denominator, &qj);
     2645  }
     2646
     2647  if (C_truep(return_remainder)) {
     2648    remainder >>= shift;
     2649    remainder = C_truep(remainder_negp) ? -remainder : remainder;
     2650
     2651    if (C_truep(return_quotient))
     2652      C_values(4, C_SCHEME_UNDEFINED, k, quotient, C_fix(remainder));
     2653    else
     2654      C_kontinue(k, C_fix(remainder));
     2655  } else {
     2656    assert(C_truep(return_quotient));
     2657    C_kontinue(k, quotient);
     2658  }
     2659}
     2660
     2661static void
     2662bignum_destructive_normalize(C_word target, C_word source, C_word shift_left)
     2663{
     2664  C_word digit, carry = 0,
     2665         *scan_source = C_bignum_digits(source),
     2666         *scan_target = C_bignum_digits(target),
     2667         *end_source = scan_source + C_bignum_size(source),
     2668         *end_target = scan_target + C_bignum_size(target),
     2669         shift_right = C_BIGNUM_DIGIT_LENGTH - shift_left,
     2670         mask = (1L << shift_right) - 1;
     2671
     2672  while (scan_source < end_source) {
     2673    digit = (*scan_source++);
     2674    (*scan_target++) = (((digit & mask) << shift_left) | carry);
     2675    carry = (digit >> shift_right);
     2676  }
     2677  if (scan_target < end_target)
     2678    (*scan_target) = carry;
     2679  else
     2680    assert(carry == 0);
     2681}
     2682
     2683/* TODO: This pointer stuff here is suspicious: it's probably slow */
     2684static C_word
     2685bignum_divide_digit_subtract(C_word v1, C_word v2, C_word guess, C_word *u)
     2686{
     2687  C_word product, diff, sum, carry;
     2688
     2689  product = v2 * guess;
     2690  diff = u[2] - C_BIGNUM_DIGIT_LO_HALF(product);
     2691  if (diff < 0) {
     2692    u[2] = diff + ((C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH);
     2693    carry = C_BIGNUM_DIGIT_HI_HALF(product) + 1;
     2694  } else {
     2695    u[2] = diff;
     2696    carry = C_BIGNUM_DIGIT_HI_HALF(product);
     2697  }
     2698 
     2699  product = v1 * guess + carry;
     2700  diff = u[1] - C_BIGNUM_DIGIT_LO_HALF(product);
     2701  if (diff < 0) {
     2702    u[1] = diff + ((C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH);
     2703    carry = C_BIGNUM_DIGIT_HI_HALF(product) + 1;
     2704  } else {
     2705    u[1] = diff;
     2706    carry = C_BIGNUM_DIGIT_HI_HALF(product);
     2707    if (carry == 0) return guess; /* DONE */
     2708  }
     2709
     2710  diff = u[0] - carry;
     2711  if (diff < 0) {
     2712    u[0] = diff + ((C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH);
     2713  } else {
     2714    u[0] = diff;
     2715    return guess; /* DONE */
     2716  }
     2717
     2718  sum = v2 + u[2];
     2719  u[2] = sum & C_BIGNUM_HALF_DIGIT_MASK;
     2720  carry = C_fitsinbignumhalfdigitp(sum);
     2721
     2722  sum = v1 + u[1] + carry;
     2723  u[1] = sum & C_BIGNUM_HALF_DIGIT_MASK;
     2724  carry = C_fitsinbignumhalfdigitp(sum);
     2725
     2726  u[0] += carry;
     2727
     2728  return guess - 1;
     2729}
     2730
     2731/* This is a reduced version of the division algorithm, applied to the
     2732   case of dividing two bignum digits by one bignum digit.  It is
     2733   assumed that the numerator, denominator are normalized. */
     2734
     2735#define BDD_STEP(qn, j)                                                 \
     2736{                                                                       \
     2737  uj = u[j];                                                            \
     2738  if (uj != v1) {                                                       \
     2739    uj_uj1 = C_BIGNUM_DIGIT_COMBINE(uj, u[j + 1]);                      \
     2740    guess = uj_uj1 / v1;                                                \
     2741    comparand = C_BIGNUM_DIGIT_COMBINE(uj_uj1 % v1, u[j + 2]);          \
     2742  } else {                                                              \
     2743    guess = C_BIGNUM_HALF_DIGIT_MASK;                                   \
     2744    comparand = C_BIGNUM_DIGIT_COMBINE(u[j + 1] + v1, u[j + 2]);        \
     2745  }                                                                     \
     2746  while ((guess * v2) > comparand) {                                    \
     2747    guess -= 1;                                                         \
     2748    comparand += v1 << C_BIGNUM_HALF_DIGIT_LENGTH;                      \
     2749    if (!C_fitsinbignumdigitp(comparand))                               \
     2750      break;                                                            \
     2751  }                                                                     \
     2752  qn = bignum_divide_digit_subtract(v1, v2, guess, &u[j]);              \
     2753}
     2754
     2755static C_word
     2756bignum_divide_digit(C_word uh, C_word ul, C_word v, C_word *q)
     2757{
     2758  C_word guess, comparand, v1, v2, uj, uj_uj1, q1, q2, u[4];
     2759
     2760  if (uh == 0) {
     2761    if (ul < v) {
     2762      *q = 0;
     2763      return ul;
     2764    } else if (ul == v) {
     2765      *q = 1;
     2766      return 0;
     2767    }
     2768  }
     2769  u[0] = C_BIGNUM_DIGIT_HI_HALF(uh);
     2770  u[1] = C_BIGNUM_DIGIT_LO_HALF(uh);
     2771  u[2] = C_BIGNUM_DIGIT_HI_HALF(ul);
     2772  u[3] = C_BIGNUM_DIGIT_LO_HALF(ul);
     2773  v1 = C_BIGNUM_DIGIT_HI_HALF(v);
     2774  v2 = C_BIGNUM_DIGIT_LO_HALF(v);
     2775  BDD_STEP(q1, 0);
     2776  BDD_STEP(q2, 1);
     2777  *q = C_BIGNUM_DIGIT_COMBINE(q1, q2);
     2778  return C_BIGNUM_DIGIT_COMBINE(u[2], u[3]);
     2779}
     2780
     2781#undef BDD_STEP
  • release/4/numbers/trunk/numbers-c.h

    r31379 r31390  
    149149
    150150#ifdef C_SIXTY_FOUR
    151 # define C_BIGNUM_DIGIT_LENGTH     62
    152 # define C_BIGNUM_HEADER_SIGN_BIT  0x4000000000000000L
    153 # define C_BIGNUM_HEADER_SIZE_MASK 0x3fffffffffffffffL
    154 # define C_BIGNUM_DIGIT_MASK       0x3fffffffffffffffL
    155 # define C_BIGNUM_HALF_DIGIT_MASK  0x000000007fffffffL
    156 # define C_BIGNUM_HALF_DIGIT_SHIFT 31
     151# define C_BIGNUM_DIGIT_LENGTH      62
     152# define C_BIGNUM_HEADER_SIGN_BIT   0x4000000000000000L
     153# define C_BIGNUM_HEADER_SIZE_MASK  0x3fffffffffffffffL
     154# define C_BIGNUM_DIGIT_MASK        0x3fffffffffffffffL
     155# define C_BIGNUM_HALF_DIGIT_MASK   0x000000007fffffffL
     156# define C_BIGNUM_HALF_DIGIT_LENGTH 31
    157157#else
    158 # define C_BIGNUM_DIGIT_LENGTH     30
    159 # define C_BIGNUM_HEADER_SIGN_BIT  0x40000000
    160 # define C_BIGNUM_HEADER_SIZE_MASK 0x3fffffff
    161 # define C_BIGNUM_DIGIT_MASK       0x3fffffff
    162 # define C_BIGNUM_HALF_DIGIT_MASK  0x00007fff
    163 # define C_BIGNUM_HALF_DIGIT_SHIFT 15
     158# define C_BIGNUM_DIGIT_LENGTH      30
     159# define C_BIGNUM_HEADER_SIGN_BIT   0x40000000
     160# define C_BIGNUM_HEADER_SIZE_MASK  0x3fffffff
     161# define C_BIGNUM_DIGIT_MASK        0x3fffffff
     162# define C_BIGNUM_HALF_DIGIT_MASK   0x00007fff
     163# define C_BIGNUM_HALF_DIGIT_LENGTH 15
    164164#endif
    165165
     
    168168
    169169#define C_BIGNUM_DIGIT_LO_HALF(d)       ((d) & C_BIGNUM_HALF_DIGIT_MASK)
    170 #define C_BIGNUM_DIGIT_HI_HALF(d)       ((d) >> C_BIGNUM_HALF_DIGIT_SHIFT)
    171 #define C_BIGNUM_DIGIT_COMBINE(h,l)     ((h) << C_BIGNUM_HALF_DIGIT_SHIFT|(l))
     170#define C_BIGNUM_DIGIT_HI_HALF(d)       ((d) >> C_BIGNUM_HALF_DIGIT_LENGTH)
     171#define C_BIGNUM_DIGIT_COMBINE(h,l)     ((h) << C_BIGNUM_HALF_DIGIT_LENGTH|(l))
    172172
    173173#define C_fitsinbignumdigitp(n)         ((C_uword)(n) == ((C_uword)(n) & C_BIGNUM_DIGIT_MASK))
     
    204204void C_ccall C_u_2_bignum_times(C_word c, C_word self, C_word k, C_word x, C_word y);
    205205
     206void C_ccall C_u_bignum_quotient_fixnum(C_word c, C_word self, C_word k, C_word x, C_word y);
     207
    206208void C_ccall C_u_fixnum_minus_bignum(C_word c, C_word self, C_word k, C_word x, C_word y);
    207209void C_ccall C_u_bignum_minus_fixnum(C_word c, C_word self, C_word k, C_word x, C_word y);
  • release/4/numbers/trunk/numbers.scm

    r31376 r31390  
    151151(define %big*big (##core#primitive "C_u_2_bignum_times"))
    152152
    153 (define %big-quotient-fix (##core#primitive "big_quotient_fix"))
     153(define %big-quotient-fix (##core#primitive "C_u_bignum_quotient_fixnum"))
    154154(define %big-quotient-big (##core#primitive "big_quotient_big"))
    155155
Note: See TracChangeset for help on using the changeset viewer.