Problem with numerical evaluation of QQBar.

[PeterLuschny]

Consider the program below. It runs until n = 44 and then silently gives up for no apparent reason.

Tested with version Nemo: 0.47.3, Julia 1.11.1
on Windows 11/23H2, and Mac Sonoma 14.4.1/M2.

using Nemo

RR = ArbField(1000)

function F(m)
    sum(RR(tanpi(QQBar(n) / (1 + 2 * m))^(2 * m)) for n in 0:m)
end

function aList(upto)
    for m in 0:upto
        i = unique_integer(F(m))
        if i[1] == true
            println(m, " ", i[2])
        else
            println(m, "Error")
        end
    end
end
aList(50)

[fingolfin]

What do you mean by "silently gives up" ?

For me it indeed "slows down" when computing step 45, but it successfully completes the loop (on a MacBook with M1 Max, macOS 12, with 64 GB of RAM).

[PeterLuschny]

OK, it did not complete the loop on my systems in any reasonable time, but the systems also did not have 64 GB of Ram. What remains is the observation that even on better equipped systems a serious change of state occurs at this point, the cause of which is not clear.

…

Message ID: ***@***.***>

[thofma]

It seems that for many values of n, m the following is fine, but then for this specific values it is hanging/taking a very long time:

julia> n, m = 43, 45;

julia> s = tanpi(QQBar(n) / (1 + 2 * m))^(2 * m);

julia> ArbField(10)(s) # hanging here. Replace 10 by anything else.

This in the end calls qqbar_get_arb. I think there are some heuristics when computing enclosures, maybe they are not optimal here. Maybe @fredrik-johansson knows.

@PeterLuschny You probobably know this, but it seems to be more efficient to do all computations with approximations? If I do

 function F(m)
    sum((tanpi(RR(n) / (1 + 2 * m))^(2 * m)) for n in 0:m)
end

then everything is pretty fast (you might have to increase the precision for large m).

[fingolfin]

Great example, @thofma . So this is fast:

julia> n, m = 42, 45;

julia> s = tanpi(QQBar(n) / (1 + 2 * m))^(2 * m);

julia> ArbField(10)(s)
[2.59e+82 +/- 5.16e+79]

julia> @time ArbField(10)(s)
  0.000020 seconds (1 allocation: 64 bytes)
[2.59e+82 +/- 5.16e+79]

But this is not (just changed n):

julia> n, m = 43, 45;

julia> s = tanpi(QQBar(n) / (1 + 2 * m))^(2 * m);

julia> @time ArbField(10)(s)
494.590022 seconds (20.98 M allocations: 202.481 GiB, 0.69% gc time, 0.00% compilation time)
[4.55e+95 +/- 7.75e+92]

I also took a backtrace to get a rough idea what it is doing:

  * frame #0: 0x000000029ced59f0 libflint.19.0.dylib`sd_fft_basecase_8_1 + 992
    frame #1: 0x000000029ced54bc libflint.19.0.dylib`sd_fft_trunc + 684
    frame #2: 0x000000029cec4728 libflint.19.0.dylib`fft_worker_func + 392
    frame #3: 0x000000029cec5224 libflint.19.0.dylib`mpn_ctx_mpn_mul + 1732
    frame #4: 0x000000029cf77bf4 libflint.19.0.dylib`arf_complex_mul + 3012
    frame #5: 0x000000029d05ed84 libflint.19.0.dylib`_acb_poly_evaluate_mid + 196
    frame #6: 0x000000029d05f048 libflint.19.0.dylib`_acb_poly_refine_roots_durand_kerner + 568
    frame #7: 0x000000029d054e08 libflint.19.0.dylib`_acb_poly_find_roots + 824
    frame #8: 0x000000029d09e374 libflint.19.0.dylib`arb_fmpz_poly_complex_roots + 404
    frame #9: 0x000000029d1fde38 libflint.19.0.dylib`_qqbar_enclosure_raw + 936
    frame #10: 0x000000029d1ffd38 libflint.19.0.dylib`qqbar_get_acb + 440
    frame #11: 0x000000029d1ffffc libflint.19.0.dylib`qqbar_get_arb + 140
...

In the meantime I wonder if @fredrik-johansson or @albinahlback have any insights on this. Perhaps we should just open a FLINT issue for this?

[albinahlback]

I think Fredrik is the person for this, I don't have a lot of insight in Calcium.

[fredrik-johansson]

Some numerical heuristic isn't working properly. Without debugging this carefully, I suspect it's an instance where polynomial root-finding performs poorly due to clustered roots (we don't have cluster detection, so convergence can be extremely slow in bad cases).