It could be useful to write some naive Python code before attempting a C implementation.

A fun test problem would be to compute Gröbner bases: https://ask.sagemath.org/question/46317/groebner-basis-with-transcendental-coefficients/

Big initial todo list.

Urgent (before 0.1 release):

• use field cache with hash table
• fix hash functions
• mctx persistent pointers (add indirection)
• improve qqbar_root_ui
• cache acb enclosures
• implement all ca_check_X methods
• complex parts
• string conversion, better printing
• some example programs

Important:

• continuous integration
• check build system
• stronger tests (check that some basic simplifications really work, not just giving UNKNOWN)
• clean up all structs and macros
• fix/unity field_init methods
• change ideal to ** to avoid copies
• recycle storage in make_field_element
• move nf_t to field object
• optimize field/extension insertion/lookup (avoid duplicate allocations, copies)
• rename ca_check_X methods?
• separate fast and slow check methods
• organize the code for inserting extra field elements -- what is the right place to do this?
• some more simplifications for major transcendental functions (exp, log, some trig functions)
• use factorisation in abs, sign, ...
• convert mixed algebraic number to single number field
• decent field pruning (remove unused elements, all powers of x^(1/n) divisible by n -> x (etc...)
• transfer between context objects
• functions to decompose ca_t -> mpoly_q + ext list (etc.)
• special values for pow
• separate in-field / inert powering
• optimize subtraction

Long-term:

• reference counting for extension numbers and fields (with recycle bins)
• more intelligent adaptive numerical evaluation
• more intelligent qqbar simplification
• algorithms for mixed algebraic/transcendental simplification
• optimized fmpz_mpoly_q composition and evaluation (fmpz_mpoly_q, arb, acb)
• simplifying number fields
• number field embeddings
• lots of optimizations for roots of unity and nth roots
• more transcendentals
• polynomials, matrices
• radical denesting (outside of qqbar, also for radicals involving transcendentals)
• add a simple expression parser (useful for test code)
• Gröbner bases
• good Gröbner basis code / plug in fast external GB engine
• user-defined extension numbers

I was considering updating flint in Debian to 2.9.0, but it looks like it breaks calcium's package. Running ./configure then make is ok, but trying make check gives:

python3 pycalcium/pyca.py
Exception (nf_init). Degree must be at least 1.

and compiling the doc with sphinx also fails because of an abort trying to read pycalcium.

When I leave the rpath settings in the configure script, the Debian checkers complain about the package, so I'm patching them out. I'll attach a patch.

Version 0.4.1
FreeBSD 13.1

qqbar_hash (done)

qqbar_get_fmpq (done)
qqbar_get_fmpz (done)

qqbar_pow_fmpq (done)
qqbar_pow_fmpz (done)
qqbar_pow_si (done)

qqbar_numerator (done)
qqbar_denominator (done)

qqbar_set_fmpz_poly_root_nearest
qqbar_set_fmpq_poly_root_nearest

qqbar_set_fmpz_poly_root_nearest_acb
qqbar_set_fmpq_poly_root_nearest_acb

(or similar)

With the master of flint, arb, antic and calcium, a computation in qqbar fails:

>>> import pyca
>>> two = pyca.qqbar(2)
>>> two_sqrt = two.sqrt()
>>> two_sqrt
1.41421 (deg 2)
>>> two_sqrt + 1
2.41421 (deg 2)
>>> two_sqrt * two_sqrt
2.00000 (deg 1)
>>> (two_sqrt + 1) * two_sqrt
qqbar_enclosure_raw: root refinement not converging
IOT instruction (core dumped)  python3

I was pondering packaging calcium for Debian (experimental at first, since it isn't considered stable yet). I noticed the configure script doesn't go with the rest of the files : it comes from flint & arb, but there it is GPL-2+ -- so I think you should make it more precise in the code.

The problem seems to be that ctypes misinterprets the pointers returned by flint_malloc, ca_mat_entry_ptr, ca_vec_entry_ptr, and ca_poly_coeff_ptr

On the buchberger2 branch (commit a7caebe) is an attempt to implement the improved Buchberger algorithm GRÖBNERNEW2, p. 232 in Becker and Weispfenning. For some reason I can't figure out, it's not working as intended; random input quickly runs into cases where it slows down drastically while the ordinary Buchberger performs fine. Something like build/examples/dft -input 3 16 breaks.

At first I thought this was just due to randomness in the pair selection, but GRÖBNERNEW2 looks systematically worse, which shouldn't be the case. SymPy has an implementation of GRÖBNERNEW2 (https://github.com/sympy/sympy/blob/master/sympy/polys/groebnertools.py) and seems to do fine on a handful of examples I've tried manually.

There is a ca_hash_repr which hashes the internal representation of a ca_t. More interesting is to have a genuine hash function C -> Z to allow using hash tables to represent finite sets of numbers. Such a function must either be trivial (e.g. always return 0) or be allowed to fail (i.e. some numbers are unhashable).

Here is a possible algorithm:

1. Let z = a*x+b where a, b are some fixed complex numbers. Compute a precise enclosure of z.
2. Round the real and imaginary parts of z to 64-bit floating-point numbers with correct rounding, and return a hash of the floating-point representations.

Note:

• This will fail when the correct rounding of z cannot be deduced. This will occur with negligible probability if a, b are random.
• All sufficiently small numbers will hash to the same value, unless one sets b = 0.
• It will not work for numbers that are too large for numerical evaluation with Arb, e.g. exp(exp(exp(exp(exp(pi))))). This could be solved by fixing some huge bound H and returning 0 when |x| > H. Some some form of level-index arithmetic needs to be implemented to deal with such inequalities. Hashing will fail when the inequality cannot be evaluated.
• It will similarly not work for numbers that are too small if b = 0. As above, a lower bound 1 / H could be used.
• This will not be compatible with any conventional hash function for integers or rational numbers, e.g. hash(int) in Python. You could check if the number is an integer, but then hashing will fail for numbers where this can't be proved. Hashing in a way that is compatible with integers is perhaps better left to the user.
• Alternatively, to stay compatible with hashing for integers, you could take a = 1 and b = some number << 1, round z = a*x + b to the nearest Gaussian integer, and return a hash of the representation. You would need some cutoff for huge numbers where computing the nearest integers would be too expensive.

Python tries to free the context object before all the number objects that reference the context object have been destroyed.

Is this some ctypes-related reference counting bug in Python or am I missing something really obvious?

Simply running ./configure (successfully) and then running make resulted in:

    CC   ../build/fmpz_mpoly_q/add.lo
In file included from add.c:12:
/home/albin/cpkg/calcium/fmpz_mpoly_q.h:269:1: error: redefinition of ‘fmpz_mpoly_divexact’
  269 | fmpz_mpoly_divexact(fmpz_mpoly_t res, const fmpz_mpoly_t x, const fmpz_mpoly_t y, const fmpz_mpoly_ctx_t ctx)
      | ^~~~~~~~~~~~~~~~~~~
In file included from /home/albin/cpkg/calcium/fmpz_mpoly_q.h:25,
                 from add.c:12:
/usr/local/include/flint/fmpz_mpoly.h:877:6: note: previous definition of ‘fmpz_mpoly_divexact’ with type ‘void(fmpz_mpoly_struct *, const fmpz_mpoly_struct *, const fmpz_mpoly_struct *, const fmpz_mpoly_ctx_struct *)’
  877 | void fmpz_mpoly_divexact(fmpz_mpoly_t Q, const fmpz_mpoly_t A,
      |      ^~~~~~~~~~~~~~~~~~~
make[2]: *** [../Makefile.subdirs:60: ../build/fmpz_mpoly_q/add.lo] Error 1
make[2]: Leaving directory '/home/albin/cpkg/calcium/fmpz_mpoly_q'
make[1]: *** [Makefile:122: libcalcium.so.0.4.0] Error 2
make[1]: Leaving directory '/home/albin/cpkg/calcium'
make: *** [Makefile:142: library] Error 2

Am I doing something wrong here?

The tests target only builds them.

qqbar binary operations guess and check for a rational solution when the inputs have high degree. This could be generalized to higher degree guesses (2? 4?) when the inputs have high enough degree, as long as things don't slow down noticeably in the generic case.

>>> exp((log(2*i) - pi*i/2)/2)
1.41421 + 0e-34*I {-a*b+a where a = 0.707107 + 0.707107*I [a^4+1=0], b = I [b^2+1=0]}

It would look better if this just output sqrt(2).

Currently a ca_t allows representing any complex number together with infinities. It should be possible to initialize context objects that restrict computation to a smaller domain:

• Algebraic numbers
• Real algebraic numbers
• Real algebraic numbers + infinities
• Real numbers
• Real numbers + infinities
• Complex numbers
• ...

Anything outside the supported set should give Undefined. For example, whereas sqrt(-1) -> i by default, one should have sqrt(-1) -> Undefined when working strictly with real numbers.

Tried to launch the binder notebook, got the following error in the log:

    CC   ../build/qqbar/neg.lo
    CC   ../build/qqbar/evaluate_fmpq_poly.lo
    CC   ../build/qqbar/validate_enclosure.lo
    CC   ../build/qqbar/init.lo
In file included from evaluate_fmpq_poly.c:12:0:
/srv/conda/envs/notebook/include/antic/nf.h:40:1: error: C++ style comments are not allowed in ISO C90
 // Visibility Macros to separate API from non-API and inline template code.
 ^
/srv/conda/envs/notebook/include/antic/nf.h:40:1: error: (this will be reported only once per input file)
../Makefile.subdirs:60: recipe for target '../build/qqbar/evaluate_fmpq_poly.lo' failed
make[2]: *** [../build/qqbar/evaluate_fmpq_poly.lo] Error 1
make[2]: *** Waiting for unfinished jobs....
make[2]: Leaving directory '/home/jovyan/qqbar'
Makefile:121: recipe for target 'libcalcium.so.0.4.1' failed
make[1]: Leaving directory '/home/jovyan'
make[1]: *** [libcalcium.so.0.4.1] Error 2
make: *** [library] Error 2
Makefile:142: recipe for target 'library' failed

but looking at https://github.com/wbhart/antic/blame/trunk/nf.h shows no such comment in the current git version of antic. It looks like the binder link uses conda-forge as a source of antic, so I suspect one could either patch the build script or ask for a new release of antic?

Implement computation of the Jordan canonical form of a matrix.

Composing matrix exponentials and logarithms is a good way to find expressions that don't simplify:

>>> ca_mat([[1,-1],[-1,-1]]).exp().log()
ca_mat of size 2 x 2
[1.00000 {(-b*d*f+b*e*f)/(2*c) where a = 1.41421 [Log(4.11325 {(c+d+e)/2})], b = -1.41421 [Log(0.243117 {(-c+d+e)/2})], c = 3.87013 [Sqrt(14.9779 {d^2+e^2-2})], d = 4.11325 [Exp(1.41421 {f})], e = 0.243117 [Exp(-1.41421 {-f})], f = 1.41421 [f^2-2=0]}, -1.00000 {(b*d^2+b*e^2-2*b)/(c*d*f-c*e*f) where a = 1.41421 [Log(4.11325 {(c+d+e)/2})], b = -1.41421 [Log(0.243117 {(-c+d+e)/2})], c = 3.87013 [Sqrt(14.9779 {d^2+e^2-2})], d = 4.11325 [Exp(1.41421 {f})], e = 0.243117 [Exp(-1.41421 {-f})], f = 1.41421 [f^2-2=0]}]
[-1.00000 {(b*d*f-b*e*f)/(2*c) where a = 1.41421 [Log(4.11325 {(c+d+e)/2})], b = -1.41421 [Log(0.243117 {(-c+d+e)/2})], c = 3.87013 [Sqrt(14.9779 {d^2+e^2-2})], d = 4.11325 [Exp(1.41421 {f})], e = 0.243117 [Exp(-1.41421 {-f})], f = 1.41421 [f^2-2=0]},             -1.00000 {(b*d*f-b*e*f)/(2*c) where a = 1.41421 [Log(4.11325 {(c+d+e)/2})], b = -1.41421 [Log(0.243117 {(-c+d+e)/2})], c = 3.87013 [Sqrt(14.9779 {d^2+e^2-2})], d = 4.11325 [Exp(1.41421 {f})], e = 0.243117 [Exp(-1.41421 {-f})], f = 1.41421 [f^2-2=0]}]

This is the problem:

>>> a = exp(2*sqrt(2))
>>> b = exp(-2*sqrt(2))
>>> sqrt(a+b-2); sqrt(a+1/a-2)
3.87013 {a where a = 3.87013 [Sqrt(14.9779 {b+c-2})], b = 16.9188 [Exp(2.82843 {2*d})], c = 0.0591057 [Exp(-2.82843 {-2*d})], d = 1.41421 [d^2-2=0]}
3.87013 {(a*b-a)/(b) where a = 4.11325 [Sqrt(16.9188 {b})], b = 16.9188 [Exp(2.82843 {2*c})], c = 1.41421 [c^2-2=0]}
>>> 
>>> sqrt(a+b-2) - (a-1)/sqrt(a)
0e-1125 {(a*b-c+1)/(b) where a = 3.87013 [Sqrt(14.9779 {c+d-2})], b = 4.11325 [Sqrt(16.9188 {c})], c = 16.9188 [Exp(2.82843 {2*e})], d = 0.0591057 [Exp(-2.82843 {-2*e})], e = 1.41421 [e^2-2=0]}
>>> 
>>> sqrt(a+1/a-2) - (a-1)/sqrt(a)
0

a+b-2 does not factor, though the equivalent rational function a+1/a-2 does.

The factoring algorithm for field elements somehow needs to be made aware of this. A possibility is to look specifically for exponentials to eliminate, but perhaps there's a more general approach.

We now have a function with this name which conflicts with your version in fmpz_mpoly_q.h.

This is super slow:

>>> from pyca import *
>>> a = 417/(962*pi+80808)
>>> a**50 - a**51 * a**-1

Of course, it works well when increasing POW_LIMIT so that the powers get expanded algebraically:

>>> ctx = ca_ctx(pow_limit=100)
>>> a = ca(pi, context=ctx)
>>> a = 417 / (962 * ca(pi, context=ctx) + 80808)
>>> a**50 - a**51 * a**-1

But this shouldn't be necessary to get acceptable performance.

v = cos(acos(sqrt(2) - sqrt(3))/3)
assert 1 - 90*v**2 + 321*v**4 - 592*v**6 + 864*v**8 - 768*v**10 + 256*v**12 == 0

This equality test interestingly enough works; however, v is not simplified to an algebraic number automatically. This ought to be done.

Evaluating an fmpz_poly or fmpq_poly at a qqbar argument just uses Horner's rule at the moment, requiring repeated expensive qqbar arithmetic operations. It should be possible to do much better, i.e. finding an annihilating polynomial in one go using linear algebra (at least so I've been told -- I don't know the algorithm myself).

The same applies to powering.

ca_set_qqbar should try harder to standardize:

• Maybe sqrt(-3) -> exp(2 pi / 3) and sqrt(5) -> phi, etc. (like Pari does)
• For roots of unity, maybe standardize to exp(2 pi i /n)
• Make a choice when there is real/imaginary symmetry
• Always normalize to a monic polynomial
• After normalizing to a monic polynomial, use some more heuristics; at least try to find a linear transformation that minimizes content
• For low-degree polynomials, implement polredabs / polredbest
• Recognize cyclotomic and trigonometric fields at least in some cases (expensive at high degrees)

1. Flint is getting a fmpz_mpoly_scalar_fmma method that can be useful

2. Some uses of gcd + divexact can be replaced by gcd_cofactors

I'm pretty sure that this is not worth the trouble that it would take to fix, but OSX does not understand time -f "%e" thereby causing make check to fail.

There are two possible ways for a user to fix this.

1. Change the line in makefile.subdirs from @time -f "%e" $< to @time $<. This gives slightly messier output.
2. Download gnu-time from homebrew and change the line to @gtime -f "%e" $<.

Hopefully this note can save someone a bit of time.

Hey there!

Love this library. I'd been wondering what the status of this project was until I jumped onto the Flint mailing lists and noticed from the 3.0 release that Calcium is now incorporated into Flint. I think it'd be useful to document this on the website and repo, so users will be redirected to flint instead.

Thanks!

The bugfix 9f317f4 caused huge_expr -ca to slow down 2x. This may be because a relation is being proved twice.

Add a module for Puiseux series with ca coefficients. Could also be something more general (Hahn series? transseries) but I suppose Puiseux series would be a good start.

Hi Fredrik,

Thanks for developing calcium! Even in 0.x status it is the better choice for some developments, so I started using it.
I did find an issue though... two algebraic numbers comparing equal using qqbar_t yield T_UNKNOWN when comparing them using ca_t.

The test t-ca_equality_20221220.c reproduces the issue.
Calcium version used: commit 59a61324a9a113e269d646691a59273b7e784d04 (current latest on master)
When you find some time, could you have a look?

Thanks and regards,

Michiel

Michiel De Wilde at Keysight Technologies (using my personal GitHub account)

The function qqbar_set_fmpz_poly_root_nearest should have a fast path that tries to compute a nearby root instead of computing all roots.

It also needs to be made public.

There could also be a version that takes an acb_t input.

https://github.com/fredrik-johansson/calcium/blob/c69f42620ae9acf59737e8bf8efaba123f3de5c8/qqbar/test/t-sec_pi.c#L57

• Basic traversal functions should be implemented iteratively rather than recursively so that deeply nested expressions don't risk causing stack overflow
• Transformations such as replace could be implemented in a streaming style to avoid recursive copies
• Efficient in-place aliased operations
• Avoid temporaries in big integer operations
• Other todos in the code