← Back to branch summary

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

« back to all changes in this revision

Viewing changes to gmp3/mpn/generic/dc_divrem_n.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/generic/dc_divrem_n.c
 
1
/* mpn_dc_divrem_n and auxilliary routines.
 
2
 
 
3
   THE FUNCTIONS IN THIS FILE ARE INTERNAL FUNCTIONS WITH MUTABLE
 
4
   INTERFACES.  IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.
 
5
   IN FACT, IT IS ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A
 
6
   FUTURE GNU MP RELEASE.
 
7
 
 
8
 
 
9
Copyright 2000, 2001 Free Software Foundation, Inc.
 
10
Contributed by Paul Zimmermann.
 
11
 
 
12
This file is part of the GNU MP Library.
 
13
 
 
14
The GNU MP Library is free software; you can redistribute it and/or modify
 
15
it under the terms of the GNU Lesser General Public License as published by
 
16
the Free Software Foundation; either version 2.1 of the License, or (at your
 
17
option) any later version.
 
18
 
 
19
The GNU MP Library is distributed in the hope that it will be useful, but
 
20
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 
21
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 
22
License for more details.
 
23
 
 
24
You should have received a copy of the GNU Lesser General Public License
 
25
along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
 
26
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 
27
MA 02111-1307, USA. */
 
28
 
 
29
#include "gmp.h"
 
30
#include "gmp-impl.h"
 
31
 
 
32
/*
 
33
[1] Fast Recursive Division, by Christoph Burnikel and Joachim Ziegler,
 
34
    Technical report MPI-I-98-1-022, october 1998.
 
35
    http://www.mpi-sb.mpg.de/~ziegler/TechRep.ps.gz
 
36
*/
 
37
 
 
38
static mp_limb_t mpn_dc_div_3_halves_by_2
 
39
  _PROTO ((mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n));
 
40
 
 
41
 
 
42
#if 0
 
43
static
 
44
unused_mpn_divrem (qp, qxn, np, nn, dp, dn)
 
45
     mp_ptr qp;
 
46
     mp_size_t qxn;
 
47
     mp_ptr np;
 
48
     mp_size_t nn;
 
49
     mp_srcptr dp;
 
50
     mp_size_t dn;
 
51
{
 
52
  /* This might be useful: */
 
53
  if (qxn != 0)
 
54
    {
 
55
      mp_limb_t c;
 
56
      mp_ptr tp = alloca ((nn + qxn) * BYTES_PER_MP_LIMB);
 
57
      MPN_COPY (tp + qxn - nn, np, nn);
 
58
      MPN_ZERO (tp, qxn);
 
59
      c = mpn_divrem (qp, 0L, tp, nn + qxn, dp, dn);
 
60
      /* Maybe copy proper part of tp to np?  Documentation is unclear about
 
61
         the returned np value when qxn != 0 */
 
62
      return c;
 
63
    }
 
64
}
 
65
#endif
 
66
 
 
67
 
 
68
/* mpn_dc_divrem_n - Implements algorithm of page 8 in [1]: divides (np,2n)
 
69
   by (dp,n) and puts the quotient in (qp,n), the remainder in (np,n).
 
70
   Returns most significant limb of the quotient, which is 0 or 1.
 
71
   Requires that the most significant bit of the divisor is set.  */
 
72
 
 
73
mp_limb_t
 
74
mpn_dc_divrem_n (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n)
 
75
{
 
76
  mp_limb_t qhl, cc;
 
77
 
 
78
  if (n % 2 != 0)
 
79
    {
 
80
      qhl = mpn_dc_divrem_n (qp + 1, np + 2, dp + 1, n - 1);
 
81
      cc = mpn_submul_1 (np + 1, qp + 1, n - 1, dp[0]);
 
82
      cc = mpn_sub_1 (np + n, np + n, 1, cc);
 
83
      if (qhl) cc += mpn_sub_1 (np + n, np + n, 1, dp[0]);
 
84
      while (cc)
 
85
        {
 
86
          qhl -= mpn_sub_1 (qp + 1, qp + 1, n - 1, (mp_limb_t) 1);
 
87
          cc -= mpn_add_n (np + 1, np + 1, dp, n);
 
88
        }
 
89
      qhl += mpn_add_1 (qp + 1, qp + 1, n - 1,
 
90
                        mpn_sb_divrem_mn (qp, np, n + 1, dp, n));
 
91
    }
 
92
  else
 
93
    {
 
94
      mp_size_t n2 = n/2;
 
95
      qhl = mpn_dc_div_3_halves_by_2 (qp + n2, np + n2, dp, n2);
 
96
      qhl += mpn_add_1 (qp + n2, qp + n2, n2,
 
97
                        mpn_dc_div_3_halves_by_2 (qp, np, dp, n2));
 
98
    }
 
99
  return qhl;
 
100
}
 
101
 
 
102
 
 
103
/* divides (np, 3n) by (dp, 2n) and puts the quotient in (qp, n),
 
104
   the remainder in (np, 2n) */
 
105
 
 
106
static mp_limb_t
 
107
mpn_dc_div_3_halves_by_2 (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n)
 
108
{
 
109
  mp_size_t twon = n + n; 
 
110
  mp_limb_t qhl, cc;
 
111
  mp_ptr tmp;
 
112
  TMP_DECL (marker);
 
113
 
 
114
  TMP_MARK (marker);
 
115
  if (n < DC_THRESHOLD)
 
116
    qhl = mpn_sb_divrem_mn (qp, np + n, twon, dp + n, n);
 
117
  else
 
118
    qhl = mpn_dc_divrem_n (qp, np + n, dp + n, n);
 
119
  tmp = (mp_ptr) TMP_ALLOC (twon * BYTES_PER_MP_LIMB);
 
120
  mpn_mul_n (tmp, qp, dp, n);
 
121
  cc = mpn_sub_n (np, np, tmp, twon);
 
122
  TMP_FREE (marker);
 
123
  if (qhl) cc += mpn_sub_n (np + n, np + n, dp, n);
 
124
  while (cc)
 
125
    {
 
126
      qhl -= mpn_sub_1 (qp, qp, n, (mp_limb_t) 1);
 
127
      cc -= mpn_add_n (np, np, dp, twon);
 
128
    }
 
129
  return qhl;
 
130
}