← Back to branch summary

~ubuntu-branches/ubuntu/quantal/gclcvs/quantal

« back to all changes in this revision

Viewing changes to gmp3/mpn/pa64w/udiv_qrnnd.c

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2004-06-24 15:13:46 UTC
  • Revision ID: james.westby@ubuntu.com-20040624151346-xh0xaaktyyp7aorc
Tags: 2.7.0-26
C_GC_OFFSET is 2 on m68k-linux

Show diffs side-by-side

added added

removed removed

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