Return to divegcd.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

Annotation of OpenXM_contrib/gmp/mpz/divegcd.c, Revision 1.1.1.1

1.1       ohara       1: /* mpz_divexact_gcd -- exact division optimized for GCDs.
                      2:
                      3:    THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE AND ARE ALMOST CERTAIN TO
                      4:    BE SUBJECT TO INCOMPATIBLE CHANGES IN FUTURE GNU MP RELEASES.
                      5:
                      6: Copyright 2000 Free Software Foundation, Inc.
                      7:
                      8: This file is part of the GNU MP Library.
                      9:
                     10: The GNU MP Library is free software; you can redistribute it and/or modify
                     11: it under the terms of the GNU Lesser General Public License as published by
                     12: the Free Software Foundation; either version 2.1 of the License, or (at your
                     13: option) any later version.
                     14:
                     15: The GNU MP Library is distributed in the hope that it will be useful, but
                     16: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
                     17: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
                     18: License for more details.
                     19:
                     20: You should have received a copy of the GNU Lesser General Public License
                     21: along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
                     22: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
                     23: MA 02111-1307, USA.  */
                     24:
                     25: #include "gmp.h"
                     26: #include "gmp-impl.h"
                     27: #include "longlong.h"
                     28:
                     29:
                     30: /* Set q to a/d, expecting d to be from a GCD and therefore usually small.
                     31:
                     32:    The distribution of GCDs of random numbers can be found in Knuth volume 2
                     33:    section 4.5.2 theorem D.
                     34:
                     35:             GCD     chance
                     36:              1       60.7%
                     37:             2^k      20.2%     (1<=k<32)
                     38:            3*2^k      9.0%     (1<=k<32)
                     39:            other     10.1%
                     40:
                     41:    Only the low limb is examined for optimizations, since GCDs bigger than
                     42:    2^32 (or 2^64) will occur very infrequently.
                     43:
                     44:    Future: This could change to an mpn_divexact_gcd, possibly partly
                     45:    inlined, if/when the relevant mpq functions change to an mpn based
                     46:    implementation.  */
                     47:
                     48:
                     49: static void
                     50: mpz_divexact_by3 (mpz_ptr q, mpz_srcptr a)
                     51: {
                     52:   mp_size_t  size = SIZ(a);
                     53:   if (size == 0)
                     54:     {
                     55:       SIZ(q) = 0;
                     56:       return;
                     57:     }
                     58:   else
                     59:     {
                     60:       mp_size_t  abs_size = ABS(size);
                     61:       mp_ptr     qp;
                     62:
                     63:       MPZ_REALLOC (q, abs_size);
                     64:
                     65:       qp = PTR(q);
                     66:       mpn_divexact_by3 (qp, PTR(a), abs_size);
                     67:
                     68:       abs_size -= (qp[abs_size-1] == 0);
                     69:       SIZ(q) = (size>0 ? abs_size : -abs_size);
                     70:     }
                     71: }
                     72:
                     73: void
                     74: mpz_divexact_gcd (mpz_ptr q, mpz_srcptr a, mpz_srcptr d)
                     75: {
                     76:   ASSERT (mpz_sgn (d) > 0);
                     77:
                     78:   if (SIZ(d) == 1)
                     79:     {
                     80:       mp_limb_t  dl = PTR(d)[0];
                     81:       int        twos;
                     82:
                     83:       if (dl == 1)
                     84:         {
                     85:           if (q != a)
                     86:             mpz_set (q, a);
                     87:           return;
                     88:         }
                     89:       if (dl == 3)
                     90:         {
                     91:           mpz_divexact_by3 (q, a);
                     92:           return;
                     93:         }
                     94:
                     95:       count_trailing_zeros (twos, dl);
                     96:       dl >>= twos;
                     97:
                     98:       if (dl == 1)
                     99:         {
                    100:           mpz_tdiv_q_2exp (q, a, twos);
                    101:           return;
                    102:         }
                    103:       if (dl == 3)
                    104:         {
                    105:           mpz_tdiv_q_2exp (q, a, twos);
                    106:           mpz_divexact_by3 (q, q);
                    107:           return;
                    108:         }
                    109:     }
                    110:
                    111:   mpz_divexact (q, a, d);
                    112: }