Changeset 31390 in project
 Timestamp:
 09/11/14 21:35:58 (7 years ago)
 Location:
 release/4/numbers/trunk
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

release/4/numbers/trunk/numbersc.c
r31376 r31390 1070 1070 1071 1071 static 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) == 11090 && BIGNUM_REF(bigx, 1) == 1 && BIGNUM_REF(bigx, 0) == 0) {1091 /*1092 * Very very special case:1093 * quotient(MOST_NEGATIVE_FIXNUM, (MOST_NEGATIVE_FIXNUM)) => 11094 */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_denominator1102 (bigx, bigy, ("ient), ((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_denominator1108 (bigx, abs_y, ("ient), ((bignum_type *) 0), neg_p, 0);1109 C_return_bignum(k, quotient);1110 } else {1111 bignum_divide_unsigned_medium_denominator1112 (bigx, abs_y, ("ient), ((bignum_type *) 0), neg_p, 0);1113 C_return_bignum(k, quotient);1114 }1115 }1116 }1117 1118 static void1119 1072 big_quotient_big(C_word c, C_word self, C_word k, C_word x, C_word y) 1120 1073 { … … 1247 1200 static void bignum_posneg_bitwise_op(C_word c, C_word self, C_word result); 1248 1201 static void bignum_pospos_bitwise_op(C_word c, C_word self, C_word result); 1202 static void bignum_destructive_normalize(C_word target, C_word source, C_word shift_left); 1203 static void bignum_destructive_divide_unsigned_halfdigit(C_word c, C_word self, C_word quotient); 1204 static void bignum_destructive_divide_unsigned_digit(C_word c, C_word self, C_word quotient); 1205 static C_word bignum_divide_digit(C_word uh, C_word ul, C_word v, C_word *q); 1206 static C_word bignum_divide_digit_subtract(C_word v1, C_word v2, C_word guess, C_word *u); 1249 1207 1250 1208 /* Eventually this will probably need to be integrated into C_2_plus. */ … … 1604 1562 } 1605 1563 1564 /* XXX TODO: Maybe pass true/false/bignum as initp, to allow copying data? */ 1606 1565 void C_ccall 1607 1566 C_allocate_bignum(C_word c, C_word self, C_word k, C_word size, C_word negp, C_word initp) … … 1695 1654 /* Given (denominator > 1), it is fairly easy to show that 1696 1655 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. */ 1698 1660 static C_word 1699 1661 bignum_digits_destructive_scale_down(C_word *start, C_word *end, C_word denominator) … … 1823 1785 ((R * R) + (R * (R  2)) + (R  1)) 1824 1786 Maximum value for carry: ((R * (R  1)) + (R  1)) 1825 where R == 2^HALF_DIGIT_ SHIFT*/1787 where R == 2^HALF_DIGIT_LENGTH */ 1826 1788 static void 1827 1789 bignum_times_bignum_unsigned(C_word k, C_word x, C_word y, C_word negp) … … 2560 2522 C_kontinue(k, C_bignum_normalize(result)); 2561 2523 } 2524 2525 void C_ccall 2526 C_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), ("ient), (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 2567 static void 2568 bignum_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 */ 2602 static void 2603 bignum_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_LENGTH1))) { 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 *(scan1) = 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 2661 static void 2662 bignum_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 */ 2684 static C_word 2685 bignum_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 2755 static C_word 2756 bignum_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/numbersc.h
r31379 r31390 149 149 150 150 #ifdef C_SIXTY_FOUR 151 # define C_BIGNUM_DIGIT_LENGTH 62152 # define C_BIGNUM_HEADER_SIGN_BIT 0x4000000000000000L153 # define C_BIGNUM_HEADER_SIZE_MASK 0x3fffffffffffffffL154 # define C_BIGNUM_DIGIT_MASK 0x3fffffffffffffffL155 # define C_BIGNUM_HALF_DIGIT_MASK 0x000000007fffffffL156 # define C_BIGNUM_HALF_DIGIT_ SHIFT31151 # 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 157 157 #else 158 # define C_BIGNUM_DIGIT_LENGTH 30159 # define C_BIGNUM_HEADER_SIGN_BIT 0x40000000160 # define C_BIGNUM_HEADER_SIZE_MASK 0x3fffffff161 # define C_BIGNUM_DIGIT_MASK 0x3fffffff162 # define C_BIGNUM_HALF_DIGIT_MASK 0x00007fff163 # define C_BIGNUM_HALF_DIGIT_ SHIFT15158 # 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 164 164 #endif 165 165 … … 168 168 169 169 #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)) 172 172 173 173 #define C_fitsinbignumdigitp(n) ((C_uword)(n) == ((C_uword)(n) & C_BIGNUM_DIGIT_MASK)) … … 204 204 void C_ccall C_u_2_bignum_times(C_word c, C_word self, C_word k, C_word x, C_word y); 205 205 206 void C_ccall C_u_bignum_quotient_fixnum(C_word c, C_word self, C_word k, C_word x, C_word y); 207 206 208 void C_ccall C_u_fixnum_minus_bignum(C_word c, C_word self, C_word k, C_word x, C_word y); 207 209 void 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 151 151 (define %big*big (##core#primitive "C_u_2_bignum_times")) 152 152 153 (define %bigquotientfix (##core#primitive " big_quotient_fix"))153 (define %bigquotientfix (##core#primitive "C_u_bignum_quotient_fixnum")) 154 154 (define %bigquotientbig (##core#primitive "big_quotient_big")) 155 155
Note: See TracChangeset
for help on using the changeset viewer.