Wrapping Calcium

[fredrik-johansson]

Here are my current thoughts on a Nemo interface to Calcium.

We'd want the following Julia types, wrapping the corresponding C types:

• ca
• ca_mat
• ca_poly
• ca_vec (optional - we don't have vectors for Flint types, but perhaps we should)

Each such element needs to hold a reference to a parent object. A parent object will wrap a ca_ctx_t defining evaluation settings. A parent object will also hold some additional flags declaring what kind of mathematical structure is being represented. The idea is that you should be able to construct a genuine complex field where infinities will never show up (1 / 0 will throw a division by zero exception instead of returning unsigned infinity). Likewise, you should be able to construct a real algebraic field where you cannot have complex numbers, transcendental numbers, or infinities (actually, such a parent should allow complex field generators internally, e.g. roots of unity, but would only allow constructing elements that lie in real subfields of such fields). We'd accordingly want a few different parent constructors:

• CaComplexField()

  Exact field of complex numbers represented using Calcium elements.

• CaComplexAlgebraicField()

  Exact field of complex algebraic numbers represented using Calcium elements.

• CaComplexFieldExtended()

  Exact field of complex numbers represented using Calcium elements, extended with special values (undefined, unsigned infinity, complex signed infinities). Note that this structure is not a mathematical field, but it represents an exact field as long as calculations are restricted to finite values.

• CaRealField()

  Exact field of real numbers represented using Calcium elements.

• CaRealAlgebraicField()

  Exact field of real algebraic numbers represented using Calcium elements.

• CaRealFieldExtended()

  Exact field of real numbers represented using Calcium elements, extended with special values (undefined, unsigned infinity, real signed infinities). Note that this structure is not a mathematical field, but it represents an exact field as long as calculations are restricted to finite values.

The names can be discussed. Each constructor should accept various optional settings (https://fredrikj.net/calcium/ca.html#context-options).

In addition, we probably want to expose the qqbar type. The corresponding parent object(s) could be singletons like ZZ. Not sure about names (AbsoluteRealAlgebraicField/AbsoluteComplexAlgebraicField, MinpolyRealAlgebraicField/MinpolyComplexAlgebraicField, CanonicalQQBar, or simply QQBar...); we'd want to emphasize that this represents the field of algebraic numbers as a mathematical domain but that users generally will want to choose a Ca field for number-crunching over QQbar.

Symbolic expressions (fexpr) could be useful to wrap internally for glue code purposes, but I don't know if it makes sense to expose them publicly as part of Nemo.

[fredrik-johansson]

Some important implementation issues:

• All ca_t elements must be cleared before the parent ca_ctx_t is cleared. This should not be a problem with exact GC. For some reason I keep having problems in Python where Python deallocates the context object before it deallocates the elements, even though there are live references left, but as far as I can tell this is a Python stupidity (or a bug in my Python code); the C functions should be fine.

• Calcium elements and context objects cannot safely be accessed simultaneously by different threads. This means that all ccalls must be blocking. I'm not sure if there is a way to do automatic fine-grained locking in Julia so that Calcium calculations can be multithreaded if the user creates separate context objects for separate threads. Alternatively, we could allow the user to create unsafe parent objects that perform non-blocking ccalls, with the user taking responsibility to use separate parents in separate threads (e.g. CaComplexField(..., threadblock=false)).

[thofma]

Could you provide a link to the python wrapper?

[fredrik-johansson]

https://fredrikj.net/calcium/pycalcium.html and https://github.com/fredrik-johansson/calcium/blob/master/pycalcium/pyca.py

Note that it is not based around having parents, it's exposing context objects more directly.

[thofma]

OK. I think as a start we should try to build binaries. Does calcium work with the latest flint/arb/antic releases?

[fredrik-johansson]

Yes, it should.

[thofma]

Binaries incoming: https://github.com/JuliaPackaging/Yggdrasil/pull/2992

Regarding the finalizers. We have to be a bit careful here, as julia does not guarentee the order in which finalizers are run (as @tthsqe12 reminded me of recently). So we have to do the reference counting ourselves, similar how it is implemented in Singular.jl or https://github.com/JuliaLang/julia/pull/20124.

[fredrik-johansson]

So I take it you need to put a refcount field on the parent, increment when you initialize an element, and decrement when you destroy an element. Then, how do you ensure that the context object gets cleared when (refcount == 0 AND the Julia parent object is no longer referenced in Julia)? You can't just clear the context object when refcount == 0 since the parent object itself may still be live.

[thofma]

I probably misunderstood. I thought the parent object is the context object?

[fredrik-johansson]

What I mean is that the Julia object wrapping the C context object may be active even if there are no active C-level references to the wrapped context object. Let's say the user creates K = CaField(). Then creates x = K(1). Then destroys x. At this point K.refcount == 0 and there are no pointers to the context object inside K, but you don't want to clear the context object because the user may yet want to create y = K(2).

Edit: corrections

[thofma]

Since the object created by CaField() is still referenced by K and K did not go out of scope, the object is not marked to be finalized.

[fredrik-johansson]

Sure. The question then is then how you ensure that, once both x and K are marked to be finalized, everything occurs in order. I'm maybe missing something obvious here (like K contributing 1 to its own refcount?).

[thofma]

Sorry, I misunderstood, I see the problem. I will have to check what we did before.

[fredrik-johansson]

Isn't this an issue with some Flint types already, though? For example, fmpz_mpoly_clear needs a context object.

[wbhart]

Some of those functions that need a context object actually don't use it, so we just don't pass it, which means the context object can be cleaned up first. In Singular.jl we actually do ref counting manually because the Singular objects actually can't be cleaned up without the ctx.

This is all fine on x86_64/x86 linux at least, since the first few parameters to functions are passed in registers, so if the parameter is not used it doesn't matter if it is never supplied. The reg will just have a random value.

I am not certain, but I thought that on Windows the parameters were also passed on the stack, so I'm not sure what happens there. I don't know about other arches.

[tthsqe12]

The flint code is not a good reference on this issue because it does not ensure that parent is alive or valid when the elements are cleared. https://github.com/wbhart/flint2/issues/773

The way it is done in singular.jl is probably the way to go here. https://github.com/oscar-system/Singular.jl/blob/master/src/number/NumberTypes.jl

function _Ring_finalizer(cf)
   cf.refcount -= 1
   if cf.refcount == 0
      libSingular.nKillChar(cf.ptr)
   end
   nothing
end

function _Elem_finalizer(n)
   R = parent(n)
   libSingular.n_Delete(n.ptr, R.ptr)
   _Ring_finalizer(R)
   nothing
end

[wbhart]

That's not necessary if the parent object is also in each element. We didn't do that in Singular.jl because pointers are used rather than parent objects containing those pointers.

It's definitely worth avoiding custom reference counting where possible.

[thofma]

I think we have to do the counting here. Last time I checked with @tthsqe12, the order of the finalizers was arbitary.

[wbhart]

No, I think you are right after all. The problem is if all elements AND the ring are scheduled for finalisation, then it can run the finalisers in arbitrary order, cleaning up the ring before the elements that need to be cleaned up.

So in fact we do really rely on Flint never actually needing the finalizers to do cleanup.

[wbhart]

Yeah I remember now. Singular stores tables of function pointers, one of which is the actual cleanup function. So the ring really is needed for the cleanup. We don't do that in Flint as there are no generics in Flint.

I assume Calcium actually does use its context objects for cleanup in some non-trivial way (i.e. other than just passing them for interface uniformity)?

[wbhart]

Actually, I wonder if the problem really is finalizers being called in arbitrary order, or gc being called during finalization. It's feasible that by keeping the ring object alive during the finalization of an element that it will not get cleaned up in such a gc cycle.

This wasn't possible in Singular.jl due to the fact that only the pointer was passed to Singular for the cleanup, not the julia object. That doesn't happen for Flint.

We should look into all this to see what's actually going on.

[wbhart]

Yeah, I think this might be it. We should try adding @gc_preserve to the ring object in all of Singular's cleanup functions and ditch the reference counting. I wouldn't mind betting that actually works.

[fredrik-johansson]

I think this matters also at least for Antic nf_elem, where nf_elem_clear needs to read the storage type from the nf_t.

[wbhart]

Yep, and we don't get crashes with nf_elem, so I think I was right about it not being a problem so long as all pointers passed to C cleanup functions are Julia objects. When Singular.jl works again I'll check my hypothesis.

[wbhart]

Actually, same applies for fq_default. No problems there either.

[fredrik-johansson]

@thofma How did it go with the binaries? If the dependency is sorted out and there is some working stub code in place to access libcalcium from Nemo, maybe I can try to wrap a couple of functions.

[thofma]

Binaries are available. You can use them like the flint ones, see src/Nemo.jl. There are some minor things to do for julia version < 1.3, but I can quickly do this once you have a PR with a couple of functions.

[fredrik-johansson]

I also need to to add it to Project.toml, right? And elsewhere? How do I retriev the uuid?

[fredrik-johansson]

Nevermind, I think I got it.

[thofma]

1067

[fredrik-johansson]

I will do a basic qqbar wrapper.

[fredrik-johansson]

Super basic initial code here: https://github.com/fredrik-johansson/Nemo.jl/tree/calcium

I will do some more work before I create a PR.

[fredrik-johansson]

Regarding random generation: I have qqbar_randtest, qqbar_randtest_real, qqbar_randtest_nonreal in C, but they need Flint random numbers. I guess I need to reimplement them in Julia?

[thofma]

We are already using flint_rand_t or whatever the name is. Look for rand_ctx() and _flint_rand_states (the latter is the global variable holding the thread local contexts).

[fredrik-johansson]

@thofma Can you upgrade Calcium_jll to 0.4.0, which I just tagged?

It can wait until after #1070 is merged, but I will need it to wrap ca_t.

[thofma]

Done. It is 0.400.0.

[fredrik-johansson]

Thanks! How do I update the dependency in Nemo and force Julia to install the latest version? Just setting Calcium_jll = "~0.400.0" in Project.toml doesn't do anything.

[wbhart]

I think you also have to do ]up

[fredrik-johansson]

That did it. Thanks!

[fredrik-johansson]

A small technical problem for wrapping polynomials: the top coefficient of a ca_poly can be non-provably zero, in which case the exact length/degree isn't known. I see three possible solutions:

1. Have length() and degree() and similar methods return upper bounds, like arb_poly/acb_poly. There are obvious problems with this (namely that it's mathematically wrong).

2. Allow constructing such polynomials, but have degree() and length() and similar methods compare the top coefficient against zero and possibly throw an exception. The drawback is that each access to the length of a polynomial will involve a computation, which can be slow (unless we cache some extra flags with the polynomial, which I'd rather avoid).

3. Make it an invariant of the Julia-wrapped ca_poly that the top coefficient is nonzero: verify that the top coefficient is nonzero at the time of construction/modification of a ca_poly and throw an exception otherwise. Then degree() and length() and similar methods can just access the integer values.

I think 3. is the way to go. There will be some overhead for certain operations where an extra zero test needs to be done, but it feels like a more robust solution.

(Relative) power series should presumably have the same behavior but reversed; there, we want to check the constant coefficient, but there is no need to check the top coefficient.

[thofma]

I agree that 3. sounds like the best option.

[wbhart]

Sounds like a reasonable solution.

I'm not sure if it is possible to always figure out whether a coefficient for relative power series is zero from the beginning, e.g. if the coefficient ring is not exact. But maybe the idea helps in some cases.