← Back to branch summary

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

« back to all changes in this revision

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