Minimal bindings for optional arb package

[cheuberg]

In #16747, the optional arb package is proposed for inclusion in sage. That ticket did not have any bindings for the Sage library.

In the present ticket, rather minimal bindings are provided.

As a proof of concept, the digamma function of a RealIntervalFieldElement is implemented by converting it to an RealBall, computing the digamma function via the arb library, and converting the result back.

Depends on #16747
Depends on #17688

CC: @sagetrac-alina @fredrik-johansson @sagetrac-ktkohl @sagetrac-skropf @jdemeyer @mezzarobba

Component: numerical

Keywords: arb, digamma function, real interval field

Author: Clemens Heuberger, Marc Mezzarobba

Branch/Commit: 195a6ee

Reviewer: Marc Mezzarobba

Issue created by migration from https://trac.sagemath.org/ticket/17194

[fredrik-johansson]

comment:2

This probably does not matter as long as the bindings are not part of the sage coercion world, but I suppose then you'd want to handle precision in a similar way to RealField, RealIntervalField etc. where it is an attribute of the parent (maybe with a RealBallField type?).

[fredrik-johansson]

comment:3

Sage should arguably have a RealBallField and ComplexBallField with concrete implementations RealBallField_arb and ComplexBallField_arb (it would be useful to provide alternative implementations, for example to provide actual circular complex balls).

Anyway, I took a quick look at the code, and it seems fine to me.

[cheuberg]

comment:4

The question is how the relation to RealIntervalField and ComplexIntervalField should be.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

eb8b584	Trac #17194: Fix typo in one doctest
a6cb2ca	Trac #17194: declare return types
55bdc8b	Trac #17194: mark all doctests in real_arb.pyx as "optional - arb"

[sagetrac-git]

Changed commit from 4a34322 to 55bdc8b

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

14f599f

Trac #17194: rename attribute precision to _precision_

[sagetrac-git]

Changed commit from 55bdc8b to 14f599f

[cheuberg]

comment:8

As a follow-up, minimal bindings for the acb type (complex balls) are provided in #17218.

[fredrik-johansson]

comment:9

The code looks ok to me. A small detail is that precision must be >= 2, not just positive (this is required by RealField et al too).

[sagetrac-git]

Changed commit from 14f599f to 11df4b4

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

11df4b4

Trac #17194: Clarify that precision >=2 is required

[sagetrac-git]

Changed commit from 11df4b4 to 141aca0

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

141aca0

Trac #17194: >2 should have been >= 2 in previous commit

[jdemeyer]

comment:12

Please don't use is_package_installed(). That's extremely slow, you don't want to call it every time.

[jdemeyer]

comment:13

I disagree with

In its present minimal form, this is not part of the sage coercion world.

You should have a Parent from the start, adding it afterwards is only going to make things more difficult.

[jdemeyer]

comment:14

And Arb should probably be renamed to RealBallElement.

[cheuberg]

comment:15

Replying to @jdemeyer:

Please don't use is_package_installed(). That's extremely slow, you don't want to call it every time.

What alternative do you recommend? Calling is_package_installed() once and storing the result in a module global variable? Or just trying to call the arb code and catching the exception?

Concerning renaming Arb to RealBallElement: real_mpfi aka RealIntervalFieldElement and arb have some similarities in their goals; the implementation is different, there might be more functions in the arb library, and there are the corresponding functions for complex balls. Do we want to stress the difference between RealBallElement and RealIntervalFieldElement? Or would RealBallElement be some kind of special case of RealIntervalFieldElement.

I'll try to implement a parent, but it will take some time. The parent would then be RealBallField; what about RealBallFieldElement for the current Arb?

[jdemeyer]

comment:16

Replying to @cheuberg:

What alternative do you recommend?

Try the import and catch ImportError. Would that work?

[fredrik-johansson]

comment:17

You could possibly make it so that one can choose between two different implementations of RealIntervalField and ComplexIntervalField, just like one can choose between:

PolynomialRing(ZZ, implementation='NTL')

PolynomialRing(ZZ, implementation='FLINT')

So you might look up how those are implemented. I don't have an opinion on whether this is better than what you're doing right now.

[jdemeyer]

comment:18

Replying to @cheuberg:

Do we want to stress the difference between RealBallElement and RealIntervalFieldElement?

Yes, I would do that. The implementations are also quite different.

I'll try to implement a parent, but it will take some time. The parent would then be RealBallField; what about RealBallFieldElement for the current Arb?

I slightly prefer the shorter name RealBallElement (like RealDoubleElement), but I can live with RealBallFieldElement.

[sagetrac-git]

Changed commit from 141aca0 to a0afebd

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

dc23341	Trac #17194: Renamed 'Arb' to 'RealBallElement'
a0afebd	Trac #17194: Replace is_package_installed by try ... except ImportError

[cheuberg]

comment:20

Replying to @jdemeyer:

Try the import and catch ImportError. Would that work?

Yes, it works, thanks.

[mezzarobba]

comment:67

Replying to @cheuberg:

Replying to @mezzarobba:

Replying to @cheuberg:

The semantics of != differ from those of real intervals: A ball is considered non-zero if it does not contain 0.

So, != is not not ==.

Hmm, if I read correctly, that's also how comparisons of RealIntervalFieldElements work. Can you please be more specific about the difference you're seeing? For me the only difference is that, in the case of RealIntervalFieldElements, __nonzero__ is not consistent with != 0.

I hope that the following examples make my impression of the current behaviour clear.

sage: a = RIF(-1, 1)
sage: from sage.rings.real_arb import RealBallField
sage: RBF = RealBallField()
sage: b = RBF(a)
sage: a.__nonzero__()
True
sage: b.__nonzero__()
False
sage: bool(a)
True
sage: bool(b)
False
sage: a == 0
False
sage: b == 0
False
sage: a != 0
False
sage: b != 0
False

So I think we agree: == and != behave the same way for intervals and balls (with your implementation), except for the case of x == y when x is y. However, __nonzero__ behaves differently. (And, to be clear: I believe that your choice is the more sensible one, but I don't like the inconsistency.) Is that correct?

So if I understand you correctly, you'd be in favor of having methods RealBall.lt and leave the original methods untouched? I fear that a user may then mistakenly use the above example and get "surprising" results.

No, I guess I wasn't clear, sorry. Imho the fact that comparisons of sage objects (or at least of sage Elements) can fall back to comparing types at all is dreadful. I was saying (as an aside, not as a remark on your code) that I would be tempted to implement safer comparison methods in addition to your implementation of Python comparisons (and with the same semantics in the case where we compare balls to balls) because Sage comparisons are dangerous in general.

Perhaps we could overwrite __richcmp__ in such a way that, if no coercions to balls are available, a NotImplementedError is raised?

I don't know if it is possible, and I doubt it is a good idea to circumvent the coercion system even if it can lead to problems.

I do not understand what you mean. In my understanding, having the implication that identical objects are mathematically equal, should not lead to unexpected behavior within the rest of sage.

Forget it, I think I convinced myself that it was ok—as long as a ball is understood as representing a single real number, of course. (And never as a subset of the real line whose image by a given function needs to be bounded, say.)

For me, the behavior of RealIntervalFieldElement

sage: a = RIF(-1, 1)
sage: from sage.rings.real_arb import RealBallField
sage: RBF = RealBallField()
sage: b = RBF(a)
sage: a == a
False
sage: b == b
True

is very unexpected.

I wouldn't say that: after all, it does not hold true that for all x, y in [-1,1], x=y...

After all, equality is usually assumed to be reflexive.

Yes, but that argument would be (almost) equally valid in the case of RBF(foo) == RBF(foo), and I don't think we want that to hold if RBF(foo) is not a point ball!

(I will be offline for a few days. If anyone would like to take over the review, please feel free.)

[cheuberg]

comment:68

Replying to @mezzarobba:

So I think we agree: == and != behave the same way for intervals and balls (with your implementation), except for the case of x == y when x is y. However, __nonzero__ behaves differently. (And, to be clear: I believe that your choice is the more sensible one, but I don't like the inconsistency.) Is that correct?

This is correct.

[cheuberg]

comment:69

(Replying in several chunks due to trac timeouts when I try to reply in one go, sorry).

Replying to @mezzarobba:

No, I guess I wasn't clear, sorry. Imho the fact that comparisons of sage objects (or at least of sage Elements) can fall back to comparing types at all is dreadful.

I agree.

I was saying (as an aside, not as a remark on your code) that I would be tempted to implement safer comparison methods in addition to your implementation of Python comparisons (and with the same semantics in the case where we compare balls to balls) because Sage comparisons are dangerous in general.

We can always do that later because it would not change behaviour.

Perhaps we could overwrite __richcmp__ in such a way that, if no coercions to balls are available, a NotImplementedError is raised?

I don't know if it is possible, and I doubt it is a good idea to circumvent the coercion system even if it can lead to problems.

So let's forget this idea.

[cheuberg]

comment:70

Replying to @mezzarobba:

I do not understand what you mean. In my understanding, having the implication that identical objects are mathematically equal, should not lead to unexpected behavior within the rest of sage.

Forget it, I think I convinced myself that it was ok—as long as a ball is understood as representing a single real number, of course. (And never as a subset of the real line whose image by a given function needs to be bounded, say.)

I do not understand your point. If M is a subset of the real line (represented by a ball) and f is a function, I still believe that f(M) can be computed by RealBall.

I wouldn't say that: after all, it does not hold true that for all x, y in [-1,1], x=y...

After all, equality is usually assumed to be reflexive.

Yes, but that argument would be (almost) equally valid in the case of RBF(foo) == RBF(foo), and I don't think we want that to hold if RBF(foo) is not a point ball!

No, I do not want that.

For reference, I quote here from the Arb manual:

"Identical pointers are understood to give permission for algebraic simplification. This assumption is made to improve performance. For example, the call arb_mul(z, x, x, prec) sets z to a ball enclosing the set {t^2 : t ∈ x} and not the (generally larger) set {tu : t ∈ x, u ∈ x}."

[mezzarobba]

comment:71

Replying to @cheuberg:

I was saying (as an aside, not as a remark on your code) that I would be tempted to implement safer comparison methods in addition to your implementation of Python comparisons (and with the same semantics in the case where we compare balls to balls) because Sage comparisons are dangerous in general.

We can always do that later because it would not change behaviour.

Yes.

[mezzarobba]

comment:72

So despite what I said about being offline for a few days, I could have a look at your last changes before disappearing.

Replying to @cheuberg:

Forget it, I think I convinced myself that it was ok—as long as a ball is understood as representing a single real number, of course. (And never as a subset of the real line whose image by a given function needs to be bounded, say.)

I do not understand your point. If M is a subset of the real line (represented by a ball) and f is a function, I still believe that f(M) can be computed by RealBall.

Ok, that was very badly expressed once again. Sorry about that, I'm trying to do too many things at once :-(.

All I was trying to say essentially is that if your think of balls as ranges, not approximate reals, then code like

def f(x, y): return x-y
def range(f, *args): return f(*args)
unit_ball = RBF(0, rad=1)
range(f, unit_ball, unit_ball)

looks correct, while it violates the rule about aliasing in arb that you quoted. I guess one could come up with more subtle examples with cached_functions or additional layers of abstraction.

A few more quick comments:

• I think we will need fmpr again later on to implement rad();
• _an_element_ can be inherited, but perhaps this requires some of the coercion & category stuff I implemented on the other ticket, I don't remember;
  (so several of your last changes will probably be reverted in the follow-up ticket)
• we should add a note about aliasing similar to that which you quoted to the python package's documentation too;
• arb_to_mpfi will probably crash sage if the ball's endpoints are not representable within the exponent range of MPFR FP numbers (but it is close to impossible to construct such balls at this point), use sig_on() or sig_str() (without _do_sig) to prevent that;
• I personally find return (lst is rt) or (arb_is_exact(...) and ... and ...) more readable than your implementation of equality, but that's a matter of taste.

None of the above is really urgent or important and the current state of the code looks good to me, so I'm setting the ticket to positive review to avoid delaying it too much. But please feel free to make changes here if you want, and then I'll try to review them when I come back.

[mezzarobba]

Reviewer: Marc Mezzarobba

[cheuberg]

comment:73

Replying to @mezzarobba:

So despite what I said about being offline for a few days, I could have a look at your last changes before disappearing.

Thank you.

• we should add a note about aliasing similar to that which you quoted to the python package's documentation too;
• arb_to_mpfi will probably crash sage if the ball's endpoints are not representable within the exponent range of MPFR FP numbers (but it is close to impossible to construct such balls at this point), use sig_on() or sig_str() (without _do_sig) to prevent that;
• I personally find return (lst is rt) or (arb_is_exact(...) and ... and ...) more readable than your implementation of equality, but that's a matter of taste.

None of the above is really urgent or important and the current state of the code looks good to me, so I'm setting the ticket to positive review to avoid delaying it too much. But please feel free to make changes here if you want, and then I'll try to review them when I come back.

While I think that your points above are valid, I do not want to make changes here any more.
I'll try to do that in a separate ticket. Moreover, I plan to adapt the ComplexBalls (#17218) to the structure implemented here.

[fredrik-johansson]

comment:74

Replying to @mezzarobba:

A few more quick comments:

• I think we will need fmpr again later on to implement rad();

Anything fmpr can be done with arf module.

• arb_to_mpfi will probably crash sage if the ball's endpoints are not representable within the exponent range of MPFR FP numbers (but it is close to impossible to construct such balls at this point), use sig_on() or sig_str() (without _do_sig) to prevent that;

Yes, unfortunately I haven't added any code yet to underflow or overflow gracefully when converting to MPFR. When developing, I deliberately wanted this to crash rather than fail silently.

I should probably fix arb_get_interval_mpfr to set the endpoint to +/- 2^MPFR_MIN_EXP on underflow and +/- inf on overflow.

Rather than changing arf_get_mpfr, I might add a separate arf_get_mpfr_overflow that behaves as if you had computed an out-of-range value within MPFR (setting MPFR overflow flags and so on).

[mezzarobba]

comment:75

Replying to @cheuberg:

While I think that your points above are valid, I do not want to make changes here any more.
I'll try to do that in a separate ticket. Moreover, I plan to adapt the ComplexBalls (#17218) to the structure implemented here.

Great, thanks for your work on these bindings! I'll do what I can to help.

[mezzarobba]

comment:76

Replying to @fredrik-johansson:

Replying to @mezzarobba:

A few more quick comments:

• I think we will need fmpr again later on to implement rad();

Anything fmpr can be done with arf module.

What should we do to convert a mag_t to an mpfr_t? Using mag_get_fmprlooked like the simplest method.

[fredrik-johansson]

comment:77

Replying to @mezzarobba:

Replying to @fredrik-johansson:

Replying to @mezzarobba:

A few more quick comments:

• I think we will need fmpr again later on to implement rad();

Anything fmpr can be done with arf module.

What should we do to convert a mag_t to an mpfr_t? Using mag_get_fmprlooked like the simplest method.

arf_set_mag, so-called because the arf module depends on the mag module and not vice versa.

[mezzarobba]

comment:78

Replying to @fredrik-johansson:

arf_set_mag, so-called because the arf module depends on the mag module and not vice versa.

Ok, I missed that function, thanks!

[mezzarobba]

comment:79

Another question that we may want to think about is whether and how to expose exact (by which I mean ARF_PREC_EXACT) operations. They would be useful to me, but perhaps being able to use them directly from cython is enough.

[fredrik-johansson]

comment:80

ARF_PREC_EXACT is something of a hack, and not really safe in general. If your numbers are too large, you get undefined behavior.

[cheuberg]

comment:81

Replying to @mezzarobba:

• we should add a note about aliasing similar to that which you quoted to the python package's documentation too;

This is now #17809.

• arb_to_mpfi will probably crash sage if the ball's endpoints are not representable within the exponent range of MPFR FP numbers (but it is close to impossible to construct such balls at this point), use sig_on() or sig_str() (without _do_sig) to prevent that;

It was hard to provoke this with available functions, but it did crash sage. This is now #17811.

• I personally find return (lst is rt) or (arb_is_exact(...) and ... and ...) more readable than your implementation of equality, but that's a matter of taste.

This is #17810.

[vbraun]

Changed branch from u/cheuberg/17194-arb to 195a6ee