l4/pkg/ocaml/contrib/otherlibs/num/bng_digit.c

   1 /***********************************************************************/
   2 /*                                                                     */
   3 /*                           Objective Caml                            */
   4 /*                                                                     */
   5 /*            Xavier Leroy, projet Cristal, INRIA Rocquencourt         */
   6 /*                                                                     */
   7 /*  Copyright 2003 Institut National de Recherche en Informatique et   */
   8 /*  en Automatique.  All rights reserved.  This file is distributed    */
   9 /*  under the terms of the GNU Library General Public License, with    */
  10 /*  the special exception on linking described in file ../../LICENSE.  */
  11 /*                                                                     */
  12 /***********************************************************************/
  13
  14 /* $Id: bng_digit.c 5880 2003-10-24 09:18:01Z xleroy $ */
  15
  16 /**** Generic operations on digits ****/
  17
  18 /* These macros can be defined in the machine-specific include file.
  19    Below are the default definitions (in plain C).
  20    Except for BngMult, all macros are guaranteed to evaluate their
  21    arguments exactly once. */
  22
  23 #ifndef BngAdd2
  24 /* res = arg1 + arg2.  carryout = carry out. */
  25 #define BngAdd2(res,carryout,arg1,arg2) {                                   \
  26   bngdigit tmp1, tmp2;                                                      \
  27   tmp1 = arg1;                                                              \
  28   tmp2 = tmp1 + (arg2);                                                     \
  29   carryout = (tmp2 < tmp1);                                                 \
  30   res = tmp2;                                                               \
  31 }
  32 #endif
  33
  34 #ifndef BngAdd2Carry
  35 /* res = arg1 + arg2 + carryin.  carryout = carry out. */
  36 #define BngAdd2Carry(res,carryout,arg1,arg2,carryin) {                      \
  37   bngdigit tmp1, tmp2, tmp3;                                                \
  38   tmp1 = arg1;                                                              \
  39   tmp2 = tmp1 + (arg2);                                                     \
  40   tmp3 = tmp2 + (carryin);                                                  \
  41   carryout = (tmp2 < tmp1) + (tmp3 < tmp2);                                 \
  42   res = tmp3;                                                               \
  43 }
  44 #endif
  45
  46 #ifndef BngAdd3
  47 /* res = arg1 + arg2 + arg3.  Each carry increments carryaccu. */
  48 #define BngAdd3(res,carryaccu,arg1,arg2,arg3) {                             \
  49   bngdigit tmp1, tmp2, tmp3;                                                \
  50   tmp1 = arg1;                                                              \
  51   tmp2 = tmp1 + (arg2);                                                     \
  52   carryaccu += (tmp2 < tmp1);                                               \
  53   tmp3 = tmp2 + (arg3);                                                     \
  54   carryaccu += (tmp3 < tmp2);                                               \
  55   res = tmp3;                                                               \
  56 }
  57 #endif
  58
  59 #ifndef BngSub2
  60 /* res = arg1 - arg2.  carryout = carry out. */
  61 #define BngSub2(res,carryout,arg1,arg2) {                                   \
  62   bngdigit tmp1, tmp2;                                                      \
  63   tmp1 = arg1;                                                              \
  64   tmp2 = arg2;                                                              \
  65   res = tmp1 - tmp2;                                                        \
  66   carryout = (tmp1 < tmp2);                                                 \
  67 }
  68 #endif
  69
  70 #ifndef BngSub2Carry
  71 /* res = arg1 - arg2 - carryin.  carryout = carry out. */
  72 #define BngSub2Carry(res,carryout,arg1,arg2,carryin) {                      \
  73   bngdigit tmp1, tmp2, tmp3;                                                \
  74   tmp1 = arg1;                                                              \
  75   tmp2 = arg2;                                                              \
  76   tmp3 = tmp1 - tmp2;                                                       \
  77   res = tmp3 - (carryin);                                                   \
  78   carryout = (tmp1 < tmp2) + (tmp3 < carryin);                              \
  79 }
  80 #endif
  81
  82 #ifndef BngSub3
  83 /* res = arg1 - arg2 - arg3.  Each carry increments carryaccu. */
  84 #define BngSub3(res,carryaccu,arg1,arg2,arg3) {                             \
  85   bngdigit tmp1, tmp2, tmp3, tmp4;                                          \
  86   tmp1 = arg1;                                                              \
  87   tmp2 = arg2;                                                              \
  88   tmp3 = arg3;                                                              \
  89   tmp4 = tmp1 - tmp2;                                                       \
  90   res = tmp4 - tmp3;                                                        \
  91   carryaccu += (tmp1 < tmp2) + (tmp4 < tmp3);                               \
  92 }
  93 #endif
  94
  95 #define BngLowHalf(d) ((d) & ((1L << BNG_BITS_PER_HALF_DIGIT) - 1))
  96 #define BngHighHalf(d) ((d) >> BNG_BITS_PER_HALF_DIGIT)
  97
  98 #ifndef BngMult
  99 /* resl = low  digit of product arg1 * arg2
 100    resh = high digit of product arg1 * arg2. */
 101 #define BngMult(resh,resl,arg1,arg2) {                                      \
 102   bngdigit p11 = BngLowHalf(arg1) * BngLowHalf(arg2);                       \
 103   bngdigit p12 = BngLowHalf(arg1) * BngHighHalf(arg2);                      \
 104   bngdigit p21 = BngHighHalf(arg1) * BngLowHalf(arg2);                      \
 105   bngdigit p22 = BngHighHalf(arg1) * BngHighHalf(arg2);                     \
 106   resh = p22 + (p12 >> BNG_BITS_PER_HALF_DIGIT)                             \
 107              + (p21 >> BNG_BITS_PER_HALF_DIGIT);                            \
 108   BngAdd3(resl, resh,                                                       \
 109      p11, p12 << BNG_BITS_PER_HALF_DIGIT, p21 << BNG_BITS_PER_HALF_DIGIT);  \
 110 }
 111 #endif
 112
 113 #ifndef BngDiv
 114 /* Divide the double-width number nh:nl by d.
 115    Require d != 0 and nh < d.
 116    Store quotient in quo, remainder in rem.
 117    Can be slow if d is not normalized. */
 118 #define BngDiv(quo,rem,nh,nl,d) bng_div_aux(&(quo),&(rem),nh,nl,d)
 119 #define BngDivNeedsNormalization
 120
 121 static void bng_div_aux(bngdigit * quo, bngdigit * rem,
 122                         bngdigit nh, bngdigit nl, bngdigit d)
 123 {
 124   bngdigit dl, dh, ql, qh, pl, ph, nsaved;
 125
 126   dl = BngLowHalf(d);
 127   dh = BngHighHalf(d);
 128   /* Under-estimate the top half of the quotient (qh) */
 129   qh = nh / (dh + 1);
 130   /* Shift nh:nl right by BNG_BITS_PER_HALF_DIGIT bits,
 131      so that we focus on the top 1.5 digits of the numerator.
 132      Then, subtract (qh * d) from nh:nl. */
 133   nsaved = BngLowHalf(nl);
 134   ph = qh * dh;
 135   pl = qh * dl;
 136   nh -= ph; /* Subtract before shifting so that carry propagates for free */
 137   nl = (nl >> BNG_BITS_PER_HALF_DIGIT) | (nh << BNG_BITS_PER_HALF_DIGIT);
 138   nh = (nh >> BNG_BITS_PER_HALF_DIGIT);
 139   nh -= (nl < pl);  /* Borrow */
 140   nl -= pl;
 141   /* Adjust estimate qh until nh:nl < 0:d */
 142   while (nh != 0 || nl >= d) {
 143     nh -= (nl < d); /* Borrow */
 144     nl -= d;
 145     qh++;
 146   }
 147   /* Under-estimate the bottom half of the quotient (ql) */
 148   ql = nl / (dh + 1);
 149   /* Shift nh:nl left by BNG_BITS_PER_HALF_DIGIT bits, restoring the
 150      low bits we saved earlier, so that we focus on the bottom 1.5 digit
 151      of the numerator.  Then, subtract (ql * d) from nh:nl. */
 152   ph = ql * dh;
 153   pl = ql * dl;
 154   nl -= ph; /* Subtract before shifting so that carry propagates for free */
 155   nh = (nl >> BNG_BITS_PER_HALF_DIGIT);
 156   nl = (nl << BNG_BITS_PER_HALF_DIGIT) | nsaved;
 157   nh -= (nl < pl);  /* Borrow */
 158   nl -= pl;
 159   /* Adjust estimate ql until nh:nl < 0:d */
 160   while (nh != 0 || nl >= d) {
 161     nh -= (nl < d); /* Borrow */
 162     nl -= d;
 163     ql++;
 164   }
 165   /* We're done */
 166   *quo = (qh << BNG_BITS_PER_HALF_DIGIT) | ql;
 167   *rem = nl;
 168 }
 169
 170 #endif
 171