rts/gmp/mpn/pa64w/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 #define TWO63 9223372036854775808.0
  28
  29 mp_limb_t
  30 #if __STDC__
  31 __MPN(udiv_qrnnd) (mp_limb_t n1, mp_limb_t n0, mp_limb_t d, mp_limb_t *r)
  32 #else
  33 __MPN(udiv_qrnnd) (n1, n0, d, r)
  34      mp_limb_t n1;
  35      mp_limb_t n0;
  36      mp_limb_t d;
  37      mp_limb_t *r;
  38 #endif
  39 {
  40   mp_limb_t q1, q2, q;
  41   mp_limb_t p1, p0;
  42   double di, dq;
  43
  44   di = 1.0 / d;
  45
  46   /* Generate upper 53 bits of quotient.  Be careful here; the `double'
  47      quotient may be rounded to 2^64 which we cannot safely convert back
  48      to a 64-bit integer.  */
  49   dq = (TWO64 * (double) n1 + (double) n0) * di;
  50   if (dq >= TWO64)
  51     q1 = 0xfffffffffffff800L;
  52 #ifndef __GNUC__
  53   /* Work around HP compiler bug.  */
  54   else if (dq > TWO63)
  55     q1 = (mp_limb_t) (dq - TWO63) + 0x8000000000000000L;
  56 #endif
  57   else
  58     q1 = (mp_limb_t) dq;
  59
  60   /* Multiply back in order to compare the product to the dividend.  */
  61   umul_ppmm (p1, p0, q1, d);
  62
  63   /* Was the 53-bit quotient greater that our sought quotient?  Test the
  64      sign of the partial remainder to find out.  */
  65   if (n1 < p1 || (n1 == p1 && n0 < p0))
  66     {
  67       /* 53-bit quotient too large.  Partial remainder is negative.
  68          Compute the absolute value of the remainder in n1,,n0.  */
  69       n1 = p1 - (n1 + (p0 < n0));
  70       n0 = p0 - n0;
  71
  72       /* Now use the partial remainder as new dividend to compute more bits of
  73          quotient.  This is an adjustment for the one we got previously.  */
  74       q2 = (mp_limb_t) ((TWO64 * (double) n1 + (double) n0) * di);
  75       umul_ppmm (p1, p0, q2, d);
  76
  77       q = q1 - q2;
  78       if (n1 < p1 || (n1 == p1 && n0 <= p0))
  79         {
  80           n0 = p0 - n0;
  81         }
  82       else
  83         {
  84           n0 = p0 - n0;
  85           n0 += d;
  86           q--;
  87         }
  88     }
  89   else
  90     {
  91       n1 = n1 - (p1 + (n0 < p0));
  92       n0 = n0 - p0;
  93
  94       q2 = (mp_limb_t) ((TWO64 * (double) n1 + (double) n0) * di);
  95       umul_ppmm (p1, p0, q2, d);
  96
  97       q = q1 + q2;
  98       if (n1 < p1 || (n1 == p1 && n0 < p0))
  99         {
 100           n0 = n0 - p0;
 101           n0 += d;
 102           q--;
 103         }
 104       else
 105         {
 106           n0 = n0 - p0;
 107           if (n0 >= d)
 108             {
 109               n0 -= d;
 110               q++;
 111             }
 112         }
 113     }
 114
 115   *r = n0;
 116   return q;
 117 }