Return to karadiv.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpfr

Annotation of OpenXM_contrib/gmp/mpfr/karadiv.c, Revision 1.1.1.1

1.1       maekawa     1: /*  mpn_divrem_n -- Karatsuba division in 2*K(n) limb operations
                      2:
                      3: Copyright (C) 1999-2000 PolKA project, Inria Lorraine and Loria
                      4:
                      5: This file is part of the MPFR Library.
                      6:
                      7: The MPFR Library is free software; you can redistribute it and/or modify
                      8: it under the terms of the GNU Library General Public License as published by
                      9: the Free Software Foundation; either version 2 of the License, or (at your
                     10: option) any later version.
                     11:
                     12: The MPFR Library is distributed in the hope that it will be useful, but
                     13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
                     14: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
                     15: License for more details.
                     16:
                     17: You should have received a copy of the GNU Library General Public License
                     18: along with the MPFR Library; see the file COPYING.LIB.  If not, write to
                     19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
                     20: MA 02111-1307, USA. */
                     21:
                     22: /*
                     23: [1] Fast Recursive Division, by Christoph Burnikel and Joachim Ziegler,
                     24:     Technical report MPI-I-98-1-022, october 1998,
                     25:     cf http://www.mpi-sb.mpg.de/~ziegler/TechRep.ps.gz.
                     26:     Implemented by Paul Zimmermann, 1999.
                     27: */
                     28:
                     29: #include "gmp.h"
                     30: #include "gmp-impl.h"
                     31: #include "mpfr.h"
                     32:
                     33: extern void mpn_divrem_n2 (mp_limb_t *, mp_limb_t *, mp_limb_t *, mp_size_t, mp_limb_t *);
                     34: extern void mpn_divrem_3by2 (mp_limb_t *, mp_limb_t *, mp_limb_t *, mp_size_t, mp_limb_t *);
                     35:
                     36: /* mpn_divrem_n(n) calls 2*mul(n/2)+2*div(n/2), thus to be faster then
                     37:    div(n)=4*div(n/2), we need mul(n/2) to be faster than the classic way,
                     38:    i.e. n/2 >= KARATSUBA_MUL_THRESHOLD */
                     39:
                     40: #define DIV_LIMIT (7*KARATSUBA_MUL_THRESHOLD)
                     41:
                     42: static void mpn_decr(mp_limb_t *Q)
                     43: {
                     44:   while ((*Q++)-- == 0);
                     45: }
                     46:
                     47: /* implements algorithm of page 8 in [1]: divides (A,2n) by (B,n) and puts the
                     48:    quotient in (Q,n), the remainder in (A,n).
                     49:    Returns most significant limb of the quotient, which is 0 or 1.
                     50: */
                     51: mp_limb_t
                     52: mpn_divrem_n(mp_limb_t *Q, mp_limb_t *A, mp_limb_t *B, mp_size_t n)
                     53: {
                     54:   if (n<DIV_LIMIT) return mpn_divrem(Q, 0, A, 2*n, B, n);
                     55:   else {
                     56:     mp_limb_t cc=0;
                     57:     if (mpn_cmp(A+n, B, n)>=0) {
                     58:       cc=1;
                     59:       mpn_sub_n(A+n, A+n, B, n);
                     60:     }
                     61:     if (n%2) {
                     62:        /* divide (2n-2) most significant limbs from A by those (n-1) from B */
                     63:        mpn_divrem_n(Q+1, A+2, B+1, n-1);
                     64:        /* now (Q+1, n-1) contains the quotient of (A+2,2n-2) by (B+1,n-1)
                     65:          and (A+2, n-1) contains the remainder */
                     66:        if (mpn_sub_1(A+n, A+n, 1, mpn_submul_1(A+1, Q+1, n-1, B[0]))) {
                     67:         /* quotient two large */
                     68:         mpn_decr(Q+1);
                     69:         if (mpn_add_n(A+1, A+1, B, n)==0) {
                     70:           mpn_decr(Q+1); mpn_add_n(A+1, A+1, B, n);
                     71:         }
                     72:        }
                     73:        /* now divide (A,n+1) by (B,n) */
                     74:        mpn_divrem(Q, 0, A, n+1, B, n);
                     75:     }
                     76:     else {
                     77:        mp_limb_t *tmp; int n2=n/2;
                     78:        TMP_DECL (marker);
                     79:
                     80:        TMP_MARK (marker);
                     81:        tmp = (mp_limb_t*) TMP_ALLOC(n*sizeof(mp_limb_t));
                     82:        mpn_divrem_3by2(Q+n2, A+n2, B, n2, tmp);
                     83:        mpn_divrem_3by2(Q, A, B, n2, tmp);
                     84:        TMP_FREE (marker);
                     85:     }
                     86:     return cc;
                     87:   }
                     88: }
                     89:
                     90: /* inner procedure, with no memory allocation
                     91:    assumes mpn_cmp(A+n, B, n) < 0
                     92: */
                     93: void mpn_divrem_n2(mp_limb_t *Q, mp_limb_t *A, mp_limb_t *B, mp_size_t n,
                     94:                        mp_limb_t *tmp)
                     95: {
                     96:   if (n%2) {
                     97:     /* divide (2n-2) most significant limbs from A by those (n-1) from B */
                     98:     mpn_divrem_n2(Q+1, A+2, B+1, n-1, tmp);
                     99:     /* now (Q+1, n-1) contains the quotient of (A+2,2n-2) by (B+1,n-1)
                    100:        and (A+2, n-1) contains the remainder */
                    101:     if (mpn_sub_1(A+n, A+n, 1, mpn_submul_1(A+1, Q+1, n-1, B[0]))) {
                    102:       /* quotient two large */
                    103:       mpn_decr(Q+1);
                    104:       if (mpn_add_n(A+1, A+1, B, n)==0) { /* still too large */
                    105:        mpn_decr(Q+1); mpn_add_n(A+1, A+1, B, n);
                    106:       }
                    107:     }
                    108:     /* now divide (A,n+1) by (B,n) */
                    109:     mpn_divrem(Q, 0, A, n+1, B, n);
                    110:   }
                    111:   else {
                    112:     int n2=n/2;
                    113:     mpn_divrem_3by2(Q+n2, A+n2, B, n2, tmp);
                    114:     mpn_divrem_3by2(Q, A, B, n2, tmp);
                    115:   }
                    116: }
                    117:
                    118: /* divides (A,3n) by (B,2n) and puts the quotient in (Q,n),
                    119:    the remainder in (A,2n) */
                    120: void
                    121: mpn_divrem_3by2(mp_limb_t *Q, mp_limb_t *A, mp_limb_t *B,
                    122:                             mp_size_t n, mp_limb_t *tmp)
                    123: {
                    124:   int twon = n+n;
                    125:
                    126:   if (n<DIV_LIMIT) mpn_divrem(Q, 0, A+n, twon, B+n, n);
                    127:   else mpn_divrem_n2(Q, A+n, B+n, n, tmp);
                    128:   /* q=(Q,n), c=(A+n,n) with the notations of [1] */
                    129:   mpn_mul_n(tmp, Q, B, n);
                    130:   if (mpn_sub_n(A, A, tmp, twon)) /* R=(A,2n) */ {
                    131:     mpn_decr(Q);
                    132:     if (mpn_add_n(A, A, B, twon)==0) { /* Q still too large */
                    133:       mpn_decr(Q); mpn_add_n(A, A, B, twon);
                    134:     }
                    135:   }
                    136:   return;
                    137: }