rts/gmp/mpn/pa64/udiv_qrnnd.c

   1 /*
   2 Copyright (C) 1999, 2000 Free Software Foundation, Inc.
   3
   4 This file is part of the GNU MP Library.
   5
   6 The GNU MP Library is free software; you can redistribute it and/or modify
   7 it under the terms of the GNU Lesser General Public License as published by
   8 the Free Software Foundation; either version 2.1 of the License, or (at your
   9 option) any later version.
  10
  11 The GNU MP Library is distributed in the hope that it will be useful, but
  12 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  13 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  14 License for more details.
  15
  16 You should have received a copy of the GNU Lesser General Public License
  17 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  18 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  19 MA 02111-1307, USA.
  20 */
  21
  22 #include "gmp.h"
  23 #include "gmp-impl.h"
  24 #include "longlong.h"
  25
  26 #define TWO64 18446744073709551616.0
  27
  28 mp_limb_t
  29 #if __STDC__
  30 __MPN(udiv_qrnnd) (mp_limb_t n1, mp_limb_t n0, mp_limb_t d, mp_limb_t *r)
  31 #else
  32 __MPN(udiv_qrnnd) (n1, n0, d, r)
  33      mp_limb_t n1;
  34      mp_limb_t n0;
  35      mp_limb_t d;
  36      mp_limb_t *r;
  37 #endif
  38 {
  39   mp_limb_t q1, q2, q;
  40   mp_limb_t p1, p0;
  41   double di, dq;
  42
  43   di = 1.0 / d;
  44
  45   /* Generate upper 53 bits of quotient.  Be careful here; the `double'
  46      quotient may be rounded to 2^64 which we cannot safely convert back
  47      to a 64-bit integer.  */
  48   dq = (TWO64 * (double) n1 + (double) n0) * di;
  49   if (dq >= TWO64)
  50     q1 = 0xfffffffffffff800LL;
  51   else
  52     q1 = (mp_limb_t) dq;
  53
  54   /* Multiply back in order to compare the product to the dividend.  */
  55   umul_ppmm (p1, p0, q1, d);
  56
  57   /* Was the 53-bit quotient greater that our sought quotient?  Test the
  58      sign of the partial remainder to find out.  */
  59   if (n1 < p1 || (n1 == p1 && n0 < p0))
  60     {
  61       /* 53-bit quotient too large.  Partial remainder is negative.
  62          Compute the absolute value of the remainder in n1,,n0.  */
  63       n1 = p1 - (n1 + (p0 < n0));
  64       n0 = p0 - n0;
  65
  66       /* Now use the partial remainder as new dividend to compute more bits of
  67          quotient.  This is an adjustment for the one we got previously.  */
  68       q2 = (mp_limb_t) ((TWO64 * (double) n1 + (double) n0) * di);
  69       umul_ppmm (p1, p0, q2, d);
  70
  71       q = q1 - q2;
  72       if (n1 < p1 || (n1 == p1 && n0 <= p0))
  73         {
  74           n0 = p0 - n0;
  75         }
  76       else
  77         {
  78           n0 = p0 - n0;
  79           n0 += d;
  80           q--;
  81         }
  82     }
  83   else
  84     {
  85       n1 = n1 - (p1 + (n0 < p0));
  86       n0 = n0 - p0;
  87
  88       q2 = (mp_limb_t) ((TWO64 * (double) n1 + (double) n0) * di);
  89       umul_ppmm (p1, p0, q2, d);
  90
  91       q = q1 + q2;
  92       if (n1 < p1 || (n1 == p1 && n0 < p0))
  93         {
  94           n0 = n0 - p0;
  95           n0 += d;
  96           q--;
  97         }
  98       else
  99         {
 100           n0 = n0 - p0;
 101           if (n0 >= d)
 102             {
 103               n0 -= d;
 104               q++;
 105             }
 106         }
 107     }
 108
 109   *r = n0;
 110   return q;
 111 }