Re: C99 and GMP

[Marc Glisse]

On Fri, 27 Apr 2018, paul zimmermann wrote:


https://gmplib.org/devel/lcov/shell/gmp/mini-gmp/mini-mpq.c.gcov.html


quite interesting. Why is gmp/mpn not tested in the head coverage?


It is tested. It appears as /var/tmp/lcov/gmp/mpn because it is a set of 
symlinks created at build time.


--
Marc Glisse
___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: dead code in div_q.c?

[Niels Möller]

paul zimmermann  writes:

>> It would be good to either prove that qh > 0 can't happen, or add a test
>> that will exercise the code handling that case.
>
> yes more generally it would be nice to know for the divappr functions, if the
> exact quotient fits on qn limbs (i.e., the exact qh is 0), can the
> "approximate" qh be 1?

I suspect we may have qh == 0 always in this special context, where we
have additional constraints on N and D.

But in general, I doubt we can promise that the qh return value from
divappr euals what we would get for the correct quotient. I.e., that the
qh return values from div and divappr always match. Seems particularly
unlikely for the mu divappr, which can have an error larger than 1 unit.

Not sure how to come up with an example, though.

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.
___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: C99 and GMP

[paul zimmermann]

> > quite interesting. Why is gmp/mpn not tested in the head coverage?
> 
> It is tested. It appears as /var/tmp/lcov/gmp/mpn because it is a set of 
> symlinks created at build time.

sorry I missed that. I see some of the files are not tested at all
(add_err3_n.c for example), and some have a low coverage (div_qr_1.c
for example). Are there any plans to improve that?

Paul
___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: documentation on internals not up to date

[Marco Bodrato]

from gmp-bugs.

Il Ven, 27 Aprile 2018 7:21 am, Niels Möller ha scritto:
> "Marco Bodrato"  writes:

>> Currently there's always at least one readable limb, so for instance
>> @code{mpz_get_ui} can fetch @code{_mp_d[0]} unconditionally (though its
>> value is then only wanted if @code{_mp_size} is non-zero).
>
> Maybe change "is only wanted if" to "is valid only if" or "is usable
> only if". At some point, I suffered the misunderstanding that _mp_size
> == 0 implies _mp_d[0] == 0, which isn't right.

Do you mean that we should underline that also with _mp_size == 0, any
limb above the size of the number is undefined?

> _mp_alloc == 0 and _mp_size != 0 is a read-only value, _mp_d is neither
> written, reallocated or freed by mpz functions. It must not be passed as
> destination argument to any mpz function. Should also link to docs for
> mpz_roinit_n and MPZ_ROINIT_N.

Currently, if an mpz_t is initialised with _roinit, it can be passed to
_clear or _clears with no errors. Should we document this? I think we
should.

Moreover, the various mpz_set_ functions should work smoothly too. And I'm
quite sure that any mpz function is able to "overwrite" an mpz_t that was
initialised with _roinit. I mean, the following code actually works:

int main(void)
{
  mpz_t f0, f1, fn;
  int n = 10;
  mp_limb_t dummy = 1;

  mpz_roinit_n (f0, , 1);
  mpz_roinit_n (f1, , 1);
  mpz_roinit_n (fn, , 0);

  for(;--n;) {
mpz_add (fn, f0, f1);
mpz_swap (f0, f1);
mpz_swap (f1, fn);
  }
  gmp_printf ("%Zd\n", f1);
  mpz_clears (f0, f1, fn, NULL);
}

Problems arise only when ones try to _update_ a value with an mpz function.
 mpz_add (f0, f0, f1); mpz_swap (f0, f1); does not work...

Should we document this too. I think we should not.

Ĝis,
m

-- 
http://bodrato.it/papers/

___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: C99 and GMP

[Marc Glisse]

On Fri, 27 Apr 2018, Marc Glisse wrote:


On Fri, 27 Apr 2018, paul zimmermann wrote:


quite interesting. Why is gmp/mpn not tested in the head coverage?


It is tested. It appears as /var/tmp/lcov/gmp/mpn because it is a set of
symlinks created at build time.


sorry I missed that. I see some of the files are not tested at all
(add_err3_n.c for example), and some have a low coverage (div_qr_1.c
for example). Are there any plans to improve that?


Even for the generic files, the coverage varies per target, maybe add_err3_n 
is better tested on ARM or something (not really answering your question).


Hmm, no, mpn_add_err3_n really seems completely unused. There is a refmpn 
implementation ready for comparisons...


--
Marc Glisse
___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: dead code in div_q.c?

[paul zimmermann]

   Niels,

I fully agree with your analysis. Small comments below.

> From: ni...@lysator.liu.se (Niels Möller)
> Date: Thu, 26 Apr 2018 22:12:09 +0200
> 
> paul zimmermann  writes:
> 
> > In the second case (lines 267-276) this is less obvious since you use an
> > approximate division (divappr), and I don't know what is the maximal allowed
> > difference between divappr and div. If divappr always return a quotient
> > less or equal to the true quotient, then qh cannot be zero in this case.
> > But if qh can be larger than the true quotient, even by one bit, then it
> > is possible that qh is 1, in the case {new_np, new_nn} = 0111...111 and
> > {new_dp, dn} = 1000...000.
> 
> You're right again... I thought we had more margin.
> 
> According to comments, quotient from sb and dc divappr is correct or
> one too large, while quotient from mu divappr can be at most 4 too
> large.
> 
> And the limit case, N = B^n/2 - 1, D = B^k/2 is 
> 
>floor ((B^n/2 - 1) / (B^k/2)) = floor(B^{n-k} - B^{-k}/2) 
>  = B^{n-k} - 1 
> 
> So even we only a single unit of error, we could get qh == 1.
> 
> I wonder what actual error can be produced by the divappr functions in
> the corner cases that error could spill over into qh?
> 
> In the cy != 0 case, it looks like we we make a balanced divappr call,
> (2qn+2) / (qn+1). Lets put k = qn+1. For divappr we then have
> 
>   N <= B^{2k}/2 - 1, 
>   D >= B^k/2, lets say D = B/2 + e, e >= 0

D = B/2 + e should read D = B^k/2 + e

> To have any possibility to get qh > 0 with from the approximation error,
> we must have Q >= B^k - 4, where Q denotes the correct quotient. So
> assume that is the case. Now,
> 
>   D <= N/Q <= (B^{2k}/2 - 1) / (B^k - 4)
> 
> which implies
> 
>   e <= N/Q - B^k/2 <= (2 B^k - 1) / (B^k - 4) < 3
> 
> So we need to consider only D = B^k/2 + {0, 1, 2}. Then, any inverse not
> taking the very least significant D limb into account (including mu inverse)
> must be *all* ones, right?

yes if you compute the inverse as floor((B^(2k-2)-1)/floor(D/B))

> It would be good to either prove that qh > 0 can't happen, or add a test
> that will exercise the code handling that case.

yes more generally it would be nice to know for the divappr functions, if the
exact quotient fits on qn limbs (i.e., the exact qh is 0), can the
"approximate" qh be 1?

> > beware that in the case cy <> 0, one divides {new_np, nn+1} by {new_dp, dn}
> > thus filling {qp, nn-dn+1}. 
> 
> Thanks for pointing this out. Now I understand better.
> 
> Regards,
> /Niels
> 
> -- 
> Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
> Internet email is subject to wholesale government surveillance.

Paul
___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: documentation on internals not up to date

[Marc Glisse]

On Fri, 27 Apr 2018, Niels Möller wrote:


_mp_alloc == 0 and _mp_size != 0 is a read-only value, _mp_d is neither
written, reallocated or freed by mpz functions. It must not be passed as
destination argument to any mpz function. Should also link to docs for
mpz_roinit_n and MPZ_ROINIT_N.


Currently, if an mpz_t is initialised with _roinit, it can be passed to
_clear or _clears with no errors. Should we document this? I think we
should.

Moreover, the various mpz_set_ functions should work smoothly too.


I'd prefer that we not document any way to pass _roinit values to any
mpz functions taking a non-const mpz_t, even if it happens to work in
the current implementation. Maybe as a later extension, *if* we find
some use cases where it provides a significant advantage.


There would be a significant advantage to mpq if we could have a
non-allocated 1 for the denominator. But indeed, with the current code 
where only some mpz functions would work, it seems safer to document that 
none work.


--
Marc Glisse
___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: C99 and GMP

[Marc Glisse]

On Fri, 27 Apr 2018, paul zimmermann wrote:


quite interesting. Why is gmp/mpn not tested in the head coverage?


It is tested. It appears as /var/tmp/lcov/gmp/mpn because it is a set of
symlinks created at build time.


sorry I missed that. I see some of the files are not tested at all
(add_err3_n.c for example), and some have a low coverage (div_qr_1.c
for example). Are there any plans to improve that?


Even for the generic files, the coverage varies per target, maybe 
add_err3_n is better tested on ARM or something (not really answering your 
question).


--
Marc Glisse
___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: documentation on internals not up to date

[Marc Glisse]

On Fri, 27 Apr 2018, Niels Möller wrote:


Marc Glisse  writes:


There would be a significant advantage to mpq if we could have a
non-allocated 1 for the denominator. But indeed, with the current code
where only some mpz functions would work, it seems safer to document
that none work.


We could do that internally, even if we don't advertise it for other gmp
users.


We document that users can use the denominator of a mpq_t as a regular 
mpz_t and apply pretty much any mpz_t operation to it. So it seems hard to 
handle just that case internally.



Then mpq_init wouldn't do any allocation, right?


Right.

We could have a single mpq object representing 0/1 with _mp_alloc == 0 
for both parts, and initialize with struct assignment or memcpy.


We would just ensure that mpz realloc keeps supporting thuis case.


--
Marc Glisse
___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: documentation on internals not up to date

[Niels Möller]

Marc Glisse  writes:

> There would be a significant advantage to mpq if we could have a
> non-allocated 1 for the denominator. But indeed, with the current code
> where only some mpz functions would work, it seems safer to document
> that none work.

We could do that internally, even if we don't advertise it for other gmp
users. Then mpq_init wouldn't do any allocation, right? We could have a
single mpq object representing 0/1 with _mp_alloc == 0 for both parts,
and initialize with struct assignment or memcpy.

We would just ensure that mpz realloc keeps supporting thuis case.

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.
___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: documentation on internals not up to date

[Marco Bodrato]

Ciao,

Il Ven, 27 Aprile 2018 2:46 pm, Marc Glisse ha scritto:
> On Fri, 27 Apr 2018, Niels Möller wrote:
>> Marc Glisse  writes:
>>> There would be a significant advantage to mpq if we could have a
>>> non-allocated 1 for the denominator. But indeed, with the current code
>>
>> We could do that internally, even if we don't advertise it for other gmp
>> users.
>
> We document that users can use the denominator of a mpq_t as a regular
> mpz_t and apply pretty much any mpz_t operation to it. So it seems hard to
> handle just that case internally.

I agree, mpq_denref(Q) should work as a generic mpz_t. That's why I also
pushed the "Support lazy mpq_t also in the denominator" patch.

>> We could have a single mpq object representing 0/1 with _mp_alloc == 0
>> for both parts, and initialize with struct assignment or memcpy.
>>
>> We would just ensure that mpz realloc keeps supporting thuis case.

... and it's not just a problem of supporting mpq_denref() we should also
change some mpq_ functions. If someone does init an mpq_t with:
  mpz_init (mpq_numref (accumulator));
  mpz_roinit_n (mpq_denref (accumulator), _one, 1);

then a simple, mpq_add (accumulator, accumulator, another_mpq) might fail.

Niels is right when he says we should not mention on the manual that
sometimes writing on a read-only mpz does actually work.
It does not work as a general rule and it may not be supported by future
releases, that's the only information we can document.

Ĝis,
m

-- 
http://bodrato.it/papers/

___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

divappr interface

[Niels Möller]

Hi, 

after the discussion on div_q, I had another look at divappr, or more
precisely, the mpn_sbpi1_divappr_q function.

It's not entirely obvious from the implementation which limbs of the
{np, nn} argument are actually used. But I think it should access only
the qn+2 upper limbs. For accuracy, it should be sufficient with qn+1
limbs, but I suspect qn+2 is more convenient since we lack a variant of
mpn_submul_1 which ignores the low half of the lowest limb product (we'd
need a function for R -= floor (D * q / B)).

For a call to divappr to produce qn quotient limbs (plus possibly a qh
being zero or one), inputs should be a numerator of qn+2 limbs and a
normalized divisor of at most qn+1 limbs.

As long as qn > dn - 1, we loop computing one exact quotient limb and
eliminating one high limb of the numerator per iteration, decrementing qn.
(This phase would be equivalent to a call to sbpi1_div_qr, I guess).

Once we reach qn = dn - 1, keep looping to produce quotient limbs, but
also discard one limb of dp in each interation, until we in the final
iteration have qn = 1, qn+2 = 3 numerator limbs, and qn+1 = 2 divisor
limbs, i.e., a single udiv_qr_3by2 (where we might consider skipping the
adjustment steps). I haven't done the error analysis, but I would expect
that errors are similar to the current code.

This is analogous to the discussion of mpn_bdiv_q, which we discussed
years ago, where N and Q are of the same size, and D possibly smaller.
The structure with two loops is also similar to the old mpn_divrem
function with fractional limbs; it would be cleaner if we could have a
function doing only the fractional part, i.e., require that dn == qn +
1.

Fixing this should provide a small speedup for mpn_div_q, since 

  new_nn = 2 * qn + 1;
  ...
  cy = mpn_lshift (new_np, np + nn - new_nn, new_nn, cnt);
  ...
  qh = mpn_sbpi1_divappr_q (tp, new_np, new_nn, new_dp, qn + 1, dinv.inv32);

implies that new_nn and the shift are twice as large as they need to be
(new_nn should be only qn+2).

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.

___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: dead code in div_q.c?

[Niels Möller]

ni...@lysator.liu.se (Niels Möller) writes:

> --- a/tests/mpn/t-div.c   Wed Apr 25 07:38:14 2018 +0200
> +++ b/tests/mpn/t-div.c   Wed Apr 25 22:19:17 2018 +0200
> @@ -445,6 +445,7 @@ main (int argc, char **argv)
> alloc = itch + 1;
>   }
> scratch[itch] = ran;
> +   MPN_COPY (qp, junkp, nn - dn + 1);
> mpn_div_q (qp, np, nn, dup, dn, scratch);
> ASSERT_ALWAYS (ran == scratch[itch]);
> ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == 
> qran1);

Committed.

> *** /tmp/extdiff.0aBNfO/gmp.765c2c27523b/mpn/generic/div_q.c  2018-04-25 
> 21:55:51.457769871 +0200
> --- /home/nisse/hack/gmp/mpn/generic/div_q.c  2018-04-25 21:18:01.516501993 
> +0200
> *** mpn_div_q (mp_ptr qp,
> *** 171,184 
> if (cy == 0)
>   qp[qn - 1] = qh;
> !   else if (UNLIKELY (qh != 0))
> ! {
> !   /* This happens only when the quotient is close to B^n and
> !  mpn_*_divappr_q returned B^n.  */
> !   mp_size_t i, n;
> !   n = new_nn - dn;
> !   for (i = 0; i < n; i++)
> ! qp[i] = GMP_NUMB_MAX;
> !   qh = 0;   /* currently ignored */
> ! }
>   }
> else  /* divisor is already normalised */
> --- 171,176 
> if (cy == 0)
>   qp[qn - 1] = qh;
> !   else
> ! ASSERT_ALWAYS (qh == 0);
>   }
> else  /* divisor is already normalised */

And committed this first part (but with ASSERT instead of ASSERT_ALWAYS);.

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.
___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: documentation on internals not up to date

[Niels Möller]

"Marco Bodrato"  writes:

> ... and it's not just a problem of supporting mpq_denref() we should also
> change some mpq_ functions. If someone does init an mpq_t with:
>   mpz_init (mpq_numref (accumulator));
>   mpz_roinit_n (mpq_denref (accumulator), _one, 1);
>
> then a simple, mpq_add (accumulator, accumulator, another_mpq) might fail.

I don't think we should document support for that code. If the user
initializes parts of an mpq with a read-only mpz objects, the user ought
to treat also the mpq object as readonly.

That said, it sounds reasonable to me to use an an unallocated one for
the internal initialization of the denominator. Then we have to support
that the user accesses it and assigned to it like any other mpz object.
But that doesn't mean that we have to promise that mpz objects the user
creates using ro_init behaves in the same friendly way.

E.g., in case we for some reason change the internals to break
assignment to non-zero unallocated objects, we could hack mpq_denref to
do a real allocation if needed (except that the macro is documented to
work on a const mpq_t). Or go back to using an allocated one.

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.
___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: documentation on internals not up to date

[Marco Bodrato]

Ciao,

Il Ven, 27 Aprile 2018 10:50 am, Niels Möller ha scritto:
> I'd prefer that we not document any way to pass _roinit values to any
> mpz functions taking a non-const mpz_t, even if it happens to work in
> the current implementation. Maybe as a later extension, *if* we find
> some use cases where it provides a significant advantage.

Ok, you are right, I agree.

By the way, "it happens to work" because of the lazy allocation.

> I'd also consider backing that up by having the realloc function abort
> in case _mp_alloc == 0 but _mp_size != 0.

With an ASSERT? Seems unnecessary to me.

Well, we decided what we should not write. Now, let's decide what to write
:-)

Ĝis,
m

-- 
http://bodrato.it/

___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel

Re: divappr interface

[Niels Möller]

ni...@lysator.liu.se (Niels Möller) writes:

> Once we reach qn = dn - 1, keep looping to produce quotient limbs, but
> also discard one limb of dp in each interation, until we in the final
> iteration have qn = 1, qn+2 = 3 numerator limbs, and qn+1 = 2 divisor
> limbs, i.e., a single udiv_qr_3by2 (where we might consider skipping the
> adjustment steps). I haven't done the error analysis, but I would expect
> that errors are similar to the current code.

In the current mpn_sbpi1_divappr_q, there's a curious flag (or mask)
variable used in this loop. It's initially all ones, cleared if we ever
fail to cancel the top limb of the current partial remainder. When the
flag is cleared, remaining quotient limbs are set to GMP_NUMB_MAX.

Torbjörn, Paul, do you remember the analysis behind this? (Code is since
2009...).

I would guess it might happen that even if higher quotient limbs are all
correct, we might get non-zero high limb in

  partial remainder - GMP_NUMB_MAX * truncated D 

If we set the rest of the quotient limbs to GMP_NUMB_MAX, then the
quotient should be large enough thanks to the low end divisor limbs
which we ignored in the truncation.

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.
___
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel