@@ -218,6 +218,38 @@ STATIC uint mpn_and(mpz_dig_t *idig, const mpz_dig_t *jdig, uint jlen, const mpz
218218 return idig + 1 - oidig ;
219219}
220220
221+ /* computes i = j & -k = j & (~k + 1)
222+ returns number of digits in i
223+ assumes enough memory in i; assumes normalised j, k
224+ can have i, j, k pointing to same memory
225+ */
226+ STATIC uint mpn_and_neg (mpz_dig_t * idig , const mpz_dig_t * jdig , uint jlen , const mpz_dig_t * kdig , uint klen ) {
227+ mpz_dig_t * oidig = idig ;
228+ mpz_dbl_dig_t carry = 1 ;
229+
230+ for (; jlen > 0 && klen > 0 ; -- jlen , -- klen , ++ idig , ++ jdig , ++ kdig ) {
231+ carry += * kdig ^ DIG_MASK ;
232+ * idig = (* jdig & carry ) & DIG_MASK ;
233+ carry >>= DIG_SIZE ;
234+ }
235+
236+ for (; jlen > 0 ; -- jlen , ++ idig , ++ jdig ) {
237+ carry += DIG_MASK ;
238+ * idig = (* jdig & carry ) & DIG_MASK ;
239+ carry >>= DIG_SIZE ;
240+ }
241+
242+ if (carry != 0 ) {
243+ * idig = carry ;
244+ } else {
245+ // remove trailing zeros
246+ for (-- idig ; idig >= oidig && * idig == 0 ; -- idig ) {
247+ }
248+ }
249+
250+ return idig + 1 - oidig ;
251+ }
252+
221253/* computes i = j | k
222254 returns number of digits in i
223255 assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
@@ -896,24 +928,35 @@ void mpz_sub_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
896928 can have dest, lhs, rhs the same
897929*/
898930void mpz_and_inpl (mpz_t * dest , const mpz_t * lhs , const mpz_t * rhs ) {
899- // make sure lhs has the most digits
900- if (lhs -> len < rhs -> len ) {
901- const mpz_t * temp = lhs ;
902- lhs = rhs ;
903- rhs = temp ;
904- }
905-
906931 if (lhs -> neg == rhs -> neg ) {
907- mpz_need_dig (dest , rhs -> len );
908- dest -> len = mpn_and (dest -> dig , lhs -> dig , lhs -> len , rhs -> dig , rhs -> len );
932+ if (lhs -> neg == 0 ) {
933+ // make sure lhs has the most digits
934+ if (lhs -> len < rhs -> len ) {
935+ const mpz_t * temp = lhs ;
936+ lhs = rhs ;
937+ rhs = temp ;
938+ }
939+ // do the and'ing
940+ mpz_need_dig (dest , rhs -> len );
941+ dest -> len = mpn_and (dest -> dig , lhs -> dig , lhs -> len , rhs -> dig , rhs -> len );
942+ dest -> neg = 0 ;
943+ } else {
944+ // TODO both args are negative
945+ assert (0 );
946+ }
909947 } else {
910- mpz_need_dig (dest , lhs -> len );
911- // TODO
912- assert (0 );
913- // dest->len = mpn_and_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
948+ // args have different sign
949+ // make sure lhs is the positive arg
950+ if (rhs -> neg == 0 ) {
951+ const mpz_t * temp = lhs ;
952+ lhs = rhs ;
953+ rhs = temp ;
954+ }
955+ mpz_need_dig (dest , lhs -> len + 1 );
956+ dest -> len = mpn_and_neg (dest -> dig , lhs -> dig , lhs -> len , rhs -> dig , rhs -> len );
957+ assert (dest -> len <= dest -> alloc );
958+ dest -> neg = 0 ;
914959 }
915-
916- dest -> neg = lhs -> neg ;
917960}
918961
919962/* computes dest = lhs | rhs
0 commit comments