home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / mitsch75.zip / scheme-7_5_17-src.zip / scheme-7.5.17 / src / microcode / bignum.c < prev    next >
C/C++ Source or Header  |  2000-12-05  |  58KB  |  1,937 lines
  1. /* -*-C-*-
  2.  
  3. $Id: bignum.c,v 9.49 2000/12/05 21:23:43 cph Exp $
  4.  
  5. Copyright (c) 1989-2000 Massachusetts Institute of Technology
  6.  
  7. This program is free software; you can redistribute it and/or modify
  8. it under the terms of the GNU General Public License as published by
  9. the Free Software Foundation; either version 2 of the License, or (at
  10. your option) any later version.
  11.  
  12. This program is distributed in the hope that it will be useful, but
  13. WITHOUT ANY WARRANTY; without even the implied warranty of
  14. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15. General Public License for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with this program; if not, write to the Free Software
  19. Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20. */
  21.  
  22. /* Implementation of Bignums (unlimited precision integers) */
  23.  
  24. #ifdef MIT_SCHEME
  25. #include "scheme.h"
  26. #undef CHAR_BIT /* redefined in "limits.h" */
  27. #else
  28. #include "bignum.h"
  29. #endif
  30.  
  31. #include "bignmint.h"
  32. #include "limits.h"
  33. #include <math.h>
  34.  
  35. #ifndef MIT_SCHEME
  36.  
  37. static bignum_type
  38. DEFUN (bignum_malloc, (length), bignum_length_type length)
  39. {
  40.   extern char * malloc ();
  41.   char * result = (malloc ((length + 1) * (sizeof (bignum_digit_type))));
  42.   BIGNUM_ASSERT (result != ((char *) 0));
  43.   return ((bignum_type) result);
  44. }
  45.  
  46. static bignum_type
  47. DEFUN (bignum_realloc, (bignum, length),
  48.        bignum_type bignum AND bignum_length_type length)
  49. {
  50.   extern char * realloc ();
  51.   char * result =
  52.     (realloc (((char *) bignum),
  53.           ((length + 1) * (sizeof (bignum_digit_type)))));
  54.   BIGNUM_ASSERT (result != ((char *) 0));
  55.   return ((bignum_type) result);
  56. }
  57.  
  58. #endif /* not MIT_SCHEME */
  59.  
  60. /* Forward references */
  61. static int EXFUN (bignum_equal_p_unsigned,
  62.           (bignum_type, bignum_type));
  63. static enum bignum_comparison EXFUN (bignum_compare_unsigned,
  64.                      (bignum_type, bignum_type));
  65. static bignum_type EXFUN (bignum_add_unsigned,
  66.               (bignum_type, bignum_type, int));
  67. static bignum_type EXFUN (bignum_subtract_unsigned,
  68.               (bignum_type, bignum_type));
  69. static bignum_type EXFUN (bignum_multiply_unsigned,
  70.               (bignum_type, bignum_type, int));
  71. static bignum_type EXFUN (bignum_multiply_unsigned_small_factor,
  72.               (bignum_type, bignum_digit_type, int));
  73. static void EXFUN (bignum_destructive_scale_up,
  74.            (bignum_type, bignum_digit_type));
  75. static void EXFUN (bignum_destructive_add,
  76.            (bignum_type, bignum_digit_type));
  77. static void EXFUN (bignum_divide_unsigned_large_denominator,
  78.            (bignum_type, bignum_type, bignum_type *, bignum_type *,
  79.             int, int));
  80. static void EXFUN (bignum_destructive_normalization,
  81.            (bignum_type, bignum_type, int));
  82. static void EXFUN (bignum_destructive_unnormalization,
  83.            (bignum_type, int));
  84. static void EXFUN (bignum_divide_unsigned_normalized,
  85.            (bignum_type, bignum_type, bignum_type));
  86. static bignum_digit_type EXFUN (bignum_divide_subtract,
  87.                 (bignum_digit_type *, bignum_digit_type *,
  88.                  bignum_digit_type, bignum_digit_type *));
  89. static void EXFUN (bignum_divide_unsigned_medium_denominator,
  90.            (bignum_type, bignum_digit_type, bignum_type *,
  91.             bignum_type *, int, int));
  92. static bignum_digit_type EXFUN (bignum_digit_divide,
  93.                 (bignum_digit_type, bignum_digit_type,
  94.                  bignum_digit_type, bignum_digit_type *));
  95. static bignum_digit_type EXFUN (bignum_digit_divide_subtract,
  96.                 (bignum_digit_type, bignum_digit_type,
  97.                  bignum_digit_type, bignum_digit_type *));
  98. static void EXFUN (bignum_divide_unsigned_small_denominator,
  99.            (bignum_type, bignum_digit_type, bignum_type *,
  100.             bignum_type *, int, int));
  101. static bignum_digit_type EXFUN (bignum_destructive_scale_down,
  102.                 (bignum_type, bignum_digit_type));
  103. static bignum_type EXFUN (bignum_remainder_unsigned_small_denominator,
  104.               (bignum_type, bignum_digit_type, int));
  105. static bignum_type EXFUN (bignum_digit_to_bignum,
  106.               (bignum_digit_type, int));
  107. static bignum_type EXFUN (bignum_allocate,
  108.               (bignum_length_type, int));
  109. static bignum_type EXFUN (bignum_allocate_zeroed,
  110.               (bignum_length_type, int));
  111. static bignum_type EXFUN (bignum_shorten_length,
  112.               (bignum_type, bignum_length_type));
  113. static bignum_type EXFUN (bignum_trim,
  114.               (bignum_type));
  115. static bignum_type EXFUN (bignum_copy,
  116.               (bignum_type));
  117. static bignum_type EXFUN (bignum_new_sign,
  118.               (bignum_type, int));
  119. static bignum_type EXFUN (bignum_maybe_new_sign,
  120.               (bignum_type, int));
  121. static void EXFUN (bignum_destructive_copy,
  122.            (bignum_type, bignum_type));
  123.  
  124. #define ULONG_LENGTH_IN_BITS(digit, len)        \
  125. do {                            \
  126.   fast unsigned long w = digit;                \
  127.   len = 0;                        \
  128.   while (w > 0xff) { len += 8; w >>= 8; }        \
  129.   while (w > 0)    { len += 1; w >>= 1; }        \
  130. } while (0)
  131.  
  132. /* Exports */
  133.  
  134. bignum_type
  135. DEFUN_VOID (bignum_make_zero)
  136. {
  137.   fast bignum_type result = (BIGNUM_ALLOCATE (0));
  138.   BIGNUM_SET_HEADER (result, 0, 0);
  139.   return (result);
  140. }
  141.  
  142. bignum_type
  143. DEFUN (bignum_make_one, (negative_p), int negative_p)
  144. {
  145.   fast bignum_type result = (BIGNUM_ALLOCATE (1));
  146.   BIGNUM_SET_HEADER (result, 1, negative_p);
  147.   (BIGNUM_REF (result, 0)) = 1;
  148.   return (result);
  149. }
  150.  
  151. int
  152. DEFUN (bignum_equal_p, (x, y),
  153.        fast bignum_type x AND fast bignum_type y)
  154. {
  155.   return
  156.     ((BIGNUM_ZERO_P (x))
  157.      ? (BIGNUM_ZERO_P (y))
  158.      : ((! (BIGNUM_ZERO_P (y)))
  159.     && ((BIGNUM_NEGATIVE_P (x))
  160.         ? (BIGNUM_NEGATIVE_P (y))
  161.         : (! (BIGNUM_NEGATIVE_P (y))))
  162.     && (bignum_equal_p_unsigned (x, y))));
  163. }
  164.  
  165. enum bignum_comparison
  166. DEFUN (bignum_test, (bignum), fast bignum_type bignum)
  167. {
  168.   return
  169.     ((BIGNUM_ZERO_P (bignum))
  170.      ? bignum_comparison_equal
  171.      : (BIGNUM_NEGATIVE_P (bignum))
  172.      ? bignum_comparison_less
  173.      : bignum_comparison_greater);
  174. }
  175.  
  176. enum bignum_comparison
  177. DEFUN (bignum_compare, (x, y),
  178.        fast bignum_type x AND fast bignum_type y)
  179. {
  180.   return
  181.     ((BIGNUM_ZERO_P (x))
  182.      ? ((BIGNUM_ZERO_P (y))
  183.     ? bignum_comparison_equal
  184.     : (BIGNUM_NEGATIVE_P (y))
  185.     ? bignum_comparison_greater
  186.     : bignum_comparison_less)
  187.      : (BIGNUM_ZERO_P (y))
  188.      ? ((BIGNUM_NEGATIVE_P (x))
  189.     ? bignum_comparison_less
  190.     : bignum_comparison_greater)
  191.      : (BIGNUM_NEGATIVE_P (x))
  192.      ? ((BIGNUM_NEGATIVE_P (y))
  193.     ? (bignum_compare_unsigned (y, x))
  194.     : (bignum_comparison_less))
  195.      : ((BIGNUM_NEGATIVE_P (y))
  196.     ? (bignum_comparison_greater)
  197.     : (bignum_compare_unsigned (x, y))));
  198. }
  199.  
  200. bignum_type
  201. DEFUN (bignum_add, (x, y),
  202.        fast bignum_type x AND fast bignum_type y)
  203. {
  204.   return
  205.     ((BIGNUM_ZERO_P (x))
  206.      ? (BIGNUM_MAYBE_COPY (y))
  207.      : (BIGNUM_ZERO_P (y))
  208.      ? (BIGNUM_MAYBE_COPY (x))
  209.      : ((BIGNUM_NEGATIVE_P (x))
  210.     ? ((BIGNUM_NEGATIVE_P (y))
  211.        ? (bignum_add_unsigned (x, y, 1))
  212.        : (bignum_subtract_unsigned (y, x)))
  213.     : ((BIGNUM_NEGATIVE_P (y))
  214.        ? (bignum_subtract_unsigned (x, y))
  215.        : (bignum_add_unsigned (x, y, 0)))));
  216. }
  217.  
  218. bignum_type
  219. DEFUN (bignum_subtract, (x, y),
  220.        fast bignum_type x AND fast bignum_type y)
  221. {
  222.   return
  223.     ((BIGNUM_ZERO_P (x))
  224.      ? ((BIGNUM_ZERO_P (y))
  225.     ? (BIGNUM_MAYBE_COPY (y))
  226.     : (bignum_new_sign (y, (! (BIGNUM_NEGATIVE_P (y))))))
  227.      : ((BIGNUM_ZERO_P (y))
  228.     ? (BIGNUM_MAYBE_COPY (x))
  229.     : ((BIGNUM_NEGATIVE_P (x))
  230.        ? ((BIGNUM_NEGATIVE_P (y))
  231.           ? (bignum_subtract_unsigned (y, x))
  232.           : (bignum_add_unsigned (x, y, 1)))
  233.        : ((BIGNUM_NEGATIVE_P (y))
  234.           ? (bignum_add_unsigned (x, y, 0))
  235.           : (bignum_subtract_unsigned (x, y))))));
  236. }
  237.  
  238. bignum_type
  239. DEFUN (bignum_negate, (x), fast bignum_type x)
  240. {
  241.   return
  242.     ((BIGNUM_ZERO_P (x))
  243.      ? (BIGNUM_MAYBE_COPY (x))
  244.      : (bignum_new_sign (x, (! (BIGNUM_NEGATIVE_P (x))))));
  245. }
  246.  
  247. bignum_type
  248. DEFUN (bignum_multiply, (x, y),
  249.        fast bignum_type x AND fast bignum_type y)
  250. {
  251.   fast bignum_length_type x_length = (BIGNUM_LENGTH (x));
  252.   fast bignum_length_type y_length = (BIGNUM_LENGTH (y));
  253.   fast int negative_p =
  254.     ((BIGNUM_NEGATIVE_P (x))
  255.      ? (! (BIGNUM_NEGATIVE_P (y)))
  256.      : (BIGNUM_NEGATIVE_P (y)));
  257.   if (BIGNUM_ZERO_P (x))
  258.     return (BIGNUM_MAYBE_COPY (x));
  259.   if (BIGNUM_ZERO_P (y))
  260.     return (BIGNUM_MAYBE_COPY (y));
  261.   if (x_length == 1)
  262.     {
  263.       bignum_digit_type digit = (BIGNUM_REF (x, 0));
  264.       if (digit == 1)
  265.     return (bignum_maybe_new_sign (y, negative_p));
  266.       if (digit < BIGNUM_RADIX_ROOT)
  267.     return (bignum_multiply_unsigned_small_factor (y, digit, negative_p));
  268.     }
  269.   if (y_length == 1)
  270.     {
  271.       bignum_digit_type digit = (BIGNUM_REF (y, 0));
  272.       if (digit == 1)
  273.     return (bignum_maybe_new_sign (x, negative_p));
  274.       if (digit < BIGNUM_RADIX_ROOT)
  275.     return (bignum_multiply_unsigned_small_factor (x, digit, negative_p));
  276.     }
  277.   return (bignum_multiply_unsigned (x, y, negative_p));
  278. }
  279.  
  280. int
  281. DEFUN (bignum_divide, (numerator, denominator, quotient, remainder),
  282.        bignum_type numerator AND bignum_type denominator
  283.        AND bignum_type * quotient AND bignum_type * remainder)
  284. {
  285.   if (BIGNUM_ZERO_P (denominator))
  286.     return (1);
  287.   if (BIGNUM_ZERO_P (numerator))
  288.     {
  289.       (*quotient) = (BIGNUM_MAYBE_COPY (numerator));
  290.       (*remainder) = (BIGNUM_MAYBE_COPY (numerator));
  291.     }
  292.   else
  293.     {
  294.       int r_negative_p = (BIGNUM_NEGATIVE_P (numerator));
  295.       int q_negative_p =
  296.     ((BIGNUM_NEGATIVE_P (denominator)) ? (! r_negative_p) : r_negative_p);
  297.       switch (bignum_compare_unsigned (numerator, denominator))
  298.     {
  299.     case bignum_comparison_equal:
  300.       {
  301.         (*quotient) = (BIGNUM_ONE (q_negative_p));
  302.         (*remainder) = (BIGNUM_ZERO ());
  303.         break;
  304.       }
  305.     case bignum_comparison_less:
  306.       {
  307.         (*quotient) = (BIGNUM_ZERO ());
  308.         (*remainder) = (BIGNUM_MAYBE_COPY (numerator));
  309.         break;
  310.       }
  311.     case bignum_comparison_greater:
  312.       {
  313.         if ((BIGNUM_LENGTH (denominator)) == 1)
  314.           {
  315.         bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
  316.         if (digit == 1)
  317.           {
  318.             (*quotient) =
  319.               (bignum_maybe_new_sign (numerator, q_negative_p));
  320.             (*remainder) = (BIGNUM_ZERO ());
  321.             break;
  322.           }
  323.         else if (digit < BIGNUM_RADIX_ROOT)
  324.           {
  325.             bignum_divide_unsigned_small_denominator
  326.               (numerator, digit,
  327.                quotient, remainder,
  328.                q_negative_p, r_negative_p);
  329.             break;
  330.           }
  331.         else
  332.           {
  333.             bignum_divide_unsigned_medium_denominator
  334.               (numerator, digit,
  335.                quotient, remainder,
  336.                q_negative_p, r_negative_p);
  337.             break;
  338.           }
  339.           }
  340.         bignum_divide_unsigned_large_denominator
  341.           (numerator, denominator,
  342.            quotient, remainder,
  343.            q_negative_p, r_negative_p);
  344.         break;
  345.       }
  346.     }
  347.     }
  348.   return (0);
  349. }
  350.  
  351. bignum_type
  352. DEFUN (bignum_quotient, (numerator, denominator),
  353.        bignum_type numerator AND bignum_type denominator)
  354. {
  355.   if (BIGNUM_ZERO_P (denominator))
  356.     return (BIGNUM_OUT_OF_BAND);
  357.   if (BIGNUM_ZERO_P (numerator))
  358.     return (BIGNUM_MAYBE_COPY (numerator));
  359.   {
  360.     int q_negative_p =
  361.       ((BIGNUM_NEGATIVE_P (denominator))
  362.        ? (! (BIGNUM_NEGATIVE_P (numerator)))
  363.        : (BIGNUM_NEGATIVE_P (numerator)));
  364.     switch (bignum_compare_unsigned (numerator, denominator))
  365.       {
  366.       case bignum_comparison_equal:
  367.     return (BIGNUM_ONE (q_negative_p));
  368.       case bignum_comparison_less:
  369.     return (BIGNUM_ZERO ());
  370.       case bignum_comparison_greater:
  371.     {
  372.       bignum_type quotient;
  373.       if ((BIGNUM_LENGTH (denominator)) == 1)
  374.         {
  375.           bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
  376.           if (digit == 1)
  377.         return (bignum_maybe_new_sign (numerator, q_negative_p));
  378.           if (digit < BIGNUM_RADIX_ROOT)
  379.         bignum_divide_unsigned_small_denominator
  380.           (numerator, digit,
  381.            ("ient), ((bignum_type *) 0),
  382.            q_negative_p, 0);
  383.           else
  384.         bignum_divide_unsigned_medium_denominator
  385.           (numerator, digit,
  386.            ("ient), ((bignum_type *) 0),
  387.            q_negative_p, 0);
  388.         }
  389.       else
  390.         bignum_divide_unsigned_large_denominator
  391.           (numerator, denominator,
  392.            ("ient), ((bignum_type *) 0),
  393.            q_negative_p, 0);
  394.       return (quotient);
  395.     }
  396.       default:
  397.     /*NOTREACHED*/
  398.     return (0);
  399.       }
  400.   }
  401. }
  402.  
  403. bignum_type
  404. DEFUN (bignum_remainder, (numerator, denominator),
  405.        bignum_type numerator AND bignum_type denominator)
  406. {
  407.   if (BIGNUM_ZERO_P (denominator))
  408.     return (BIGNUM_OUT_OF_BAND);
  409.   if (BIGNUM_ZERO_P (numerator))
  410.     return (BIGNUM_MAYBE_COPY (numerator));
  411.   switch (bignum_compare_unsigned (numerator, denominator))
  412.     {
  413.     case bignum_comparison_equal:
  414.       return (BIGNUM_ZERO ());
  415.     case bignum_comparison_less:
  416.       return (BIGNUM_MAYBE_COPY (numerator));
  417.     case bignum_comparison_greater:
  418.       {
  419.     bignum_type remainder;
  420.     if ((BIGNUM_LENGTH (denominator)) == 1)
  421.       {
  422.         bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
  423.         if ((digit & (digit-1)) == 0) /* i.e. digit = 2^k, including 1 */
  424.           {
  425.         bignum_digit_type unsigned_remainder;
  426.         if (BIGNUM_LENGTH (numerator) == 0)
  427.           return (BIGNUM_ZERO ());
  428.         unsigned_remainder = (digit-1) & (BIGNUM_REF (numerator, 0));
  429.         if (unsigned_remainder == 0)
  430.           return (BIGNUM_ZERO ());
  431.         return (bignum_digit_to_bignum
  432.             (unsigned_remainder, BIGNUM_NEGATIVE_P (numerator)));
  433.           }
  434.         if (digit < BIGNUM_RADIX_ROOT)
  435.           return
  436.         (bignum_remainder_unsigned_small_denominator
  437.          (numerator, digit, (BIGNUM_NEGATIVE_P (numerator))));
  438.         bignum_divide_unsigned_medium_denominator
  439.           (numerator, digit,
  440.            ((bignum_type *) 0), (&remainder),
  441.            0, (BIGNUM_NEGATIVE_P (numerator)));
  442.       }
  443.     else
  444.       bignum_divide_unsigned_large_denominator
  445.         (numerator, denominator,
  446.          ((bignum_type *) 0), (&remainder),
  447.          0, (BIGNUM_NEGATIVE_P (numerator)));
  448.     return (remainder);
  449.       }
  450.     default:
  451.       /*NOTREACHED*/
  452.       return (0);
  453.     }
  454. }
  455.  
  456. /* These procedures depend on the non-portable type `unsigned long'.
  457.    If your compiler doesn't support this type, either define the
  458.    switch `BIGNUM_NO_ULONG' to disable them (in "bignum.h"), or write
  459.    new versions that don't use this type. */
  460.  
  461. #ifndef BIGNUM_NO_ULONG
  462.  
  463. bignum_type
  464. DEFUN (long_to_bignum, (n), long n)
  465. {
  466.   int negative_p;
  467.   bignum_digit_type result_digits [BIGNUM_DIGITS_FOR_LONG];
  468.   fast bignum_digit_type * end_digits = result_digits;
  469.   /* Special cases win when these small constants are cached. */
  470.   if (n == 0) return (BIGNUM_ZERO ());
  471.   if (n == 1) return (BIGNUM_ONE (0));
  472.   if (n == -1) return (BIGNUM_ONE (1));
  473.   {
  474.     fast unsigned long accumulator = ((negative_p = (n < 0)) ? (-n) : n);
  475.     do
  476.       {
  477.     (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK);
  478.     accumulator >>= BIGNUM_DIGIT_LENGTH;
  479.       }
  480.     while (accumulator != 0);
  481.   }
  482.   {
  483.     bignum_type result =
  484.       (bignum_allocate ((end_digits - result_digits), negative_p));
  485.     fast bignum_digit_type * scan_digits = result_digits;
  486.     fast bignum_digit_type * scan_result = (BIGNUM_START_PTR (result));
  487.     while (scan_digits < end_digits)
  488.       (*scan_result++) = (*scan_digits++);
  489.     return (result);
  490.   }
  491. }
  492.  
  493. long
  494. DEFUN (bignum_to_long, (bignum), bignum_type bignum)
  495. {
  496.   if (BIGNUM_ZERO_P (bignum))
  497.     return (0);
  498.   {
  499.     fast unsigned long accumulator = 0;
  500.     fast bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  501.     fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
  502.     while (start < scan)
  503.       accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan));
  504.     return
  505.       ((BIGNUM_NEGATIVE_P (bignum))
  506.        ? (- ((long) accumulator))
  507.        : ((long) accumulator));
  508.   }
  509. }
  510.  
  511. bignum_type
  512. DEFUN (ulong_to_bignum, (n), unsigned long n)
  513. {
  514.   bignum_digit_type result_digits [BIGNUM_DIGITS_FOR_LONG];
  515.   fast bignum_digit_type * end_digits = result_digits;
  516.   /* Special cases win when these small constants are cached. */
  517.   if (n == 0) return (BIGNUM_ZERO ());
  518.   if (n == 1) return (BIGNUM_ONE (0));
  519.   {
  520.     fast unsigned long accumulator = n;
  521.     do
  522.       {
  523.     (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK);
  524.     accumulator >>= BIGNUM_DIGIT_LENGTH;
  525.       }
  526.     while (accumulator != 0);
  527.   }
  528.   {
  529.     bignum_type result =
  530.       (bignum_allocate ((end_digits - result_digits), 0));
  531.     fast bignum_digit_type * scan_digits = result_digits;
  532.     fast bignum_digit_type * scan_result = (BIGNUM_START_PTR (result));
  533.     while (scan_digits < end_digits)
  534.       (*scan_result++) = (*scan_digits++);
  535.     return (result);
  536.   }
  537. }
  538.  
  539. unsigned long
  540. DEFUN (bignum_to_ulong, (bignum), bignum_type bignum)
  541. {
  542.   if (BIGNUM_ZERO_P (bignum))
  543.     return (0);
  544.   {
  545.     fast unsigned long accumulator = 0;
  546.     fast bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  547.     fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
  548.     while (start < scan)
  549.       accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan));
  550.     return (accumulator);
  551.   }
  552. }
  553.  
  554. #endif /* not BIGNUM_NO_ULONG */
  555.  
  556. #define DTB_WRITE_DIGIT(factor)                        \
  557. {                                    \
  558.   significand *= (factor);                        \
  559.   digit = ((bignum_digit_type) significand);                \
  560.   (*--scan) = digit;                            \
  561.   significand -= ((double) digit);                    \
  562. }
  563.  
  564. bignum_type
  565. DEFUN (double_to_bignum, (x), double x)
  566. {
  567.   extern double frexp ();
  568.   int exponent;
  569.   fast double significand = (frexp (x, (&exponent)));
  570.   if (exponent <= 0) return (BIGNUM_ZERO ());
  571.   if (exponent == 1) return (BIGNUM_ONE (x < 0));
  572.   if (significand < 0) significand = (-significand);
  573.   {
  574.     bignum_length_type length = (BIGNUM_BITS_TO_DIGITS (exponent));
  575.     bignum_type result = (bignum_allocate (length, (x < 0)));
  576.     bignum_digit_type * start = (BIGNUM_START_PTR (result));
  577.     fast bignum_digit_type * scan = (start + length);
  578.     fast bignum_digit_type digit;
  579.     int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH);
  580.     if (odd_bits > 0)
  581.       DTB_WRITE_DIGIT (1L << odd_bits);
  582.     while (start < scan)
  583.       {
  584.     if (significand == 0)
  585.       {
  586.         while (start < scan)
  587.           (*--scan) = 0;
  588.         break;
  589.       }
  590.     DTB_WRITE_DIGIT (BIGNUM_RADIX);
  591.       }
  592.     return (result);
  593.   }
  594. }
  595.  
  596. #undef DTB_WRITE_DIGIT
  597.  
  598. /*  This version sometimes gets the answer wrong due to rounding errors.
  599.     It would be slightly better if the digits were accumulated lsb to msb.
  600.  */
  601. /*
  602. double
  603. DEFUN (bignum_to_double, (bignum), bignum_type bignum)
  604. {
  605.   if (BIGNUM_ZERO_P (bignum))
  606.     return (0);
  607.   {
  608.     fast double accumulator = 0;
  609.     fast bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  610.     fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
  611.     while (start < scan)
  612.       accumulator = ((accumulator * BIGNUM_RADIX) + (*--scan));
  613.     return ((BIGNUM_NEGATIVE_P (bignum)) ? (-accumulator) : accumulator);
  614.   }
  615. }
  616. */
  617.  
  618. double
  619. DEFUN (bignum_to_double, (bignum), bignum_type bignum)
  620. {
  621.   if (BIGNUM_ZERO_P (bignum))
  622.     return (0.0);
  623.  
  624.   {
  625.     bignum_length_type length = (BIGNUM_LENGTH (bignum));
  626.     bignum_length_type index = length - 1;
  627.     bignum_length_type scale_words = length - 1;
  628.     bignum_digit_type msd = (BIGNUM_REF (bignum, (index)));
  629. #if (FLT_RADIX == 2)
  630.     int bits_to_get = DBL_MANT_DIG; /* includes implicit 1 */
  631. #else
  632. #include "error: must have FLT_RADIX==2"
  633. #endif
  634.     double value = 0;
  635.     bignum_digit_type mask = 0;
  636.     bignum_digit_type guard_bit_mask = BIGNUM_RADIX>>1;
  637.     bignum_digit_type rounding_correction = 0;
  638.     int current_digit_bit_count = 0;
  639.  
  640.     ULONG_LENGTH_IN_BITS (msd, current_digit_bit_count);
  641.     mask = (1 << (current_digit_bit_count)) - 1;
  642.  
  643.     while (1) {
  644.       if (current_digit_bit_count > bits_to_get) {
  645.     guard_bit_mask = (1 << (current_digit_bit_count - bits_to_get - 1));
  646.     mask &= ~((guard_bit_mask << 1) - 1);
  647.     current_digit_bit_count = bits_to_get;
  648.     bits_to_get = 0;
  649.       } else {
  650.     bits_to_get -= current_digit_bit_count;
  651.       }
  652.  
  653.       value = (value * BIGNUM_RADIX) + ((BIGNUM_REF (bignum, index)) & mask);
  654.  
  655.       if (bits_to_get == 0) {
  656.     scale_words = index;
  657.     if (current_digit_bit_count == BIGNUM_DIGIT_LENGTH) {
  658.       if (index == 0) /* there is no guard bit */
  659.         goto finished;
  660.       guard_bit_mask = (1 << (BIGNUM_DIGIT_LENGTH - 1));
  661.       rounding_correction = 1;
  662.       index -= 1;
  663.     } else {
  664.       rounding_correction = (guard_bit_mask << 1);
  665.     }
  666.     break;
  667.       }
  668.       if (index == 0)  /* fewer than DBL_MANT_DIG bits */
  669.     goto finished;
  670.  
  671.       index -= 1;
  672.       current_digit_bit_count = BIGNUM_DIGIT_LENGTH;
  673.       mask = BIGNUM_DIGIT_MASK;
  674.     }
  675.  
  676.     /* round-to-even depending on lsb, guard and following bits: lgfffff */
  677.  
  678.     if ((BIGNUM_REF(bignum,index) & guard_bit_mask) == 0) /* case x0xxxx */
  679.       goto round_down;
  680.  
  681.     if ((BIGNUM_REF(bignum,index) & (guard_bit_mask-1)) != 0) /* case x1xx1x */
  682.       goto round_up;
  683.  
  684.     /* cases 110000 and 1101xx: test "odd?", i.e. round-to-even rounds up */
  685.     if ((guard_bit_mask << 1) == BIGNUM_RADIX) {
  686.       if (((BIGNUM_REF (bignum, index+1)) & 1) != 0)  /* "odd?" */
  687.     goto round_up;
  688.     } else {
  689.       if (((BIGNUM_REF (bignum, index)) & (guard_bit_mask << 1)) != 0)
  690.     goto round_up;
  691.     }
  692.  
  693.     if (index==0)   /* case 010000, no more words of following bits */
  694.       goto finished;
  695.  
  696.     { /* distinguish between cases 0100...00 and 0100..1xx, multiple words */
  697.       bignum_length_type index2 = index - 1;
  698.       while (index2 >= 0) {
  699.     if ((BIGNUM_REF (bignum, index2)) != 0)
  700.       goto round_up;
  701.     index2--;
  702.       }
  703.       goto round_down;
  704.     }
  705.  
  706.   round_up:
  707.     value += rounding_correction;
  708.   round_down:
  709.     /* note, ldexp `sticks' at the maximal non-infinity value, which
  710.        is a reasonable compromise for numbers with DBL_MAX_EXP bits
  711.        that round up */
  712.     if (scale_words > 0)
  713.       value = ldexp (value, scale_words * BIGNUM_DIGIT_LENGTH);
  714.  
  715.   finished:
  716.     return ((BIGNUM_NEGATIVE_P (bignum)) ? (-value) : value);
  717.   }
  718. }
  719.  
  720. /*
  721. ;; Code to test bignum_to_double
  722.  
  723. (declare (usual-integrations))
  724.  
  725. (define integer->flonum (make-primitive-procedure 'integer->flonum 2))
  726.  
  727. (define (check n)
  728.   (let ((f1 (integer->flonum n #b10))
  729.     (f2 (exact->inexact n)))
  730.     (if (and f1 (= f1 f2))
  731.     (number->string f1 2)
  732.     (begin
  733.       (pp n)
  734.       (pp (number->string f1 2))
  735.       (pp (number->string f2 2))))))
  736.  
  737. (define (test)
  738.   (define n 0)
  739.   (do ((i 0 (+ i 1)))            ; guard bit zone patterns
  740.       ((= i 256))
  741.     (do ((j 0 (+ j 1)))            ; general word alignment
  742.     ((= j 30))
  743.       (do ((e 0 (+ e 1)))        ; random `insignificant' bit
  744.       ((= e 2))
  745.     (do ((up 0 (+ up 1)))        ; test potential for carry
  746.         ((= up 2))
  747.       (do ((end-pad 0 (+ end-pad 100)))
  748.           ((> end-pad 100))
  749.         (set! n (+ n 1)) (if (= (remainder n 1000) 0) (pp `(test ,n)))
  750.         (let ((s (string-append "1"
  751.                    (make-string 48 (vector-ref '#(#\0 #\1) up))
  752.                    (string-pad-left (number->string i 2) 8 #\0)
  753.                     (make-string (* j 23) #\0) ; gcd(23,30)=1
  754.                     (number->string e)
  755.                     (make-string end-pad #\0))))
  756.           (check (string->number s 2)))))))))
  757. */
  758.  
  759. int
  760. DEFUN (bignum_fits_in_word_p, (bignum, word_length, twos_complement_p),
  761.        bignum_type bignum AND long word_length AND int twos_complement_p)
  762. {
  763.   unsigned int n_bits = (twos_complement_p ? (word_length - 1) : word_length);
  764.   BIGNUM_ASSERT (n_bits > 0);
  765.   {
  766.     fast bignum_length_type length = (BIGNUM_LENGTH (bignum));
  767.     fast unsigned int max_digits = (BIGNUM_BITS_TO_DIGITS (n_bits));
  768.     if (((unsigned int) length) < max_digits)
  769.       return (1);
  770.     if (((unsigned int) length) > max_digits)
  771.       return (0);
  772.     {
  773.       bignum_digit_type msd = (BIGNUM_REF (bignum, (length - 1)));
  774.       bignum_digit_type max
  775.     = (1L << (n_bits - ((length - 1) * BIGNUM_DIGIT_LENGTH)));
  776.       return
  777.     (((msd < max)
  778.       || (twos_complement_p
  779.           && (msd == max)
  780.           && (BIGNUM_NEGATIVE_P (bignum)))));
  781.     }
  782.   }
  783. }
  784.  
  785. bignum_type
  786. DEFUN (bignum_length_in_bits, (bignum), bignum_type bignum)
  787. {
  788.   if (BIGNUM_ZERO_P (bignum))
  789.     return (BIGNUM_ZERO ());
  790.   {
  791.     bignum_length_type index = ((BIGNUM_LENGTH (bignum)) - 1);
  792.     fast bignum_digit_type digit = (BIGNUM_REF (bignum, index));
  793.     fast bignum_digit_type delta = 0;
  794.     fast bignum_type result = (bignum_allocate (2, 0));
  795.     (BIGNUM_REF (result, 0)) = index;
  796.     (BIGNUM_REF (result, 1)) = 0;
  797.     bignum_destructive_scale_up (result, BIGNUM_DIGIT_LENGTH);
  798.     ULONG_LENGTH_IN_BITS (digit, delta);
  799.     bignum_destructive_add (result, ((bignum_digit_type) delta));
  800.     return (bignum_trim (result));
  801.   }
  802. }
  803.  
  804. bignum_type
  805. DEFUN_VOID (bignum_length_upper_limit)
  806. {
  807.   fast bignum_type result = (bignum_allocate (2, 0));
  808.   (BIGNUM_REF (result, 0)) = 0;
  809.   (BIGNUM_REF (result, 1)) = BIGNUM_DIGIT_LENGTH;
  810.   return (result);
  811. }
  812.  
  813. bignum_type
  814. DEFUN (bignum_shift_left, (n, m), bignum_type n AND unsigned long m)
  815. {
  816.   unsigned long ln = (BIGNUM_LENGTH (n));
  817.   unsigned long delta = 0;
  818.   if (m == 0)
  819.     return (n);
  820.  
  821.   ULONG_LENGTH_IN_BITS ((BIGNUM_REF (n, (ln - 1))), delta);
  822.  
  823.   {
  824.     unsigned long zeroes = (m / BIGNUM_DIGIT_LENGTH);
  825.     unsigned long shift = (m % BIGNUM_DIGIT_LENGTH);
  826.     unsigned long ln2
  827.       = (((ln - 1) + ((delta + m) / BIGNUM_DIGIT_LENGTH))
  828.      + (((delta + m) % BIGNUM_DIGIT_LENGTH) != 0));
  829.     bignum_type result = (bignum_allocate (ln2, (BIGNUM_NEGATIVE_P (n))));
  830.     bignum_digit_type * scan_n = (BIGNUM_START_PTR (n));
  831.     bignum_digit_type * end_n = (scan_n + ln);
  832.     bignum_digit_type * scan_result = (BIGNUM_START_PTR (result));
  833.     while ((zeroes--) > 0)
  834.       (*scan_result++) = 0;
  835.     if (shift == 0)
  836.       while (scan_n < end_n)
  837.     (*scan_result++) = (*scan_n++);
  838.     else
  839.       {
  840.     unsigned long temp = 0;
  841.     while (scan_n < end_n)
  842.       {
  843.         bignum_digit_type digit = (*scan_n++);
  844.         (*scan_result++) = (((digit << shift) & BIGNUM_DIGIT_MASK) | temp);
  845.         temp = (digit >> (BIGNUM_DIGIT_LENGTH - shift));
  846.       }
  847.     if (temp != 0)
  848.       (*scan_result) = temp;
  849.       }
  850.     return (result);
  851.   }
  852. }
  853.  
  854. bignum_type
  855. DEFUN (unsigned_long_to_shifted_bignum, (n, m, sign),
  856.        unsigned long n AND
  857.        unsigned long m AND
  858.        int sign)
  859. {
  860.   unsigned long delta = 0;
  861.   if (n == 0)
  862.     return (BIGNUM_ZERO ());
  863.  
  864.   ULONG_LENGTH_IN_BITS (n, delta);
  865.  
  866.   {
  867.     unsigned long zeroes = (m / BIGNUM_DIGIT_LENGTH);
  868.     unsigned long shift = (m % BIGNUM_DIGIT_LENGTH);
  869.     unsigned long ln
  870.       = (((delta + m) / BIGNUM_DIGIT_LENGTH)
  871.      + (((delta + m) % BIGNUM_DIGIT_LENGTH) != 0));
  872.     bignum_type result = (bignum_allocate (ln, sign));
  873.     bignum_digit_type * scan_result = (BIGNUM_START_PTR (result));
  874.     while ((zeroes--) > 0)
  875.       (*scan_result++) = 0;
  876.     (*scan_result++) = ((n << shift) & BIGNUM_DIGIT_MASK);
  877.     n >>= (BIGNUM_DIGIT_LENGTH - shift);
  878.     while (n > 0)
  879.       {
  880.     (*scan_result++) = (n & BIGNUM_DIGIT_MASK);
  881.     n >>= BIGNUM_DIGIT_LENGTH;
  882.       }
  883.     return (result);
  884.   }
  885. }
  886.  
  887. bignum_type
  888. DEFUN (digit_stream_to_bignum,
  889.        (n_digits, producer, context, radix, negative_p),
  890.        fast unsigned int n_digits
  891.        AND unsigned int EXFUN ((*producer), (bignum_procedure_context))
  892.        AND bignum_procedure_context context
  893.        AND fast unsigned int radix
  894.        AND int negative_p)
  895. {
  896.   BIGNUM_ASSERT ((radix > 1) && (radix <= BIGNUM_RADIX_ROOT));
  897.   if (n_digits == 0)
  898.     return (BIGNUM_ZERO ());
  899.   if (n_digits == 1)
  900.     {
  901.       long digit = ((long) ((*producer) (context)));
  902.       return (long_to_bignum (negative_p ? (- digit) : digit));
  903.     }
  904.   {
  905.     bignum_length_type length;
  906.     {
  907.       fast unsigned int log_radix = 0;
  908.       ULONG_LENGTH_IN_BITS (radix, log_radix);
  909.       /* This length will be at least as large as needed. */
  910.       length = (BIGNUM_BITS_TO_DIGITS (n_digits * log_radix));
  911.     }
  912.     {
  913.       fast bignum_type result = (bignum_allocate_zeroed (length, negative_p));
  914.       while ((n_digits--) > 0)
  915.     {
  916.       bignum_destructive_scale_up (result, ((bignum_digit_type) radix));
  917.       bignum_destructive_add
  918.         (result, ((bignum_digit_type) ((*producer) (context))));
  919.     }
  920.       return (bignum_trim (result));
  921.     }
  922.   }
  923. }
  924.  
  925. void
  926. DEFUN (bignum_to_digit_stream, (bignum, radix, consumer, context),
  927.        bignum_type bignum
  928.        AND unsigned int radix
  929.        AND void EXFUN ((*consumer), (bignum_procedure_context, long))
  930.        AND bignum_procedure_context context)
  931. {
  932.   BIGNUM_ASSERT ((radix > 1) && (radix <= BIGNUM_RADIX_ROOT));
  933.   if (! (BIGNUM_ZERO_P (bignum)))
  934.     {
  935.       fast bignum_type working_copy = (bignum_copy (bignum));
  936.       fast bignum_digit_type * start = (BIGNUM_START_PTR (working_copy));
  937.       fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (working_copy)));
  938.       while (start < scan)
  939.     {
  940.       if ((scan[-1]) == 0)
  941.         scan -= 1;
  942.       else
  943.         (*consumer)
  944.           (context, (bignum_destructive_scale_down (working_copy, radix)));
  945.     }
  946.       BIGNUM_DEALLOCATE (working_copy);
  947.     }
  948.   return;
  949. }
  950.  
  951. long
  952. DEFUN_VOID (bignum_max_digit_stream_radix)
  953. {
  954.   return (BIGNUM_RADIX_ROOT);
  955. }
  956.  
  957. /* Comparisons */
  958.  
  959. static int
  960. DEFUN (bignum_equal_p_unsigned, (x, y),
  961.        bignum_type x AND bignum_type y)
  962. {
  963.   bignum_length_type length = (BIGNUM_LENGTH (x));
  964.   if (length != ((bignum_length_type) (BIGNUM_LENGTH (y))))
  965.     return (0);
  966.   else
  967.     {
  968.       fast bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
  969.       fast bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
  970.       fast bignum_digit_type * end_x = (scan_x + length);
  971.       while (scan_x < end_x)
  972.     if ((*scan_x++) != (*scan_y++))
  973.       return (0);
  974.       return (1);
  975.     }
  976. }
  977.  
  978. static enum bignum_comparison
  979. DEFUN (bignum_compare_unsigned, (x, y),
  980.        bignum_type x AND bignum_type y)
  981. {
  982.   bignum_length_type x_length = (BIGNUM_LENGTH (x));
  983.   bignum_length_type y_length = (BIGNUM_LENGTH (y));
  984.   if (x_length < y_length)
  985.     return (bignum_comparison_less);
  986.   if (x_length > y_length)
  987.     return (bignum_comparison_greater);
  988.   {
  989.     fast bignum_digit_type * start_x = (BIGNUM_START_PTR (x));
  990.     fast bignum_digit_type * scan_x = (start_x + x_length);
  991.     fast bignum_digit_type * scan_y = ((BIGNUM_START_PTR (y)) + y_length);
  992.     while (start_x < scan_x)
  993.       {
  994.     fast bignum_digit_type digit_x = (*--scan_x);
  995.     fast bignum_digit_type digit_y = (*--scan_y);
  996.     if (digit_x < digit_y)
  997.       return (bignum_comparison_less);
  998.     if (digit_x > digit_y)
  999.       return (bignum_comparison_greater);
  1000.       }
  1001.   }
  1002.   return (bignum_comparison_equal);
  1003. }
  1004.  
  1005. /* Addition */
  1006.  
  1007. static bignum_type
  1008. DEFUN (bignum_add_unsigned, (x, y, negative_p),
  1009.        bignum_type x AND bignum_type y AND int negative_p)
  1010. {
  1011.   if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
  1012.     {
  1013.       bignum_type z = x;
  1014.       x = y;
  1015.       y = z;
  1016.     }
  1017.   {
  1018.     bignum_length_type x_length = (BIGNUM_LENGTH (x));
  1019.     bignum_type r = (bignum_allocate ((x_length + 1), negative_p));
  1020.     fast bignum_digit_type sum;
  1021.     fast bignum_digit_type carry = 0;
  1022.     fast bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
  1023.     fast bignum_digit_type * scan_r = (BIGNUM_START_PTR (r));
  1024.     {
  1025.       fast bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
  1026.       fast bignum_digit_type * end_y = (scan_y + (BIGNUM_LENGTH (y)));
  1027.       while (scan_y < end_y)
  1028.     {
  1029.       sum = ((*scan_x++) + (*scan_y++) + carry);
  1030.       if (sum < BIGNUM_RADIX)
  1031.         {
  1032.           (*scan_r++) = sum;
  1033.           carry = 0;
  1034.         }
  1035.       else
  1036.         {
  1037.           (*scan_r++) = (sum - BIGNUM_RADIX);
  1038.           carry = 1;
  1039.         }
  1040.     }
  1041.     }
  1042.     {
  1043.       fast bignum_digit_type * end_x = ((BIGNUM_START_PTR (x)) + x_length);
  1044.       if (carry != 0)
  1045.     while (scan_x < end_x)
  1046.       {
  1047.         sum = ((*scan_x++) + 1);
  1048.         if (sum < BIGNUM_RADIX)
  1049.           {
  1050.         (*scan_r++) = sum;
  1051.         carry = 0;
  1052.         break;
  1053.           }
  1054.         else
  1055.           (*scan_r++) = (sum - BIGNUM_RADIX);
  1056.       }
  1057.       while (scan_x < end_x)
  1058.     (*scan_r++) = (*scan_x++);
  1059.     }
  1060.     if (carry != 0)
  1061.       {
  1062.     (*scan_r) = 1;
  1063.     return (r);
  1064.       }
  1065.     return (bignum_shorten_length (r, x_length));
  1066.   }
  1067. }
  1068.  
  1069. /* Subtraction */
  1070.  
  1071. static bignum_type
  1072. DEFUN (bignum_subtract_unsigned, (x, y),
  1073.        bignum_type x AND bignum_type y)
  1074. {
  1075.   int negative_p = 0;
  1076.   switch (bignum_compare_unsigned (x, y))
  1077.     {
  1078.     case bignum_comparison_equal:
  1079.       return (BIGNUM_ZERO ());
  1080.     case bignum_comparison_less:
  1081.       {
  1082.     bignum_type z = x;
  1083.     x = y;
  1084.     y = z;
  1085.       }
  1086.       negative_p = 1;
  1087.       break;
  1088.     case bignum_comparison_greater:
  1089.       negative_p = 0;
  1090.       break;
  1091.     }
  1092.   {
  1093.     bignum_length_type x_length = (BIGNUM_LENGTH (x));
  1094.     bignum_type r = (bignum_allocate (x_length, negative_p));
  1095.     fast bignum_digit_type difference;
  1096.     fast bignum_digit_type borrow = 0;
  1097.     fast bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
  1098.     fast bignum_digit_type * scan_r = (BIGNUM_START_PTR (r));
  1099.     {
  1100.       fast bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
  1101.       fast bignum_digit_type * end_y = (scan_y + (BIGNUM_LENGTH (y)));
  1102.       while (scan_y < end_y)
  1103.     {
  1104.       difference = (((*scan_x++) - (*scan_y++)) - borrow);
  1105.       if (difference < 0)
  1106.         {
  1107.           (*scan_r++) = (difference + BIGNUM_RADIX);
  1108.           borrow = 1;
  1109.         }
  1110.       else
  1111.         {
  1112.           (*scan_r++) = difference;
  1113.           borrow = 0;
  1114.         }
  1115.     }
  1116.     }
  1117.     {
  1118.       fast bignum_digit_type * end_x = ((BIGNUM_START_PTR (x)) + x_length);
  1119.       if (borrow != 0)
  1120.     while (scan_x < end_x)
  1121.       {
  1122.         difference = ((*scan_x++) - borrow);
  1123.         if (difference < 0)
  1124.           (*scan_r++) = (difference + BIGNUM_RADIX);
  1125.         else
  1126.           {
  1127.         (*scan_r++) = difference;
  1128.         borrow = 0;
  1129.         break;
  1130.           }
  1131.       }
  1132.       BIGNUM_ASSERT (borrow == 0);
  1133.       while (scan_x < end_x)
  1134.     (*scan_r++) = (*scan_x++);
  1135.     }
  1136.     return (bignum_trim (r));
  1137.   }
  1138. }
  1139.  
  1140. /* Multiplication
  1141.    Maximum value for product_low or product_high:
  1142.     ((R * R) + (R * (R - 2)) + (R - 1))
  1143.    Maximum value for carry: ((R * (R - 1)) + (R - 1))
  1144.     where R == BIGNUM_RADIX_ROOT */
  1145.  
  1146. static bignum_type
  1147. DEFUN (bignum_multiply_unsigned, (x, y, negative_p),
  1148.        bignum_type x AND bignum_type y AND int negative_p)
  1149. {
  1150.   if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
  1151.     {
  1152.       bignum_type z = x;
  1153.       x = y;
  1154.       y = z;
  1155.     }
  1156.   {
  1157.     fast bignum_digit_type carry;
  1158.     fast bignum_digit_type y_digit_low;
  1159.     fast bignum_digit_type y_digit_high;
  1160.     fast bignum_digit_type x_digit_low;
  1161.     fast bignum_digit_type x_digit_high;
  1162.     bignum_digit_type product_low;
  1163.     fast bignum_digit_type * scan_r;
  1164.     fast bignum_digit_type * scan_y;
  1165.     bignum_length_type x_length = (BIGNUM_LENGTH (x));
  1166.     bignum_length_type y_length = (BIGNUM_LENGTH (y));
  1167.     bignum_type r =
  1168.       (bignum_allocate_zeroed ((x_length + y_length), negative_p));
  1169.     bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
  1170.     bignum_digit_type * end_x = (scan_x + x_length);
  1171.     bignum_digit_type * start_y = (BIGNUM_START_PTR (y));
  1172.     bignum_digit_type * end_y = (start_y + y_length);
  1173.     bignum_digit_type * start_r = (BIGNUM_START_PTR (r));
  1174. #define x_digit x_digit_high
  1175. #define y_digit y_digit_high
  1176. #define product_high carry
  1177.     while (scan_x < end_x)
  1178.       {
  1179.     x_digit = (*scan_x++);
  1180.     x_digit_low = (HD_LOW (x_digit));
  1181.     x_digit_high = (HD_HIGH (x_digit));
  1182.     carry = 0;
  1183.     scan_y = start_y;
  1184.     scan_r = (start_r++);
  1185.     while (scan_y < end_y)
  1186.       {
  1187.         y_digit = (*scan_y++);
  1188.         y_digit_low = (HD_LOW (y_digit));
  1189.         y_digit_high = (HD_HIGH (y_digit));
  1190.         product_low =
  1191.           ((*scan_r) +
  1192.            (x_digit_low * y_digit_low) +
  1193.            (HD_LOW (carry)));
  1194.         product_high =
  1195.           ((x_digit_high * y_digit_low) +
  1196.            (x_digit_low * y_digit_high) +
  1197.            (HD_HIGH (product_low)) +
  1198.            (HD_HIGH (carry)));
  1199.         (*scan_r++) =
  1200.           (HD_CONS ((HD_LOW (product_high)), (HD_LOW (product_low))));
  1201.         carry =
  1202.           ((x_digit_high * y_digit_high) +
  1203.            (HD_HIGH (product_high)));
  1204.       }
  1205.     (*scan_r) += carry;
  1206.       }
  1207.     return (bignum_trim (r));
  1208. #undef x_digit
  1209. #undef y_digit
  1210. #undef product_high
  1211.   }
  1212. }
  1213.  
  1214. static bignum_type
  1215. DEFUN (bignum_multiply_unsigned_small_factor, (x, y, negative_p),
  1216.        bignum_type x AND bignum_digit_type y AND int negative_p)
  1217. {
  1218.   bignum_length_type length_x = (BIGNUM_LENGTH (x));
  1219.   bignum_type p = (bignum_allocate ((length_x + 1), negative_p));
  1220.   bignum_destructive_copy (x, p);
  1221.   (BIGNUM_REF (p, length_x)) = 0;
  1222.   bignum_destructive_scale_up (p, y);
  1223.   return (bignum_trim (p));
  1224. }
  1225.  
  1226. static void
  1227. DEFUN (bignum_destructive_scale_up, (bignum, factor),
  1228.        bignum_type bignum AND bignum_digit_type factor)
  1229. {
  1230.   fast bignum_digit_type carry = 0;
  1231.   fast bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
  1232.   fast bignum_digit_type two_digits;
  1233.   fast bignum_digit_type product_low;
  1234. #define product_high carry
  1235.   bignum_digit_type * end = (scan + (BIGNUM_LENGTH (bignum)));
  1236.   BIGNUM_ASSERT ((factor > 1) && (factor < BIGNUM_RADIX_ROOT));
  1237.   while (scan < end)
  1238.     {
  1239.       two_digits = (*scan);
  1240.       product_low = ((factor * (HD_LOW (two_digits))) + (HD_LOW (carry)));
  1241.       product_high =
  1242.     ((factor * (HD_HIGH (two_digits))) +
  1243.      (HD_HIGH (product_low)) +
  1244.      (HD_HIGH (carry)));
  1245.       (*scan++) = (HD_CONS ((HD_LOW (product_high)), (HD_LOW (product_low))));
  1246.       carry = (HD_HIGH (product_high));
  1247.     }
  1248.   /* A carry here would be an overflow, i.e. it would not fit.
  1249.      Hopefully the callers allocate enough space that this will
  1250.      never happen.
  1251.    */
  1252.   BIGNUM_ASSERT (carry == 0);
  1253.   return;
  1254. #undef product_high
  1255. }
  1256.  
  1257. static void
  1258. DEFUN (bignum_destructive_add, (bignum, n),
  1259.        bignum_type bignum AND bignum_digit_type n)
  1260. {
  1261.   fast bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
  1262.   fast bignum_digit_type digit;
  1263.   digit = ((*scan) + n);
  1264.   if (digit < BIGNUM_RADIX)
  1265.     {
  1266.       (*scan) = digit;
  1267.       return;
  1268.     }
  1269.   (*scan++) = (digit - BIGNUM_RADIX);
  1270.   while (1)
  1271.     {
  1272.       digit = ((*scan) + 1);
  1273.       if (digit < BIGNUM_RADIX)
  1274.     {
  1275.       (*scan) = digit;
  1276.       return;
  1277.     }
  1278.       (*scan++) = (digit - BIGNUM_RADIX);
  1279.     }
  1280. }
  1281.  
  1282. /* Division */
  1283.  
  1284. /* For help understanding this algorithm, see:
  1285.    Knuth, Donald E., "The Art of Computer Programming",
  1286.    volume 2, "Seminumerical Algorithms"
  1287.    section 4.3.1, "Multiple-Precision Arithmetic". */
  1288.  
  1289. static void
  1290. DEFUN (bignum_divide_unsigned_large_denominator, (numerator, denominator,
  1291.                           quotient, remainder,
  1292.                           q_negative_p, r_negative_p),
  1293.        bignum_type numerator
  1294.        AND bignum_type denominator
  1295.        AND bignum_type * quotient
  1296.        AND bignum_type * remainder
  1297.        AND int q_negative_p
  1298.        AND int r_negative_p)
  1299. {
  1300.   bignum_length_type length_n = ((BIGNUM_LENGTH (numerator)) + 1);
  1301.   bignum_length_type length_d = (BIGNUM_LENGTH (denominator));
  1302.   bignum_type q =
  1303.     ((quotient != ((bignum_type *) 0))
  1304.      ? (bignum_allocate ((length_n - length_d), q_negative_p))
  1305.      : BIGNUM_OUT_OF_BAND);
  1306.   bignum_type u = (bignum_allocate (length_n, r_negative_p));
  1307.   int shift = 0;
  1308.   BIGNUM_ASSERT (length_d > 1);
  1309.   {
  1310.     fast bignum_digit_type v1 = (BIGNUM_REF ((denominator), (length_d - 1)));
  1311.     while (v1 < (BIGNUM_RADIX / 2))
  1312.       {
  1313.     v1 <<= 1;
  1314.     shift += 1;
  1315.       }
  1316.   }
  1317.   if (shift == 0)
  1318.     {
  1319.       bignum_destructive_copy (numerator, u);
  1320.       (BIGNUM_REF (u, (length_n - 1))) = 0;
  1321.       bignum_divide_unsigned_normalized (u, denominator, q);
  1322.     }
  1323.   else
  1324.     {
  1325.       bignum_type v = (bignum_allocate (length_d, 0));
  1326.       bignum_destructive_normalization (numerator, u, shift);
  1327.       bignum_destructive_normalization (denominator, v, shift);
  1328.       bignum_divide_unsigned_normalized (u, v, q);
  1329.       BIGNUM_DEALLOCATE (v);
  1330.       if (remainder != ((bignum_type *) 0))
  1331.     bignum_destructive_unnormalization (u, shift);
  1332.     }
  1333.   if (quotient != ((bignum_type *) 0))
  1334.     (*quotient) = (bignum_trim (q));
  1335.   if (remainder != ((bignum_type *) 0))
  1336.     (*remainder) = (bignum_trim (u));
  1337.   else
  1338.     BIGNUM_DEALLOCATE (u);
  1339.   return;
  1340. }
  1341.  
  1342. static void
  1343. DEFUN (bignum_divide_unsigned_normalized, (u, v, q),
  1344.        bignum_type u AND bignum_type v AND bignum_type q)
  1345. {
  1346.   bignum_length_type u_length = (BIGNUM_LENGTH (u));
  1347.   bignum_length_type v_length = (BIGNUM_LENGTH (v));
  1348.   bignum_digit_type * u_start = (BIGNUM_START_PTR (u));
  1349.   bignum_digit_type * u_scan = (u_start + u_length);
  1350.   bignum_digit_type * u_scan_limit = (u_start + v_length);
  1351.   bignum_digit_type * u_scan_start = (u_scan - v_length);
  1352.   bignum_digit_type * v_start = (BIGNUM_START_PTR (v));
  1353.   bignum_digit_type * v_end = (v_start + v_length);
  1354.   bignum_digit_type * q_scan = 0;
  1355.   bignum_digit_type v1 = (v_end[-1]);
  1356.   bignum_digit_type v2 = (v_end[-2]);
  1357.   fast bignum_digit_type ph;    /* high half of double-digit product */
  1358.   fast bignum_digit_type pl;    /* low half of double-digit product */
  1359.   fast bignum_digit_type guess;
  1360.   fast bignum_digit_type gh;    /* high half-digit of guess */
  1361.   fast bignum_digit_type ch;    /* high half of double-digit comparand */
  1362.   fast bignum_digit_type v2l = (HD_LOW (v2));
  1363.   fast bignum_digit_type v2h = (HD_HIGH (v2));
  1364.   fast bignum_digit_type cl;    /* low half of double-digit comparand */
  1365. #define gl ph            /* low half-digit of guess */
  1366. #define uj pl
  1367. #define qj ph
  1368.   bignum_digit_type gm;        /* memory loc for reference parameter */
  1369.   if (q != BIGNUM_OUT_OF_BAND)
  1370.     q_scan = ((BIGNUM_START_PTR (q)) + (BIGNUM_LENGTH (q)));
  1371.   while (u_scan_limit < u_scan)
  1372.     {
  1373.       uj = (*--u_scan);
  1374.       if (uj != v1)
  1375.     {
  1376.       /* comparand =
  1377.          (((((uj * BIGNUM_RADIX) + uj1) % v1) * BIGNUM_RADIX) + uj2);
  1378.          guess = (((uj * BIGNUM_RADIX) + uj1) / v1); */
  1379.       cl = (u_scan[-2]);
  1380.       ch = (bignum_digit_divide (uj, (u_scan[-1]), v1, (&gm)));
  1381.       guess = gm;
  1382.     }
  1383.       else
  1384.     {
  1385.       cl = (u_scan[-2]);
  1386.       ch = ((u_scan[-1]) + v1);
  1387.       guess = (BIGNUM_RADIX - 1);
  1388.     }
  1389.       while (1)
  1390.     {
  1391.       /* product = (guess * v2); */
  1392.       gl = (HD_LOW (guess));
  1393.       gh = (HD_HIGH (guess));
  1394.       pl = (v2l * gl);
  1395.       ph = ((v2l * gh) + (v2h * gl) + (HD_HIGH (pl)));
  1396.       pl = (HD_CONS ((HD_LOW (ph)), (HD_LOW (pl))));
  1397.       ph = ((v2h * gh) + (HD_HIGH (ph)));
  1398.       /* if (comparand >= product) */
  1399.       if ((ch > ph) || ((ch == ph) && (cl >= pl)))
  1400.         break;
  1401.       guess -= 1;
  1402.       /* comparand += (v1 << BIGNUM_DIGIT_LENGTH) */
  1403.       ch += v1;
  1404.       /* if (comparand >= (BIGNUM_RADIX * BIGNUM_RADIX)) */
  1405.       if (ch >= BIGNUM_RADIX)
  1406.         break;
  1407.     }
  1408.       qj = (bignum_divide_subtract (v_start, v_end, guess, (--u_scan_start)));
  1409.       if (q != BIGNUM_OUT_OF_BAND)
  1410.     (*--q_scan) = qj;
  1411.     }
  1412.   return;
  1413. #undef gl
  1414. #undef uj
  1415. #undef qj
  1416. }
  1417.  
  1418. static bignum_digit_type
  1419. DEFUN (bignum_divide_subtract, (v_start, v_end, guess, u_start),
  1420.        bignum_digit_type * v_start
  1421.        AND bignum_digit_type * v_end
  1422.        AND bignum_digit_type guess
  1423.        AND bignum_digit_type * u_start)
  1424. {
  1425.   bignum_digit_type * v_scan = v_start;
  1426.   bignum_digit_type * u_scan = u_start;
  1427.   fast bignum_digit_type carry = 0;
  1428.   if (guess == 0) return (0);
  1429.   {
  1430.     bignum_digit_type gl = (HD_LOW (guess));
  1431.     bignum_digit_type gh = (HD_HIGH (guess));
  1432.     fast bignum_digit_type v;
  1433.     fast bignum_digit_type pl;
  1434.     fast bignum_digit_type vl;
  1435. #define vh v
  1436. #define ph carry
  1437. #define diff pl
  1438.     while (v_scan < v_end)
  1439.       {
  1440.     v = (*v_scan++);
  1441.     vl = (HD_LOW (v));
  1442.     vh = (HD_HIGH (v));
  1443.     pl = ((vl * gl) + (HD_LOW (carry)));
  1444.     ph = ((vl * gh) + (vh * gl) + (HD_HIGH (pl)) + (HD_HIGH (carry)));
  1445.     diff = ((*u_scan) - (HD_CONS ((HD_LOW (ph)), (HD_LOW (pl)))));
  1446.     if (diff < 0)
  1447.       {
  1448.         (*u_scan++) = (diff + BIGNUM_RADIX);
  1449.         carry = ((vh * gh) + (HD_HIGH (ph)) + 1);
  1450.       }
  1451.     else
  1452.       {
  1453.         (*u_scan++) = diff;
  1454.         carry = ((vh * gh) + (HD_HIGH (ph)));
  1455.       }
  1456.       }
  1457.     if (carry == 0)
  1458.       return (guess);
  1459.     diff = ((*u_scan) - carry);
  1460.     if (diff < 0)
  1461.       (*u_scan) = (diff + BIGNUM_RADIX);
  1462.     else
  1463.       {
  1464.     (*u_scan) = diff;
  1465.     return (guess);
  1466.       }
  1467. #undef vh
  1468. #undef ph
  1469. #undef diff
  1470.   }
  1471.   /* Subtraction generated carry, implying guess is one too large.
  1472.      Add v back in to bring it back down. */
  1473.   v_scan = v_start;
  1474.   u_scan = u_start;
  1475.   carry = 0;
  1476.   while (v_scan < v_end)
  1477.     {
  1478.       bignum_digit_type sum = ((*v_scan++) + (*u_scan) + carry);
  1479.       if (sum < BIGNUM_RADIX)
  1480.     {
  1481.       (*u_scan++) = sum;
  1482.       carry = 0;
  1483.     }
  1484.       else
  1485.     {
  1486.       (*u_scan++) = (sum - BIGNUM_RADIX);
  1487.       carry = 1;
  1488.     }
  1489.     }
  1490.   if (carry == 1)
  1491.     {
  1492.       bignum_digit_type sum = ((*u_scan) + carry);
  1493.       (*u_scan) = ((sum < BIGNUM_RADIX) ? sum : (sum - BIGNUM_RADIX));
  1494.     }
  1495.   return (guess - 1);
  1496. }
  1497.  
  1498. static void
  1499. DEFUN (bignum_divide_unsigned_medium_denominator, (numerator, denominator,
  1500.                            quotient, remainder,
  1501.                            q_negative_p, r_negative_p),
  1502.        bignum_type numerator
  1503.        AND bignum_digit_type denominator
  1504.        AND bignum_type * quotient
  1505.        AND bignum_type * remainder
  1506.        AND int q_negative_p
  1507.        AND int r_negative_p)
  1508. {
  1509.   bignum_length_type length_n = (BIGNUM_LENGTH (numerator));
  1510.   bignum_length_type length_q;
  1511.   bignum_type q;
  1512.   int shift = 0;
  1513.   /* Because `bignum_digit_divide' requires a normalized denominator. */
  1514.   while (denominator < (BIGNUM_RADIX / 2))
  1515.     {
  1516.       denominator <<= 1;
  1517.       shift += 1;
  1518.     }
  1519.   if (shift == 0)
  1520.     {
  1521.       length_q = length_n;
  1522.       q = (bignum_allocate (length_q, q_negative_p));
  1523.       bignum_destructive_copy (numerator, q);
  1524.     }
  1525.   else
  1526.     {
  1527.       length_q = (length_n + 1);
  1528.       q = (bignum_allocate (length_q, q_negative_p));
  1529.       bignum_destructive_normalization (numerator, q, shift);
  1530.     }
  1531.   {
  1532.     fast bignum_digit_type r = 0;
  1533.     fast bignum_digit_type * start = (BIGNUM_START_PTR (q));
  1534.     fast bignum_digit_type * scan = (start + length_q);
  1535.     bignum_digit_type qj;
  1536.     if (quotient != ((bignum_type *) 0))
  1537.       {
  1538.     while (start < scan)
  1539.       {
  1540.         r = (bignum_digit_divide (r, (*--scan), denominator, (&qj)));
  1541.         (*scan) = qj;
  1542.       }
  1543.     (*quotient) = (bignum_trim (q));
  1544.       }
  1545.     else
  1546.       {
  1547.     while (start < scan)
  1548.       r = (bignum_digit_divide (r, (*--scan), denominator, (&qj)));
  1549.     BIGNUM_DEALLOCATE (q);
  1550.       }
  1551.     if (remainder != ((bignum_type *) 0))
  1552.       {
  1553.     if (shift != 0)
  1554.       r >>= shift;
  1555.     (*remainder) = (bignum_digit_to_bignum (r, r_negative_p));
  1556.       }
  1557.   }
  1558.   return;
  1559. }
  1560.  
  1561. static void
  1562. DEFUN (bignum_destructive_normalization, (source, target, shift_left),
  1563.        bignum_type source AND bignum_type target AND int shift_left)
  1564. {
  1565.   fast bignum_digit_type digit;
  1566.   fast bignum_digit_type * scan_source = (BIGNUM_START_PTR (source));
  1567.   fast bignum_digit_type carry = 0;
  1568.   fast bignum_digit_type * scan_target = (BIGNUM_START_PTR (target));
  1569.   bignum_digit_type * end_source = (scan_source + (BIGNUM_LENGTH (source)));
  1570.   bignum_digit_type * end_target = (scan_target + (BIGNUM_LENGTH (target)));
  1571.   int shift_right = (BIGNUM_DIGIT_LENGTH - shift_left);
  1572.   bignum_digit_type mask = ((1L << shift_right) - 1);
  1573.   while (scan_source < end_source)
  1574.     {
  1575.       digit = (*scan_source++);
  1576.       (*scan_target++) = (((digit & mask) << shift_left) | carry);
  1577.       carry = (digit >> shift_right);
  1578.     }
  1579.   if (scan_target < end_target)
  1580.     (*scan_target) = carry;
  1581.   else
  1582.     BIGNUM_ASSERT (carry == 0);
  1583.   return;
  1584. }
  1585.  
  1586. static void
  1587. DEFUN (bignum_destructive_unnormalization, (bignum, shift_right),
  1588.        bignum_type bignum AND int shift_right)
  1589. {
  1590.   bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  1591.   fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
  1592.   fast bignum_digit_type digit;
  1593.   fast bignum_digit_type carry = 0;
  1594.   int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
  1595.   bignum_digit_type mask = ((1L << shift_right) - 1);
  1596.   while (start < scan)
  1597.     {
  1598.       digit = (*--scan);
  1599.       (*scan) = ((digit >> shift_right) | carry);
  1600.       carry = ((digit & mask) << shift_left);
  1601.     }
  1602.   BIGNUM_ASSERT (carry == 0);
  1603.   return;
  1604. }
  1605.  
  1606. /* This is a reduced version of the division algorithm, applied to the
  1607.    case of dividing two bignum digits by one bignum digit.  It is
  1608.    assumed that the numerator and denominator are normalized. */
  1609.  
  1610. #define BDD_STEP(qn, j)                            \
  1611. {                                    \
  1612.   uj = (u[j]);                                \
  1613.   if (uj != v1)                                \
  1614.     {                                    \
  1615.       uj_uj1 = (HD_CONS (uj, (u[j + 1])));                \
  1616.       guess = (uj_uj1 / v1);                        \
  1617.       comparand = (HD_CONS ((uj_uj1 % v1), (u[j + 2])));        \
  1618.     }                                    \
  1619.   else                                    \
  1620.     {                                    \
  1621.       guess = (BIGNUM_RADIX_ROOT - 1);                    \
  1622.       comparand = (HD_CONS (((u[j + 1]) + v1), (u[j + 2])));        \
  1623.     }                                    \
  1624.   while ((guess * v2) > comparand)                    \
  1625.     {                                    \
  1626.       guess -= 1;                            \
  1627.       comparand += (v1 << BIGNUM_HALF_DIGIT_LENGTH);            \
  1628.       if (comparand >= BIGNUM_RADIX)                    \
  1629.     break;                                \
  1630.     }                                    \
  1631.   qn = (bignum_digit_divide_subtract (v1, v2, guess, (&u[j])));        \
  1632. }
  1633.  
  1634. static bignum_digit_type
  1635. DEFUN (bignum_digit_divide, (uh, ul, v, q),
  1636.        bignum_digit_type uh AND bignum_digit_type ul
  1637.        AND bignum_digit_type v AND bignum_digit_type * q) /* return value */
  1638. {
  1639.   fast bignum_digit_type guess;
  1640.   fast bignum_digit_type comparand;
  1641.   fast bignum_digit_type v1 = (HD_HIGH (v));
  1642.   fast bignum_digit_type v2 = (HD_LOW (v));
  1643.   fast bignum_digit_type uj;
  1644.   fast bignum_digit_type uj_uj1;
  1645.   bignum_digit_type q1;
  1646.   bignum_digit_type q2;
  1647.   bignum_digit_type u [4];
  1648.   if (uh == 0)
  1649.     {
  1650.       if (ul < v)
  1651.     {
  1652.       (*q) = 0;
  1653.       return (ul);
  1654.     }
  1655.       else if (ul == v)
  1656.     {
  1657.       (*q) = 1;
  1658.       return (0);
  1659.     }
  1660.     }
  1661.   (u[0]) = (HD_HIGH (uh));
  1662.   (u[1]) = (HD_LOW (uh));
  1663.   (u[2]) = (HD_HIGH (ul));
  1664.   (u[3]) = (HD_LOW (ul));
  1665.   v1 = (HD_HIGH (v));
  1666.   v2 = (HD_LOW (v));
  1667.   BDD_STEP (q1, 0);
  1668.   BDD_STEP (q2, 1);
  1669.   (*q) = (HD_CONS (q1, q2));
  1670.   return (HD_CONS ((u[2]), (u[3])));
  1671. }
  1672.  
  1673. #undef BDD_STEP
  1674.  
  1675. #define BDDS_MULSUB(vn, un, carry_in)                    \
  1676. {                                    \
  1677.   product = ((vn * guess) + carry_in);                    \
  1678.   diff = (un - (HD_LOW (product)));                    \
  1679.   if (diff < 0)                                \
  1680.     {                                    \
  1681.       un = (diff + BIGNUM_RADIX_ROOT);                    \
  1682.       carry = ((HD_HIGH (product)) + 1);                \
  1683.     }                                    \
  1684.   else                                    \
  1685.     {                                    \
  1686.       un = diff;                            \
  1687.       carry = (HD_HIGH (product));                    \
  1688.     }                                    \
  1689. }
  1690.  
  1691. #define BDDS_ADD(vn, un, carry_in)                    \
  1692. {                                    \
  1693.   sum = (vn + un + carry_in);                        \
  1694.   if (sum < BIGNUM_RADIX_ROOT)                        \
  1695.     {                                    \
  1696.       un = sum;                                \
  1697.       carry = 0;                            \
  1698.     }                                    \
  1699.   else                                    \
  1700.     {                                    \
  1701.       un = (sum - BIGNUM_RADIX_ROOT);                    \
  1702.       carry = 1;                            \
  1703.     }                                    \
  1704. }
  1705.  
  1706. static bignum_digit_type
  1707. DEFUN (bignum_digit_divide_subtract, (v1, v2, guess, u),
  1708.        bignum_digit_type v1 AND bignum_digit_type v2
  1709.        AND bignum_digit_type guess AND bignum_digit_type * u)
  1710. {
  1711.   {
  1712.     fast bignum_digit_type product;
  1713.     fast bignum_digit_type diff;
  1714.     fast bignum_digit_type carry;
  1715.     BDDS_MULSUB (v2, (u[2]), 0);
  1716.     BDDS_MULSUB (v1, (u[1]), carry);
  1717.     if (carry == 0)
  1718.       return (guess);
  1719.     diff = ((u[0]) - carry);
  1720.     if (diff < 0)
  1721.       (u[0]) = (diff + BIGNUM_RADIX);
  1722.     else
  1723.       {
  1724.     (u[0]) = diff;
  1725.     return (guess);
  1726.       }
  1727.   }
  1728.   {
  1729.     fast bignum_digit_type sum;
  1730.     fast bignum_digit_type carry;
  1731.     BDDS_ADD(v2, (u[2]), 0);
  1732.     BDDS_ADD(v1, (u[1]), carry);
  1733.     if (carry == 1)
  1734.       (u[0]) += 1;
  1735.   }
  1736.   return (guess - 1);
  1737. }
  1738.  
  1739. #undef BDDS_MULSUB
  1740. #undef BDDS_ADD
  1741.  
  1742. static void
  1743. DEFUN (bignum_divide_unsigned_small_denominator, (numerator, denominator,
  1744.                           quotient, remainder,
  1745.                           q_negative_p, r_negative_p),
  1746.        bignum_type numerator
  1747.        AND bignum_digit_type denominator
  1748.        AND bignum_type * quotient
  1749.        AND bignum_type * remainder
  1750.        AND int q_negative_p
  1751.        AND int r_negative_p)
  1752. {
  1753.   bignum_type q = (bignum_new_sign (numerator, q_negative_p));
  1754.   bignum_digit_type r = (bignum_destructive_scale_down (q, denominator));
  1755.   (*quotient) = (bignum_trim (q));
  1756.   if (remainder != ((bignum_type *) 0))
  1757.     (*remainder) = (bignum_digit_to_bignum (r, r_negative_p));
  1758.   return;
  1759. }
  1760.  
  1761. /* Given (denominator > 1), it is fairly easy to show that
  1762.    (quotient_high < BIGNUM_RADIX_ROOT), after which it is easy to see
  1763.    that all digits are < BIGNUM_RADIX. */
  1764.  
  1765. static bignum_digit_type
  1766. DEFUN (bignum_destructive_scale_down, (bignum, denominator),
  1767.        bignum_type bignum AND fast bignum_digit_type denominator)
  1768. {
  1769.   fast bignum_digit_type numerator;
  1770.   fast bignum_digit_type remainder = 0;
  1771.   fast bignum_digit_type two_digits;
  1772. #define quotient_high remainder
  1773.   bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  1774.   bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
  1775.   BIGNUM_ASSERT ((denominator > 1) && (denominator < BIGNUM_RADIX_ROOT));
  1776.   while (start < scan)
  1777.     {
  1778.       two_digits = (*--scan);
  1779.       numerator = (HD_CONS (remainder, (HD_HIGH (two_digits))));
  1780.       quotient_high = (numerator / denominator);
  1781.       numerator = (HD_CONS ((numerator % denominator), (HD_LOW (two_digits))));
  1782.       (*scan) = (HD_CONS (quotient_high, (numerator / denominator)));
  1783.       remainder = (numerator % denominator);
  1784.     }
  1785.   return (remainder);
  1786. #undef quotient_high
  1787. }
  1788.  
  1789. static bignum_type
  1790. DEFUN (bignum_remainder_unsigned_small_denominator, (n, d, negative_p),
  1791.        bignum_type n AND bignum_digit_type d AND int negative_p)
  1792. {
  1793.   fast bignum_digit_type two_digits;
  1794.   bignum_digit_type * start = (BIGNUM_START_PTR (n));
  1795.   fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (n)));
  1796.   fast bignum_digit_type r = 0;
  1797.   BIGNUM_ASSERT ((d > 1) && (d < BIGNUM_RADIX_ROOT));
  1798.   while (start < scan)
  1799.     {
  1800.       two_digits = (*--scan);
  1801.       r =
  1802.     ((HD_CONS (((HD_CONS (r, (HD_HIGH (two_digits)))) % d),
  1803.            (HD_LOW (two_digits))))
  1804.      % d);
  1805.     }
  1806.   return (bignum_digit_to_bignum (r, negative_p));
  1807. }
  1808.  
  1809. static bignum_type
  1810. DEFUN (bignum_digit_to_bignum, (digit, negative_p),
  1811.        fast bignum_digit_type digit AND int negative_p)
  1812. {
  1813.   if (digit == 0)
  1814.     return (BIGNUM_ZERO ());
  1815.   else
  1816.     {
  1817.       fast bignum_type result = (bignum_allocate (1, negative_p));
  1818.       (BIGNUM_REF (result, 0)) = digit;
  1819.       return (result);
  1820.     }
  1821. }
  1822.  
  1823. /* Allocation */
  1824.  
  1825. static bignum_type
  1826. DEFUN (bignum_allocate, (length, negative_p),
  1827.        fast bignum_length_type length AND int negative_p)
  1828. {
  1829.   BIGNUM_ASSERT ((length >= 0) || (length < BIGNUM_RADIX));
  1830.   {
  1831.     fast bignum_type result = (BIGNUM_ALLOCATE (length));
  1832.     BIGNUM_SET_HEADER (result, length, negative_p);
  1833.     return (result);
  1834.   }
  1835. }
  1836.  
  1837. static bignum_type
  1838. DEFUN (bignum_allocate_zeroed, (length, negative_p),
  1839.        fast bignum_length_type length AND int negative_p)
  1840. {
  1841.   BIGNUM_ASSERT ((length >= 0) || (length < BIGNUM_RADIX));
  1842.   {
  1843.     fast bignum_type result = (BIGNUM_ALLOCATE (length));
  1844.     fast bignum_digit_type * scan = (BIGNUM_START_PTR (result));
  1845.     fast bignum_digit_type * end = (scan + length);
  1846.     BIGNUM_SET_HEADER (result, length, negative_p);
  1847.     while (scan < end)
  1848.       (*scan++) = 0;
  1849.     return (result);
  1850.   }
  1851. }
  1852.  
  1853. static bignum_type
  1854. DEFUN (bignum_shorten_length, (bignum, length),
  1855.        fast bignum_type bignum AND fast bignum_length_type length)
  1856. {
  1857.   fast bignum_length_type current_length = (BIGNUM_LENGTH (bignum));
  1858.   BIGNUM_ASSERT ((length >= 0) || (length <= current_length));
  1859.   if (length < current_length)
  1860.     {
  1861.       BIGNUM_SET_HEADER
  1862.     (bignum, length, ((length != 0) && (BIGNUM_NEGATIVE_P (bignum))));
  1863.       BIGNUM_REDUCE_LENGTH (bignum, bignum, length);
  1864.     }
  1865.   return (bignum);
  1866. }
  1867.  
  1868. static bignum_type
  1869. DEFUN (bignum_trim, (bignum), bignum_type bignum)
  1870. {
  1871.   fast bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  1872.   fast bignum_digit_type * end = (start + (BIGNUM_LENGTH (bignum)));
  1873.   fast bignum_digit_type * scan = end;
  1874.   while ((start <= scan) && ((*--scan) == 0))
  1875.     ;
  1876.   scan += 1;
  1877.   if (scan < end)
  1878.     {
  1879.       fast bignum_length_type length = (scan - start);
  1880.       BIGNUM_SET_HEADER
  1881.     (bignum, length, ((length != 0) && (BIGNUM_NEGATIVE_P (bignum))));
  1882.       BIGNUM_REDUCE_LENGTH (bignum, bignum, length);
  1883.     }
  1884.   return (bignum);
  1885. }
  1886.  
  1887. /* Copying */
  1888.  
  1889. static bignum_type
  1890. DEFUN (bignum_copy, (source), fast bignum_type source)
  1891. {
  1892.   fast bignum_type target =
  1893.     (bignum_allocate ((BIGNUM_LENGTH (source)), (BIGNUM_NEGATIVE_P (source))));
  1894.   bignum_destructive_copy (source, target);
  1895.   return (target);
  1896. }
  1897.  
  1898. static bignum_type
  1899. DEFUN (bignum_new_sign, (bignum, negative_p),
  1900.        fast bignum_type bignum AND int negative_p)
  1901. {
  1902.   fast bignum_type result =
  1903.     (bignum_allocate ((BIGNUM_LENGTH (bignum)), negative_p));
  1904.   bignum_destructive_copy (bignum, result);
  1905.   return (result);
  1906. }
  1907.  
  1908. static bignum_type
  1909. DEFUN (bignum_maybe_new_sign, (bignum, negative_p),
  1910.        fast bignum_type bignum AND int negative_p)
  1911. {
  1912. #ifndef BIGNUM_FORCE_NEW_RESULTS
  1913.   if ((BIGNUM_NEGATIVE_P (bignum)) ? negative_p : (! negative_p))
  1914.     return (bignum);
  1915.   else
  1916. #endif /* not BIGNUM_FORCE_NEW_RESULTS */
  1917.     {
  1918.       fast bignum_type result =
  1919.     (bignum_allocate ((BIGNUM_LENGTH (bignum)), negative_p));
  1920.       bignum_destructive_copy (bignum, result);
  1921.       return (result);
  1922.     }
  1923. }
  1924.  
  1925. static void
  1926. DEFUN (bignum_destructive_copy, (source, target),
  1927.        bignum_type source AND bignum_type target)
  1928. {
  1929.   fast bignum_digit_type * scan_source = (BIGNUM_START_PTR (source));
  1930.   fast bignum_digit_type * end_source =
  1931.     (scan_source + (BIGNUM_LENGTH (source)));
  1932.   fast bignum_digit_type * scan_target = (BIGNUM_START_PTR (target));
  1933.   while (scan_source < end_source)
  1934.     (*scan_target++) = (*scan_source++);
  1935.   return;
  1936. }
  1937.