[flint] 01/03: backport gmp support from flint git

felix salfelder felix-guest at alioth.debian.org
Sat Aug 31 19:20:04 UTC 2013


This is an automated email from the git hooks/post-receive script.

felix-guest pushed a commit to branch master-gmp
in repository flint.

commit 1bccef7df68bbf6cffb594dca8863a91e7b1aa7e
Author: Felix Salfelder <felix at salfelder.org>
Date:   Sat Aug 31 16:57:37 2013 +0200

    backport gmp support from flint git
---
 fft.h                               |   42 +++++++++++--
 fft/combine_bits.c                  |    4 +-
 fft/convolution.c                   |   51 ++++++++-------
 fft/mul_fft_main.c                  |    2 +-
 fft/mul_mfa_truncate_sqrt2.c        |    2 +-
 fft/mul_truncate_sqrt2.c            |    7 ++-
 fft/mulmod_2expp1.c                 |    7 ++-
 fft/profile/p-mul_fft_main.c        |    4 +-
 fft/split_bits.c                    |    8 +--
 flint.h                             |   24 +++++++
 mpn_extras.h                        |   55 ++++++++++++++--
 mpn_extras/mulmod_2expp1_basecase.c |  117 +++++++++++++++++++++++++++++++++++
 12 files changed, 275 insertions(+), 48 deletions(-)

diff --git a/fft.h b/fft.h
index 587cbd9..f1d383a 100644
--- a/fft.h
+++ b/fft.h
@@ -49,11 +49,43 @@ or implied, of William Hart.
  extern "C" {
 #endif
 
-#if !defined(__MPIR_RELEASE ) || __MPIR_RELEASE < 20600
-#define mpn_sumdiff_n __MPN(sumdiff_n)
-extern
-mp_limb_t mpn_sumdiff_n(mp_ptr, mp_ptr, mp_srcptr, mp_srcptr, mp_size_t);
-#endif
+/* #if !defined(__MPIR_RELEASE ) || __MPIR_RELEASE < 20600
+// #define mpn_sumdiff_n __MPN(sumdiff_n)
+// extern
+// mp_limb_t mpn_sumdiff_n(mp_ptr, mp_ptr, mp_srcptr, mp_srcptr, mp_size_t);
+// #endif
+*/
+
+static __inline__ mp_limb_t
+mpn_sumdiff_n(mp_ptr s, mp_ptr d, mp_srcptr x, mp_srcptr y, mp_size_t n)
+{
+    mp_limb_t ret;
+    mp_ptr t;
+
+    if (n == 0)
+        return 0;
+
+    if ((s == x && d == y) || (s == y && d == x))
+    {
+        t = flint_malloc(n * sizeof(mp_limb_t));
+        ret = mpn_sub_n(t, x, y, n);
+        ret += 2 * mpn_add_n(s, x, y, n);
+        flint_mpn_copyi(d, t, n);
+        flint_free(t);
+        return ret;
+    }
+
+    if (s == x || s == y)
+    {
+        ret = mpn_sub_n(d, x, y, n);
+        ret += 2 * mpn_add_n(s, x, y, n);
+        return ret;
+    }
+
+    ret = 2 * mpn_add_n(s, x, y, n);
+    ret += mpn_sub_n(d, x, y, n);
+    return ret;
+}
 
 #define fft_sumdiff(t, u, r, s, n) \
    (n == 0 ? 0 : mpn_sumdiff_n(t, u, r, s, n))
diff --git a/fft/combine_bits.c b/fft/combine_bits.c
index f9e1fbd..58f84fe 100644
--- a/fft/combine_bits.c
+++ b/fft/combine_bits.c
@@ -33,7 +33,7 @@ or implied, of William Hart.
 #include "flint.h"
 #include "fft.h"
 
-void fft_combine_limbs(mp_limb_t * res, mp_limb_t ** poly, long length, 
+void fft_combine_limbs(mp_limb_t * res, mp_limb_t ** poly, slong length, 
             mp_size_t coeff_limbs, mp_size_t output_limbs, mp_size_t total_limbs)
 {
    mp_size_t skip, i;
@@ -51,7 +51,7 @@ void fft_combine_limbs(mp_limb_t * res, mp_limb_t ** poly, long length,
    }  
 }
 
-void fft_combine_bits(mp_limb_t * res, mp_limb_t ** poly, long length, 
+void fft_combine_bits(mp_limb_t * res, mp_limb_t ** poly, slong length, 
                   mp_bitcnt_t bits, mp_size_t output_limbs, mp_size_t total_limbs)
 {
    mp_bitcnt_t shift_bits, top_bits = ((FLINT_BITS - 1) & bits);
diff --git a/fft/convolution.c b/fft/convolution.c
index dff9279..907838d 100644
--- a/fft/convolution.c
+++ b/fft/convolution.c
@@ -1,27 +1,32 @@
-/*=============================================================================
+/*
 
-    This file is part of FLINT.
+Copyright 2008-2011 William Hart. All rights reserved.
 
-    FLINT is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
+Redistribution and use in source and binary forms, with or without modification, are
+permitted provided that the following conditions are met:
 
-    FLINT is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
+   1. Redistributions of source code must retain the above copyright notice, this list of
+      conditions and the following disclaimer.
 
-    You should have received a copy of the GNU General Public License
-    along with FLINT; if not, write to the Free Software
-    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
+   2. Redistributions in binary form must reproduce the above copyright notice, this list
+      of conditions and the following disclaimer in the documentation and/or other materials
+      provided with the distribution.
 
-=============================================================================*/
-/******************************************************************************
+THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED
+WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-    Copyright (C) 2008-2011 William Hart
-    
-******************************************************************************/
+The views and conclusions contained in the software and documentation are those of the
+authors and should not be interpreted as representing official policies, either expressed
+or implied, of William Hart.
+
+*/
 
 #include <stdlib.h>
 #include <mpir.h>
@@ -31,13 +36,13 @@
 #include "fmpz_poly.h"
 #include "fft.h"
 
-void fft_convolution(mp_limb_t ** ii, mp_limb_t ** jj, long depth, 
-                              long limbs, long trunc, mp_limb_t ** t1, 
+void fft_convolution(mp_limb_t ** ii, mp_limb_t ** jj, slong depth, 
+                              slong limbs, slong trunc, mp_limb_t ** t1, 
                           mp_limb_t ** t2, mp_limb_t ** s1, mp_limb_t * tt)
 {
-   long n = (1L<<depth), j;
-   long w = (limbs*FLINT_BITS)/n;
-   long sqrt = (1L<<(depth/2));
+   slong n = (1L<<depth), j;
+   slong w = (limbs*FLINT_BITS)/n;
+   slong sqrt = (1L<<(depth/2));
    
    if (depth <= 6)
    {
diff --git a/fft/mul_fft_main.c b/fft/mul_fft_main.c
index 3487f8f..e9970c3 100644
--- a/fft/mul_fft_main.c
+++ b/fft/mul_fft_main.c
@@ -37,7 +37,7 @@ or implied, of William Hart.
 static int fft_tuning_table[5][2] = FFT_TAB;
 
 void flint_mpn_mul_fft_main(mp_limb_t * r1, mp_limb_t * i1, mp_size_t n1, 
-                        mp_limb_t * i2, mp_size_t n2)
+                                                    mp_limb_t * i2, mp_size_t n2)
 {
    mp_size_t off, depth = 6;
    mp_size_t w = 1;
diff --git a/fft/mul_mfa_truncate_sqrt2.c b/fft/mul_mfa_truncate_sqrt2.c
index e820a49..7197e20 100644
--- a/fft/mul_mfa_truncate_sqrt2.c
+++ b/fft/mul_mfa_truncate_sqrt2.c
@@ -34,7 +34,7 @@ or implied, of William Hart.
 #include "ulong_extras.h"
 
 void mul_mfa_truncate_sqrt2(mp_limb_t * r1, mp_limb_t * i1, mp_size_t n1, 
-                        mp_limb_t * i2, mp_size_t n2, mp_bitcnt_t depth, mp_bitcnt_t w)
+                  mp_limb_t * i2, mp_size_t n2, mp_bitcnt_t depth, mp_bitcnt_t w)
 {
    mp_size_t n = (1UL<<depth);
    mp_bitcnt_t bits1 = (n*w - (depth+1))/2; 
diff --git a/fft/mul_truncate_sqrt2.c b/fft/mul_truncate_sqrt2.c
index c3f6ef9..bdff77d 100644
--- a/fft/mul_truncate_sqrt2.c
+++ b/fft/mul_truncate_sqrt2.c
@@ -31,9 +31,10 @@ or implied, of William Hart.
 #include "mpir.h"
 #include "flint.h"
 #include "fft.h"
-      
+#include "mpn_extras.h"
+
 void mul_truncate_sqrt2(mp_limb_t * r1, mp_limb_t * i1, mp_size_t n1, 
-                        mp_limb_t * i2, mp_size_t n2, mp_bitcnt_t depth, mp_bitcnt_t w)
+                  mp_limb_t * i2, mp_size_t n2, mp_bitcnt_t depth, mp_bitcnt_t w)
 {
    mp_size_t n = (1UL<<depth);
    mp_bitcnt_t bits1 = (n*w - (depth+1))/2; 
@@ -93,7 +94,7 @@ void mul_truncate_sqrt2(mp_limb_t * r1, mp_limb_t * i1, mp_size_t n1,
       mpn_normmod_2expp1(ii[j], limbs);
       if (i1 != i2) mpn_normmod_2expp1(jj[j], limbs);
       c = 2*ii[j][limbs] + jj[j][limbs];
-      ii[j][limbs] = mpn_mulmod_2expp1(ii[j], ii[j], jj[j], c, n*w, tt);
+      ii[j][limbs] = flint_mpn_mulmod_2expp1_basecase(ii[j], ii[j], jj[j], c, n*w, tt);
    }
 
    ifft_truncate_sqrt2(ii, n, w, &t1, &t2, &s1, trunc);
diff --git a/fft/mulmod_2expp1.c b/fft/mulmod_2expp1.c
index 466adf5..0e0dd5c 100644
--- a/fft/mulmod_2expp1.c
+++ b/fft/mulmod_2expp1.c
@@ -34,6 +34,7 @@ or implied, of William Hart.
 #include "longlong.h"
 #include "ulong_extras.h"
 #include "fft_tuning.h"
+#include "mpn_extras.h"
 
 static mp_size_t mulmod_2expp1_table_n[FFT_N_NUM] = MULMOD_TAB;
 
@@ -121,7 +122,7 @@ void _fft_mulmod_2expp1(mp_limb_t * r1, mp_limb_t * i1, mp_limb_t * i2,
    {
       if (i1 != i2) mpn_normmod_2expp1(jj[j], limbs);
       c = 2*ii[j][limbs] + jj[j][limbs];
-      ii[j][limbs] = mpn_mulmod_2expp1(ii[j], ii[j], jj[j], c, n*w, tt);
+      ii[j][limbs] = flint_mpn_mulmod_2expp1_basecase(ii[j], ii[j], jj[j], c, n*w, tt);
    }
    
    ifft_negacyclic(ii, n, w, &t1, &t2, &s1);
@@ -203,7 +204,7 @@ void fft_mulmod_2expp1(mp_limb_t * r, mp_limb_t * i1, mp_limb_t * i2,
 
    if (limbs <= FFT_MULMOD_2EXPP1_CUTOFF) 
    {
-      r[limbs] = mpn_mulmod_2expp1(r, i1, i2, c, bits, tt);
+      r[limbs] = flint_mpn_mulmod_2expp1_basecase(r, i1, i2, c, bits, tt);
       return;
    }
    
@@ -218,7 +219,7 @@ void fft_mulmod_2expp1(mp_limb_t * r, mp_limb_t * i1, mp_limb_t * i2,
    _fft_mulmod_2expp1(r, i1, i2, limbs, depth1, w1);
 }
 
-long fft_adjust_limbs(mp_size_t limbs)
+slong fft_adjust_limbs(mp_size_t limbs)
 {
    mp_size_t bits1 = limbs*FLINT_BITS, bits2;
    mp_size_t depth = 1, limbs2, depth1 = 1, depth2 = 1, adj;
diff --git a/fft/profile/p-mul_fft_main.c b/fft/profile/p-mul_fft_main.c
index fa89d3c..c635063 100644
--- a/fft/profile/p-mul_fft_main.c
+++ b/fft/profile/p-mul_fft_main.c
@@ -62,8 +62,8 @@ main(void)
        r1 = i2 + int_limbs;
        r2 = r1 + 2*int_limbs;
    
-       mpn_urandomb(i1, state->gmp_state, int_limbs*FLINT_BITS);
-       mpn_urandomb(i2, state->gmp_state, int_limbs*FLINT_BITS);
+       flint_mpn_urandomb(i1, state->gmp_state, int_limbs*FLINT_BITS);
+       flint_mpn_urandomb(i2, state->gmp_state, int_limbs*FLINT_BITS);
   
        for (j = 0; j < iters; j++)
           //mpn_mul(r2, i1, int_limbs, i2, int_limbs);
diff --git a/fft/split_bits.c b/fft/split_bits.c
index 282adcd..46f1463 100644
--- a/fft/split_bits.c
+++ b/fft/split_bits.c
@@ -28,12 +28,12 @@ or implied, of William Hart.
 
 */
 
-#include "mpir.h"
+#include "gmp.h"
 #include "flint.h"
 #include "fft.h"
 
 mp_size_t fft_split_limbs(mp_limb_t ** poly, mp_limb_t * limbs, 
-                mp_size_t total_limbs, mp_size_t coeff_limbs, mp_size_t output_limbs)
+            mp_size_t total_limbs, mp_size_t coeff_limbs, mp_size_t output_limbs)
 {
    mp_size_t i, skip, length = (total_limbs - 1)/coeff_limbs + 1;
    
@@ -53,11 +53,11 @@ mp_size_t fft_split_limbs(mp_limb_t ** poly, mp_limb_t * limbs,
 }
 
 mp_size_t fft_split_bits(mp_limb_t ** poly, mp_limb_t * limbs, 
-               mp_size_t total_limbs, mp_bitcnt_t bits, mp_size_t output_limbs)
+                 mp_size_t total_limbs, mp_bitcnt_t bits, mp_size_t output_limbs)
 {
    mp_size_t i, coeff_limbs, limbs_left, length = (FLINT_BITS*total_limbs - 1)/bits + 1;
    mp_bitcnt_t shift_bits, top_bits = ((FLINT_BITS - 1) & bits);
-   mp_limb_t * limb_ptr;
+   mp_srcptr limb_ptr;
    mp_limb_t mask;
    
    if (top_bits == 0)
diff --git a/flint.h b/flint.h
index 92ce429..43c4c60 100644
--- a/flint.h
+++ b/flint.h
@@ -31,6 +31,8 @@
 #include "longlong.h"
 #include "config.h"
 
+typedef mp_size_t slong;
+
 #ifdef __cplusplus
  extern "C" {
 #endif
@@ -177,6 +179,28 @@ unsigned int FLINT_BIT_COUNT(mp_limb_t x)
          (xxx)[ixxx] = yyy; \
    } while (0)
 
+/* compatibility between gmp and mpir */
+#ifndef mpn_com_n
+#define mpn_com_n mpn_com
+#endif
+
+#ifndef mpn_neg_n
+#define mpn_neg_n mpn_neg
+#endif
+
+#ifndef mpn_tdiv_q
+/* substitute for mpir's mpn_tdiv_q */
+static __inline__
+void mpn_tdiv_q(mp_ptr qp,
+	   mp_srcptr np, mp_size_t nn,
+	   mp_srcptr dp, mp_size_t dn)
+    {
+    mp_ptr _scratch = (mp_ptr) flint_malloc(dn * sizeof(mp_limb_t));
+    mpn_tdiv_qr(qp, _scratch, 0, np, nn, dp, dn);
+    flint_free(_scratch);
+    }
+#endif
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/mpn_extras.h b/mpn_extras.h
index 5a75af5..c3b2a4d 100644
--- a/mpn_extras.h
+++ b/mpn_extras.h
@@ -53,7 +53,7 @@
         (bn) = __tn;           \
     } while (0)
 
-/* Not defined in mpir.h
+/* Not defined in gmp.h
 mp_limb_t  __gmpn_modexact_1_odd(mp_srcptr src, mp_size_t size,
                                  mp_limb_t divisor);
 #define mpn_modexact_1_odd __gmpn_modexact_1_odd
@@ -75,7 +75,7 @@ flint_mpn_divisible_1_p(mp_srcptr x, mp_size_t xsize, mp_limb_t d)
 static __inline__
 int flint_mpn_zero_p(mp_srcptr x, mp_size_t xsize)
 {
-    long i;
+    slong i;
     for (i = 0; i < xsize; i++)
     {
         if (x[i])
@@ -100,7 +100,7 @@ mp_size_t flint_mpn_remove_2exp(mp_ptr x, mp_size_t xsize, mp_bitcnt_t *bits);
 mp_size_t flint_mpn_remove_power_ascending(mp_ptr x, mp_size_t xsize,
                                      mp_ptr p, mp_size_t psize, ulong *exp);
 
-int flint_mpn_factor_trial(mp_srcptr x, mp_size_t xsize, long start, long stop);
+int flint_mpn_factor_trial(mp_srcptr x, mp_size_t xsize, slong start, slong stop);
 
 int flint_mpn_divides(mp_ptr q, mp_srcptr array1, 
          mp_size_t limbs1, mp_srcptr arrayg, mp_size_t limbsg, mp_ptr temp);
@@ -182,7 +182,54 @@ do {                                                         \
 void
 flint_mpn_harmonic_odd_balanced(mp_ptr t, mp_size_t * tsize,
                           mp_ptr v, mp_size_t * vsize,
-                          long a, long b, long n, int d);
+                          slong a, slong b, slong n, int d);
+
+mp_limb_t flint_mpn_preinv1(mp_limb_t d, mp_limb_t d2);
+
+mp_limb_t flint_mpn_divrem_preinv1(mp_ptr q, mp_ptr a, 
+           mp_size_t m, mp_srcptr b, mp_size_t n, mp_limb_t dinv);
+
+#define flint_mpn_divrem21_preinv(q, a_hi, a_lo, dinv) \
+   do { \
+      mp_limb_t __q2, __q3, __q4; \
+      umul_ppmm((q), __q2, (a_hi), (dinv)); \
+      umul_ppmm(__q3, __q4, (a_lo), (dinv)); \
+      add_ssaaaa((q), __q2, (q), __q2, 0, __q3); \
+      add_ssaaaa((q), __q2, (q), __q2, (a_hi), (a_lo)); \
+   } while (0)
+
+void flint_mpn_mulmod_preinv1(mp_ptr r, 
+        mp_srcptr a, mp_srcptr b, mp_size_t n, 
+        mp_srcptr d, mp_limb_t dinv, ulong norm);
+
+void flint_mpn_preinvn(mp_ptr dinv, mp_srcptr d, mp_size_t n);
+
+void flint_mpn_mulmod_preinvn(mp_ptr r, 
+        mp_srcptr a, mp_srcptr b, mp_size_t n, 
+        mp_srcptr d, mp_srcptr dinv, ulong norm);
+
+int flint_mpn_mulmod_2expp1_basecase(mp_ptr xp, mp_srcptr yp, mp_srcptr zp, int c,
+    mp_bitcnt_t b, mp_ptr tp);
+
+static __inline__
+void flint_mpn_rrandom(mp_limb_t *rp, gmp_randstate_t state, mp_size_t n)
+{
+  __mpz_struct str;
+  str._mp_d = rp;
+  str._mp_alloc = n;
+  str._mp_size =n;
+  mpz_rrandomb(&str,state,FLINT_BITS*n);
+}
+
+static __inline__
+void flint_mpn_urandomb(mp_limb_t *rp, gmp_randstate_t state, mp_bitcnt_t n)
+{
+  __mpz_struct str;
+  str._mp_d = rp;
+  str._mp_alloc = (n + FLINT_BITS - 1)/FLINT_BITS;
+  str._mp_size = (n + FLINT_BITS - 1)/FLINT_BITS;
+  mpz_rrandomb(&str,state,n);
+}
 
 #ifdef __cplusplus
 }
diff --git a/mpn_extras/mulmod_2expp1_basecase.c b/mpn_extras/mulmod_2expp1_basecase.c
new file mode 100644
index 0000000..844bb1a
--- /dev/null
+++ b/mpn_extras/mulmod_2expp1_basecase.c
@@ -0,0 +1,117 @@
+/*=============================================================================
+
+    This file is part of FLINT.
+
+    FLINT is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    FLINT is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with FLINT; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright 2009 Jason Moxham
+
+******************************************************************************/
+
+#include "gmp.h"
+#include "flint.h"
+#include "mpn_extras.h"
+
+#define BITS_TO_LIMBS(b) (((b) + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS)
+
+/* ret+(xp,n)=(yp,n)*(zp,n) % 2^b+1  
+   needs (tp,2n) temp space , everything reduced mod 2^b 
+   inputs,outputs are fully reduced
+   NOTE: 2n is not the same as 2b rounded up to nearest limb
+*/
+static __inline__ int
+mpn_mulmod_2expp1_internal(mp_ptr xp, mp_srcptr yp, mp_srcptr zp,
+    mp_bitcnt_t b, mp_ptr tp)
+{
+    mp_size_t n, k;
+    mp_limb_t c;
+
+    n = BITS_TO_LIMBS(b);
+    k = GMP_NUMB_BITS * n - b;
+
+    mpn_mul_n(tp, yp, zp, n);
+
+    if (k == 0)
+    {
+        c = mpn_sub_n(xp, tp, tp + n, n);
+        return mpn_add_1 (xp, xp, n, c);
+    }
+
+    c = tp[n - 1];
+    tp[n - 1] &= GMP_NUMB_MASK >> k;
+
+#if HAVE_NATIVE_mpn_sublsh_nc
+    c = mpn_sublsh_nc (xp, tp, tp + n, n, k, c);
+#else
+    {
+        mp_limb_t c1;
+        c1 = mpn_lshift (tp + n, tp + n, n, k);
+        tp[n] |= c >> (GMP_NUMB_BITS - k);
+        c = mpn_sub_n (xp, tp, tp + n, n) + c1;
+    }
+#endif
+    c = mpn_add_1 (xp, xp, n, c);
+    xp[n - 1] &= GMP_NUMB_MASK >> k;
+    return c;
+}
+
+/* c is the top bits of the inputs, must be fully reduced */
+int
+flint_mpn_mulmod_2expp1_basecase (mp_ptr xp, mp_srcptr yp, mp_srcptr zp, int c,
+    mp_bitcnt_t b, mp_ptr tp)
+{
+    int cy, cz;
+    mp_size_t n, k;
+
+    cy = c & 2;
+    cz = c & 1;
+    n = BITS_TO_LIMBS(b);
+    k = GMP_NUMB_BITS * n - b;
+
+    if (cy == 0)
+    {
+        if (cz == 0)
+        {
+            c = mpn_mulmod_2expp1_internal(xp, yp, zp, b, tp);
+        }
+        else
+        {
+            c = mpn_neg_n(xp, yp, n);
+            c = mpn_add_1 (xp, xp, n, c);
+            xp[n - 1] &= GMP_NUMB_MASK >> k;
+        }
+    }
+    else
+    {
+        if (cz == 0)
+	    {
+            c = mpn_neg_n(xp, zp, n);
+            c = mpn_add_1(xp, xp, n, c);
+            xp[n - 1] &= GMP_NUMB_MASK >> k;
+        }
+        else
+        {
+            c = 0;
+            xp[0] = 1;
+            flint_mpn_zero(xp + 1, n - 1);
+        }
+    }
+
+    return c;
+}
+

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/flint.git



More information about the debian-science-commits mailing list