Changeset 31411 in project


Ignore:
Timestamp:
09/12/14 17:31:51 (5 years ago)
Author:
sjamaan
Message:

numbers: Convert quotient fully to core naming conventions

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

Legend:

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

    r31407 r31411  
    10341034
    10351035static void
    1036 big_quotient_big(C_word c, C_word self, C_word k, C_word x, C_word y)
    1037 {
    1038   bignum_type numerator = big_of(x);
    1039   bignum_type denominator = big_of(y);
    1040   int neg_p = ((BIGNUM_NEGATIVE_P (denominator))
    1041                ? (! (BIGNUM_NEGATIVE_P (numerator)))
    1042                : (BIGNUM_NEGATIVE_P (numerator)));
    1043 
    1044   switch (bignum_compare_unsigned (numerator, denominator))
    1045     {
    1046     case bignum_comparison_equal:
    1047       C_kontinue(k, neg_p ? C_fix(-1) : C_fix(1));
    1048     case bignum_comparison_less:
    1049       C_kontinue(k, C_fix(0));
    1050     case bignum_comparison_greater:
    1051     default:                                  /* to appease gcc -Wall */
    1052       {
    1053         bignum_type quotient;
    1054         assert(BIGNUM_LENGTH(denominator) > 1);
    1055         bignum_divide_unsigned_large_denominator
    1056           (numerator, denominator, (&quotient), ((bignum_type *) 0), neg_p, 0);
    1057         C_return_bignum(k, quotient);
    1058       }
    1059     }
    1060 }
    1061 
    1062 static void
    10631036big_remainder_fix(C_word c, C_word self, C_word k, C_word x, C_word y)
    10641037{
     
    11681141static void bignum_destructive_divide_unsigned_digit(C_word c, C_word self, C_word quotient);
    11691142static C_word bignum_divide_digit(C_word uh, C_word ul, C_word v, C_word *q);
    1170 static C_word bignum_divide_digit_subtract(C_word v1, C_word v2, C_word guess, C_word *u);
     1143static C_word bignum_divide_and_subtract_digit(C_word v1, C_word v2, C_word guess, C_word *u);
     1144static void bignum_divide_2_unsigned(C_word c, C_word self, C_word quotient);
     1145static void bignum_divide_2_unsigned_2(C_word c, C_word self, C_word remainder);
     1146static void bignum_divide_2_unsigned_3(C_word c, C_word self, C_word tmp_big);
     1147static void bignum_destructive_divide_normalized(C_word u, C_word v, C_word q);
     1148static C_word bignum_divide_and_subtract(C_word *v_start, C_word *v_end, C_word guess, C_word *u_start);
     1149static C_word bignum_normalize_shifted(C_word bignum, C_word shift_right);
    11711150
    11721151/* Copy all the digits from source to target, obliterating what was
     
    24732452  } else if (y == C_MOST_NEGATIVE_FIXNUM) {
    24742453    /* This is the only case we need to go allocate a bignum for */
    2475     C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
    2476     bignum_type quotient;
     2454    C_word ab[C_SIZEOF_FIX_BIGNUM+C_SIZEOF_CLOSURE(9)], *a = ab, k2,
     2455           negp = C_mk_bool(!C_bignum_negativep(x)), size;
    24772456
    24782457    y = C_a_u_i_fix_to_big(&a, C_fix(C_MOST_NEGATIVE_FIXNUM));
    24792458
    2480     /* XXX TODO */
    2481     bignum_divide_unsigned_large_denominator
    2482       (big_of(x), big_of(y), (&quotient), (bignum_type)NULL, !C_bignum_negativep(x), 0);
    2483     C_return_bignum(k, quotient);
     2459    k2 = C_closure(&a, 9, (C_word)bignum_divide_2_unsigned, k, x, y,
     2460                   /* Return quotient, not remainder */
     2461                   C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE,
     2462                   /* Will be filled in later */
     2463                   C_SCHEME_UNDEFINED, C_SCHEME_UNDEFINED);
     2464    size = C_fix(C_bignum_size(x) + 1 - C_bignum_size(y));
     2465    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
    24842466  } else {
    24852467    C_word negp = C_mk_bool((y < 0) ? !C_bignum_negativep(x) :
     
    24972479    }
    24982480
    2499     k2 = C_closure(&ka, 5, func, k, x, C_fix(absy),
     2481    k2 = C_closure(&ka, 7, func, k, x, C_fix(absy),
    25002482                   /* Return quotient, not remainder */
    25012483                   C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE);
     2484    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
     2485  }
     2486}
     2487
     2488void C_ccall
     2489C_u_bignum_quotient_bignum(C_word c, C_word self, C_word k, C_word x, C_word y)
     2490{
     2491  C_word ab[C_SIZEOF_FIX_BIGNUM+C_SIZEOF_CLOSURE(9)], *a = ab, k2, size,
     2492         negp = C_mk_bool(C_bignum_negativep(x) ?
     2493                         !C_bignum_negativep(y) :
     2494                         C_bignum_negativep(y));
     2495
     2496  switch (bignum_cmp_unsigned(x, y)) {
     2497  case 0:
     2498    C_kontinue(k, C_truep(negp) ? C_fix(-1) : C_fix(1));
     2499  case -1:
     2500    C_kontinue(k, C_fix(0));
     2501  case 1:
     2502    k2 = C_closure(&a, 9, (C_word)bignum_divide_2_unsigned, k, x, y,
     2503                   /* Return quotient, not remainder */
     2504                   C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE,
     2505                   /* Will be filled in later */
     2506                   C_SCHEME_UNDEFINED, C_SCHEME_UNDEFINED);
     2507    size = C_fix(C_bignum_size(x) + 1 - C_bignum_size(y));
    25022508    C_allocate_bignum(3, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
    25032509  }
     
    26092615/* TODO: This pointer stuff here is suspicious: it's probably slow */
    26102616static C_word
    2611 bignum_divide_digit_subtract(C_word v1, C_word v2, C_word guess, C_word *u)
     2617bignum_divide_and_subtract_digit(C_word v1, C_word v2, C_word guess, C_word *u)
    26122618{
    26132619  C_word product, diff, sum, carry;
     
    26762682      break;                                                            \
    26772683  }                                                                     \
    2678   qn = bignum_divide_digit_subtract(v1, v2, guess, &u[j]);              \
    2679 }
    2680 
     2684  qn = bignum_divide_and_subtract_digit(v1, v2, guess, &u[j]);          \
     2685}
     2686
     2687/* Because the algorithm from Knuth requires combining the two highest
     2688 * digits of u (which we can't fit in a C_word), we instead do two
     2689 * such steps, a halfdigit at a time.
     2690 */
    26812691static C_word
    26822692bignum_divide_digit(C_word uh, C_word ul, C_word v, C_word *q)
     
    27062716
    27072717#undef BDD_STEP
     2718
     2719/* Full bignum division */
     2720
     2721static void
     2722bignum_divide_2_unsigned(C_word c, C_word self, C_word quotient)
     2723{
     2724  C_word numerator = C_block_item(self, 2),
     2725         remainder_negp = C_block_item(self, 6),
     2726         size = C_fix(C_bignum_size(numerator) + 1);
     2727
     2728  /* Nice: We can recycle the current closure */
     2729  C_set_block_item(self, 0, (C_word)bignum_divide_2_unsigned_2);
     2730  C_set_block_item(self, 7, quotient);
     2731  C_allocate_bignum(3, (C_word)NULL, self, size, remainder_negp, C_SCHEME_FALSE);
     2732}
     2733
     2734/* For help understanding this algorithm, see:
     2735   Knuth, Donald E., "The Art of Computer Programming",
     2736   volume 2, "Seminumerical Algorithms"
     2737   section 4.3.1, "Multiple-Precision Arithmetic". */
     2738
     2739static void
     2740bignum_divide_2_unsigned_2(C_word c, C_word self, C_word remainder)
     2741{
     2742  C_word k = C_block_item(self, 1),
     2743         numerator = C_block_item(self, 2),
     2744         denominator = C_block_item(self, 3),
     2745         return_quotient = C_block_item(self, 4),
     2746         return_remainder = C_block_item(self, 5),
     2747         /* This one may be overwritten with the remainder */
     2748         /* remainder_negp = C_block_item(self, 6), */
     2749         quotient = C_block_item(self, 7),
     2750         length_d = C_bignum_size(denominator),
     2751         d1 = *(C_bignum_digits(denominator) + length_d - 1),
     2752         shift = 0;
     2753
     2754  assert(length_d > 1);
     2755  while (d1 < (BIGNUM_RADIX / 2)) {
     2756    d1 <<= 1;
     2757    shift++;
     2758  }
     2759
     2760  if (shift == 0) { /* Already normalized */
     2761    bignum_digits_destructive_copy(remainder, numerator);
     2762    /* Set most significant digit */
     2763    *(C_bignum_digits(remainder) + C_bignum_size(numerator)) = 0;
     2764
     2765    bignum_destructive_divide_normalized(remainder, denominator, quotient);
     2766
     2767    if (C_truep(C_and(return_quotient, return_remainder))) {
     2768      C_values(4, C_SCHEME_UNDEFINED, k,
     2769               C_bignum_normalize(quotient), C_bignum_normalize(remainder));
     2770    } else if (C_truep(return_remainder)) {
     2771      C_kontinue(k, C_bignum_normalize(remainder));
     2772    } else {
     2773      assert(C_truep(return_quotient));
     2774      C_kontinue(k, C_bignum_normalize(quotient));
     2775    }
     2776  } else {
     2777    /* Requires normalisation of denominator;  Allocate temp bignum for it. */
     2778    C_set_block_item(self, 0, (C_word)bignum_divide_2_unsigned_3);
     2779    C_set_block_item(self, 6, remainder);
     2780    C_set_block_item(self, 8, C_fix(shift));
     2781    C_allocate_bignum(3, (C_word)NULL, self, C_fix(length_d),
     2782                      C_SCHEME_FALSE, C_SCHEME_FALSE);
     2783  }
     2784}
     2785
     2786static void
     2787bignum_divide_2_unsigned_3(C_word c, C_word self, C_word tmp_big)
     2788{
     2789  C_word k = C_block_item(self, 1),
     2790         numerator = C_block_item(self, 2),
     2791         denominator = C_block_item(self, 3),
     2792         return_quotient = C_block_item(self, 4),
     2793         return_remainder = C_block_item(self, 5),
     2794         remainder = C_block_item(self, 6),
     2795         quotient = C_block_item(self, 7),
     2796         shift = C_unfix(C_block_item(self, 8));
     2797
     2798  bignum_destructive_normalize(remainder, numerator, shift);
     2799  bignum_destructive_normalize(tmp_big, denominator, shift);
     2800  bignum_destructive_divide_normalized(remainder, tmp_big, quotient);
     2801
     2802  if (C_truep(return_remainder)) {
     2803    if (C_truep(return_quotient)) {
     2804      C_values(4, C_SCHEME_UNDEFINED, k,
     2805               C_bignum_normalize(quotient),
     2806               bignum_normalize_shifted(remainder, shift));
     2807    } else {
     2808      C_kontinue(k, bignum_normalize_shifted(remainder, shift));
     2809    }
     2810  } else {
     2811    assert(C_truep(return_quotient));
     2812    C_kontinue(k, C_bignum_normalize(quotient));
     2813  }
     2814}
     2815
     2816static void
     2817bignum_destructive_divide_normalized(C_word u, C_word v, C_word q)
     2818{
     2819  C_word u_length = C_bignum_size(u),
     2820         v_length = C_bignum_size(v),
     2821         *u_start = C_bignum_digits(u),
     2822         *u_scan = u_start + u_length,
     2823         *u_scan_limit = u_start + v_length,
     2824         *u_scan_start = u_scan - v_length,
     2825         *v_start = C_bignum_digits(v),
     2826         *v_end = v_start + v_length,
     2827         *q_scan = (q == C_SCHEME_UNDEFINED) ? NULL :
     2828                   (C_bignum_digits(q) + C_bignum_size(q)),
     2829         v1 = v_end[-1],
     2830         v2 = v_end[-2],
     2831         ph, /* high half of double-digit product */
     2832         pl, /* low half of double-digit product */
     2833         guess, uj, qj,
     2834         gh, /* high half-digit of guess */
     2835         ch, /* high half of double-digit comparand */
     2836         v2l = C_BIGNUM_DIGIT_LO_HALF(v2),
     2837         v2h = C_BIGNUM_DIGIT_HI_HALF(v2),
     2838         cl, /* low half of double-digit comparand */
     2839         gl, /* low half-digit of guess */
     2840         gm; /* memory loc for reference parameter */
     2841
     2842  while (u_scan_limit < u_scan) {
     2843    uj = (*--u_scan);
     2844    if (uj != v1) {
     2845      /* comparand = (((((uj * B) + uj1) % v1) * b) + uj2);
     2846         guess = (((uj * B) + uj1) / v1); */
     2847      cl = u_scan[-2];
     2848      ch = bignum_divide_digit(uj, (u_scan[-1]), v1, (&gm));
     2849      guess = gm;
     2850    } else {
     2851      cl = u_scan[-2];
     2852      ch = u_scan[-1] + v1;
     2853      guess = C_BIGNUM_DIGIT_MASK;
     2854    }
     2855    while (1) {
     2856      /* product = (guess * v2); */
     2857      gl = C_BIGNUM_DIGIT_LO_HALF(guess);
     2858      gh = C_BIGNUM_DIGIT_HI_HALF(guess);
     2859      pl = v2l * gl;
     2860      ph = v2l * gh + v2h * gl + C_BIGNUM_DIGIT_HI_HALF(pl);
     2861      pl = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(ph),
     2862                                  C_BIGNUM_DIGIT_LO_HALF(pl));
     2863      ph = v2h * gh + C_BIGNUM_DIGIT_HI_HALF(ph);
     2864      /* if (comparand >= product) */
     2865      if ((ch > ph) || ((ch == ph) && (cl >= pl)))
     2866        break;
     2867      guess--;
     2868      /* comparand += (v1 << C_BIGNUM_DIGIT_LENGTH) */
     2869      ch += v1;
     2870      /* if (comparand >= (B^2)) */
     2871      if (!C_fitsinbignumdigitp(ch))
     2872        break;
     2873    }
     2874    qj = bignum_divide_and_subtract(v_start, v_end, guess, (--u_scan_start));
     2875    if (q_scan != NULL)
     2876      (*--q_scan) = qj;
     2877  }
     2878}
     2879
     2880static C_word
     2881bignum_divide_and_subtract(C_word *v_start, C_word *v_end, C_word guess, C_word *u_start)
     2882{
     2883  if (guess == 0) {
     2884    return 0;
     2885  } else {
     2886    C_word *v_scan = v_start,
     2887           *u_scan = u_start,
     2888           carry = 0,
     2889           gl, gh, v, pl, vl, vh, ph, diff, sum;
     2890
     2891    gl = C_BIGNUM_DIGIT_LO_HALF(guess);
     2892    gh = C_BIGNUM_DIGIT_HI_HALF(guess);
     2893    while (v_scan < v_end) {
     2894      v = (*v_scan++);
     2895      vl = C_BIGNUM_DIGIT_LO_HALF(v);
     2896      vh = C_BIGNUM_DIGIT_HI_HALF(v);
     2897      pl = vl * gl + C_BIGNUM_DIGIT_LO_HALF(carry);
     2898      ph = vl * gh + vh * gl + C_BIGNUM_DIGIT_HI_HALF(pl) +
     2899                               C_BIGNUM_DIGIT_HI_HALF(carry);
     2900      diff = (*u_scan) - C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(ph),
     2901                                                C_BIGNUM_DIGIT_LO_HALF(pl));
     2902      if (diff < 0) {
     2903        (*u_scan++) = diff + ((C_word)1 << C_BIGNUM_DIGIT_LENGTH);
     2904        carry = vh * gh + C_BIGNUM_DIGIT_HI_HALF(ph) + 1;
     2905      } else {
     2906        (*u_scan++) = diff;
     2907        carry = vh * gh + C_BIGNUM_DIGIT_HI_HALF(ph);
     2908      }
     2909    }
     2910    if (carry == 0)
     2911      return guess;
     2912
     2913    diff = ((*u_scan) - carry);
     2914    if (diff < 0) {
     2915      (*u_scan) = diff + ((C_word)1 << C_BIGNUM_DIGIT_LENGTH);
     2916    } else {
     2917      (*u_scan) = diff;
     2918      return guess;
     2919    }
     2920 
     2921    /* Subtraction generated carry, implying guess is one too large.
     2922       Add v back in to bring it back down. */
     2923    v_scan = v_start;
     2924    u_scan = u_start;
     2925    carry = 0;
     2926    while (v_scan < v_end) {
     2927      sum = ((*v_scan++) + (*u_scan) + carry);
     2928      if (sum < BIGNUM_RADIX) {
     2929        (*u_scan++) = sum;
     2930        carry = 0;
     2931      } else {
     2932        (*u_scan++) = sum & C_BIGNUM_DIGIT_MASK;
     2933        carry = 1;
     2934      }
     2935    }
     2936    if (carry) {
     2937      sum = (*u_scan) + carry;
     2938      (*u_scan) = C_fitsinbignumdigitp(sum) ? sum : (sum  & C_BIGNUM_DIGIT_MASK);
     2939    }
     2940    return guess - 1;
     2941  }
     2942}
     2943
     2944/* Like bignum_normalize, but this also shifts division-normalized numbers */
     2945static C_word
     2946bignum_normalize_shifted(C_word big, C_word shift_right)
     2947{
     2948  C_word length = C_bignum_size(big),
     2949        *start = C_bignum_digits(big),
     2950        *scan = start + length,
     2951        digit, carry = 0,
     2952        shift_left = C_BIGNUM_DIGIT_LENGTH - shift_right,
     2953        mask = (1L << shift_right) - 1;
     2954 
     2955  while (!(*--scan)) {
     2956    if (start == scan) { /* Don't bother with anything else */
     2957      return C_fix(0);
     2958    }
     2959    --length;
     2960  }
     2961
     2962  digit = (*scan);
     2963  (*scan) = (digit >> shift_right);
     2964  length -= (*scan == 0); /* Add 1 or 0 */
     2965  carry = ((digit & mask) << shift_left);
     2966 
     2967  while (start < scan) {
     2968    digit = (*--scan);
     2969    (*scan) = ((digit >> shift_right) | carry);
     2970    carry = ((digit & mask) << shift_left);
     2971  }
     2972  assert(carry == 0);
     2973  assert(C_bignum_size(big) != length);
     2974  assert(length != 1 || *C_bignum_digits(big) != 0);
     2975
     2976  switch (length) {
     2977  case 0:
     2978    return C_fix(0);
     2979  case 1:
     2980    return C_fix(C_bignum_negativep(big) ? -*start : *start);
     2981  case 2:
     2982    if (C_bignum_negativep(big) && *scan == 1 && *start == 0)
     2983      return C_fix(C_MOST_NEGATIVE_FIXNUM);
     2984    /* FALLTHROUGH */
     2985  default:
     2986    /* Mutate vector size of internal bignum vector. */
     2987    C_block_header(C_internal_bignum(big)) = (C_STRING_TYPE | C_wordstobytes(length+1));
     2988    /* Set internal header. */
     2989    C_bignum_header(big) = (C_bignum_header(big) & C_BIGNUM_HEADER_SIGN_BIT) | length;
     2990    return big;
     2991  }
     2992}
  • release/4/numbers/trunk/numbers-c.h

    r31392 r31411  
    196196
    197197void C_ccall C_u_bignum_quotient_fixnum(C_word c, C_word self, C_word k, C_word x, C_word y);
     198void C_ccall C_u_bignum_quotient_bignum(C_word c, C_word self, C_word k, C_word x, C_word y);
    198199
    199200void C_ccall C_u_fixnum_minus_bignum(C_word c, C_word self, C_word k, C_word x, C_word y);
  • release/4/numbers/trunk/numbers.scm

    r31390 r31411  
    152152
    153153(define %big-quotient-fix (##core#primitive "C_u_bignum_quotient_fixnum"))
    154 (define %big-quotient-big (##core#primitive "big_quotient_big"))
     154(define %big-quotient-big (##core#primitive "C_u_bignum_quotient_bignum"))
    155155
    156156(define %big-remainder-fix (##core#primitive "big_remainder_fix"))
     
    469469              (ratnum (fx/ x g) (fx/ y g)))]
    470470       [FLO (fp/ (%fix->flo x) y)]
     471       ;; XXX TODO: Get rid of the gcd call here, which also divides!
    471472       [BIG (let ((g (%gcd-0 '/ x y)))
    472473              (ratnum (%quotient '/ x g) (%quotient '/ y g)))]
Note: See TracChangeset for help on using the changeset viewer.