Question on the real_root.c example

[reneruhr]

Hi,
I am trying to extend the real_root.c example to finding the root of f(x)=sinh(2Jx)sinh(2Kx)-1 for J=K=1.
To create this function, I do

    _arb_poly_sinh_series(out, x, xlen, order, prec);
    _arb_poly_sinh_series(factor2, x2, xlen, order, prec);
    _arb_poly_mullow(out, out, order, factor2, order, order, prec);
    _arb_poly_sub(out, out, order, one->coeffs, one->length, prec);

where
arb_poly_one(one),
x2 is setup exactly like x,
and arb_struct factor2[2].
It results in a segfault after a couple of bisection steps.

The segfault is caused by _arb_poly_sinh_series(factor2, x2, xlen, order, prec) from which I deduce that I have not setup factor2 (since it persists when setting x = x2). Taking arb_t factor2 also gives a segfault.

Another issue:
Simplifying the function to f(x)=sinh(2*x)-1, given by

    _arb_poly_sinh_series(out, x, xlen, order, prec);
    _arb_poly_sub(out, out, order, one->coeffs, one->length, prec);

will find the correct root for the interval bounds a=0, b=5:

after bisection 2: [0.439453125000000, 0.444335937500000]
refined root (0/1):
[0.4418938 +/- 3.77e-10]

but will fail for a=0, b=10 which is surprising to me. It yells

after bisection 2: [0.439453125000000, 0.449218750000000]
warning: some newton steps failed! refined root (0/1):
0

I have further removed the if(arb_calc_verbose) clauses appearing in real_root.c since they cause the following compile error, when compiling with c++ with the gcc 11 on an ARM Mac.

Undefined symbols for architecture arm64:
  "___emutls_v.arb_calc_verbose", referenced from:
      __ZN7Onsager3runEv in cc2Pzdb0.o
ld: symbol(s) not found for architecture arm64

Here is the modified example, I have commented code that I added (compared to real_root.c) with ^'s. Any help is appreciated, I hope this is an appropiate place to request for such.

#include "arb.h"
#include "arb_calc.h"

  slong eval_count = 0;

  struct Interactions
  {
    arb_t J;
    arb_t K;
  };

  Interactions interactions;
  
  int critical_temp(arb_ptr out, const arb_t inp, void * params, slong order, slong prec)
  {
    arb_ptr x;
    int xlen = FLINT_MIN(2, order);
    x = _arb_vec_init(xlen);
    arb_set(x, inp);
    if (xlen > 1)
      arb_one(x + 1); 

    arb_ptr x2;                                         //^^^^^^
    x2 = _arb_vec_init(xlen);                                         //^^^^^^
    arb_set(x2, inp);                                         //^^^^^^
    if (xlen > 1)                                         //^^^^^^
      arb_one(x2 + 1);                                         //^^^^^^

    auto inter_ptr = static_cast<Interactions*>(params);                                         //^^^^^^
    auto J = inter_ptr->J;                                         //^^^^^^
    auto K = inter_ptr->K;                                         //^^^^^^
    arb_mul(x, x, J, prec);                                         //^^^^^^
    arb_mul(x2, x2, K, prec);                                         //^^^^^^

    arb_t two;                                         //^^^^^^
    arb_init(two);                                         //^^^^^^
    arb_set_ui(two,2);                                         //^^^^^^
    arb_mul(x, x, two, prec);                                         //^^^^^^
    arb_mul(x2, x, two, prec);                                         //^^^^^^
    arb_clear(two);                                         //^^^^^^

    arb_poly_t one;                                         //^^^^^^
    arb_poly_init(one);                                         //^^^^^^
    arb_poly_one(one);                                         //^^^^^^

    arb_struct factor2[2];  // as in arb_calc::check_block                                          //^^^^^^
    arb_init(factor2);                                         //^^^^^^
    arb_init(factor2+1);                                         //^^^^^^
    //or  arb_t factor2; // as in arb_calc::arb_calc_isolate_roots


    _arb_poly_sinh_series(out, x, xlen, order, prec);                                         //^^^^^^
    _arb_poly_sinh_series(factor2, x2, xlen, order, prec);                                         //^^^^^^
     _arb_poly_mullow(out, out, order, factor2, order, order, prec);                                         //^^^^^^
    _arb_poly_sub(out, out, order, one->coeffs, one->length, prec);                                         //^^^^^^


    _arb_vec_clear(x, xlen);
    _arb_vec_clear(x2, xlen);                                         //^^^^^^
    arb_poly_clear(one);                                         //^^^^^^
    arb_clear(factor2);                                         //^^^^^^
    arb_clear(factor2+1);                                         //^^^^^^
    eval_count++;
    return 0;
  }


  void run()
  {
    arf_interval_ptr blocks;
    arb_calc_func_t function;
    int * info;
    void * params;
    int param1;
    slong digits, low_prec, high_prec, i, num, found_roots, found_unknown;
    slong maxdepth, maxeval, maxfound;
    int refine;
    double a, b;
    arf_t C;
    arf_interval_t t, interval;
    arb_t v, w, z;

    auto inter = Interactions();                                         //^^^^^^
    arb_init(inter.J);                                         //^^^^^^
    arb_init(inter.K);                                         //^^^^^^
    arb_one(inter.J);                                         //^^^^^^
    arb_one(inter.K);                                         //^^^^^^

    params = (void*)(&inter);                                         //^^^^^^
    
    function = critical_temp;

    a = 0;
    b = 5;  // fails for b=10: after bisection 2: [0.439453125000000, 0.449218750000000]   warning: some newton steps failed! refined root (0/1): 0

    refine = 1;
    digits = 5;
    maxdepth = 20;
    maxeval = 1000;
    maxfound = 1;
    low_prec = 60;

    high_prec = digits * 3.32192809488736 + 10;
    found_roots = 0;
    found_unknown = 0;

    arf_init(C);
    arf_interval_init(t);
    arf_interval_init(interval);
    arb_init(v);
    arb_init(w);
    arb_init(z);

    arf_set_d(&interval->a, a);
    arf_set_d(&interval->b, b);

    flint_printf("interval: "); arf_interval_printd(interval, 15); flint_printf("\n");
    flint_printf("maxdepth = %wd, maxeval = %wd, maxfound = %wd, low_prec = %wd\n",
		 maxdepth, maxeval, maxfound, low_prec);

    num = arb_calc_isolate_roots(&blocks, &info, function,
				 params, interval, maxdepth, maxeval, maxfound, low_prec);

    for (i = 0; i < num; i++) {
	if (info[i] != 1) {
	  if (1) //(arb_calc_verbose)
	    // else the following compile error: Undefined symbols for architecture arm64: "___emutls_v.arb_calc_verbose", referenced from: __ZN7Onsager3runEv in ccSUd03Y.o ld: symbol(s) not found for architecture arm64
	      {
		flint_printf("unable to count roots in ");
		arf_interval_printd(blocks + i, 15);
		flint_printf("\n");
	    }
	    found_unknown++;
	    continue;
	}

	found_roots++;

	if (!refine)
	  continue;

	if (arb_calc_refine_root_bisect(t, function, params, blocks + i, 5, low_prec) != ARB_CALC_SUCCESS) {
	    flint_printf("warning: some bisection steps failed!\n");
	  }

	if (1) //(arb_calc_verbose)
	  {
	    flint_printf("after bisection 1: ");
	    arf_interval_printd(t, 15);
	    flint_printf("\n");
	  }

	if (arb_calc_refine_root_bisect(blocks + i, function, params, t, 5, low_prec) != ARB_CALC_SUCCESS) {
	    flint_printf("warning: some bisection steps failed!\n");
	  }

	if (1)// (arb_calc_verbose)
	  {
	    flint_printf("after bisection 2: ");
	    arf_interval_printd(blocks + i, 15);
	    flint_printf("\n");
	  }

	arf_interval_get_arb(v, t, high_prec);
	arb_calc_newton_conv_factor(C, function, params, v, low_prec);

	arf_interval_get_arb(w, blocks + i, high_prec);
	if (arb_calc_refine_root_newton(z, function, params, w, v, C, 10, high_prec) != ARB_CALC_SUCCESS) {
	    flint_printf("warning: some newton steps failed!\n");
	  }

	flint_printf("refined root (%wd/%wd):\n", i, num);
	arb_printn(z, digits + 2, 0);
	flint_printf("\n\n");
      }

    flint_printf("---------------------------------------------------------------\n");
    flint_printf("Found roots: %wd\n", found_roots);
    flint_printf("Subintervals possibly containing undetected roots: %wd\n", found_unknown);
    flint_printf("Function evaluations: %wd\n", eval_count);

    for (i = 0; i < num; i++)
      arf_interval_clear(blocks + i);
    flint_free(blocks);
    flint_free(info);

    arb_clear(interactions.J);                                         //^^^^^^
    arb_clear(interactions.K);                                         //^^^^^^

    
    arf_interval_clear(t);
    arf_interval_clear(interval);
    arf_clear(C);
    arb_clear(v);
    arb_clear(w);
    arb_clear(z);
    flint_cleanup();
  }

[saraedum]

Any help is appreciated, I hope this is an appropiate place to request for such.

Yes, I think this is an appropriate place.

It seems you are describing three separate issues here. It's maybe helpful to split this into three GitHub issues with shorter self-contained reproducers.

Segfault

I don't understand which program produces the segfault you are describing? The complete program you posted at the bottom? Could you post the entire source code for the program that segfaults for you?

As a general comment, I see that you are using functions starting with an _ in some places. Why was that necessary? These usually have more complicated contracts (and I imagine that that contract is not always spelled out entirely in the documentation.)

Undefined Symbols

How did you install Arb? Are you sure that your headers match your library build?

some newton steps failed

I would have to debug into this to see what's going on. The code you posted shows that behavior if I perfom the following replacement?

 a = 0;
-b = 5;
+b = 5;

Which version of Arb is this?

---

btw., if you wrap your code in

```c
C soruce code here
```

it's a bit easier to read since we get syntax highlighting.

[reneruhr]

To differentiate the problems, I installed the current version of arb (and flint,mpfr,gmp) on a different machine so that "if(arb_calc_verbose)" no longer fails.

I would have to debug into this to see what's going on. The code you posted shows that behavior if I perfom the following replacement?
a = 0;
-b = 5;
+b = 5;

I don't see the difference.

On the new machine, b can be arbitrary, but a needs to be bigger than 0.3 and no roots are found. So I presume that the code doesn't really do what it is supposed to do thanks to lack of understanding of the underlying types.

The original code is the whole code except for the main program

int main(){
   run();
}

I rewrote the program without underscore as below.

The problem is that the arb example "real_roots" uses the root-finding function slong arb_calc_isolate_roots( function )
gets as input a function of type int func(arb_ptr out, const arb_t inp, void * param, slong order, slong prec)

Since

An arb_poly_t is defined as an array of length one of type arb_poly_struct, permitting an arb_poly_t to be passed by reference.

type arb_ptr¶ Alias for arb_struct *, used for vectors of numbers.

I define a global arb_poly_struct out_inner and take out = out_inner[0].coeffs.
This gives the following code, which does not work. (it doesn't find any zeros, even if I change the polynomial to "x-1".)

#include "arb.h"
#include "arb_calc.h"

  slong eval_count = 0;

  struct Interactions
  {
    arb_t J;
    arb_t K;
  };

  Interactions interactions;
  
  arb_poly_t out_inner;

   int critical_temp(arb_ptr out, const arb_t inp, void * params, slong order, slong prec)
  {
    auto inter_ptr = static_cast<Interactions*>(params);
    auto J = inter_ptr->J;
    auto K = inter_ptr->K;

    arb_poly_init(out_inner);
    arb_poly_t x,y;
    arb_poly_init(x);
    arb_poly_init(y);
    arb_poly_one(x);
    arb_poly_one(y);
    arb_poly_shift_left(x,x,1); 
    arb_poly_shift_left(y,y,1);

    arb_poly_scalar_mul(x,x,J,prec);
    arb_poly_scalar_mul(y,y,K,prec);
    
    arb_t two;
    arb_init(two);
    arb_set_ui(two,2);
    arb_poly_scalar_mul(x,x,two,prec);
    arb_poly_scalar_mul(y,y,two,prec);
    arb_clear(two);

    arb_poly_t one;
    arb_poly_init(one);
    arb_poly_one(one);

    arb_poly_sinh_series(x, x, order, prec);
    arb_poly_sinh_series(y, y, order, prec);
    arb_poly_mullow(out_inner, x, y, order, prec);
    arb_poly_sub(out_inner, out_inner, one, prec); // sub = "difference between A and B. Is this A-B or B-A?

    out = out_inner[0].coeffs;
    
    arb_poly_clear(one); 
    arb_poly_clear(x);
    arb_poly_clear(y);
    return 0;
  }

  void run()
  {
    arf_interval_ptr blocks;
    arb_calc_func_t function;
    int * info;
    void * params;
    int param1;
    slong digits, low_prec, high_prec, i, num, found_roots, found_unknown;
    slong maxdepth, maxeval, maxfound;
    int refine;
    double a, b;
    arf_t C;
    arf_interval_t t, interval;
    arb_t v, w, z;

    auto inter = Interactions();
    arb_init(inter.J);
    arb_init(inter.K);
    arb_one(inter.J);
    arb_one(inter.K);

    params = (void*)(&inter);
    
    function = critical_temp;

    a = -10;
    b = 10;  

    refine = 1;
    digits = 10;
    maxdepth = 30;
    maxeval = 100000;
    maxfound = 100000;
    low_prec = 30;

    high_prec = digits * 3.32192809488736 + 10;
    found_roots = 0;
    found_unknown = 0;

    arf_init(C);
    arf_interval_init(t);
    arf_interval_init(interval);
    arb_init(v);
    arb_init(w);
    arb_init(z);

    arf_set_d(&interval->a, a);
    arf_set_d(&interval->b, b);

    flint_printf("interval: "); arf_interval_printd(interval, 15); flint_printf("\n");
    flint_printf("maxdepth = %wd, maxeval = %wd, maxfound = %wd, low_prec = %wd\n",
		 maxdepth, maxeval, maxfound, low_prec);

    num = arb_calc_isolate_roots(&blocks, &info, function,
				 params, interval, maxdepth, maxeval, maxfound, low_prec);

    for (i = 0; i < num; i++) {
	if (info[i] != 1) {
	  if (arb_calc_verbose)
	      {
		flint_printf("unable to count roots in ");
		arf_interval_printd(blocks + i, 15);
		flint_printf("\n");
	    }
	    found_unknown++;
	    continue;
	}

	found_roots++;

	if (!refine)
	  continue;

	if (arb_calc_refine_root_bisect(t, function, params, blocks + i, 5, low_prec) != ARB_CALC_SUCCESS) {
	    flint_printf("warning: some bisection steps failed!\n");
	  }

	if(arb_calc_verbose)
	  {
	    flint_printf("after bisection 1: ");
	    arf_interval_printd(t, 15);
	    flint_printf("\n");
	  }

	if (arb_calc_refine_root_bisect(blocks + i, function, params, t, 5, low_prec) != ARB_CALC_SUCCESS) {
	    flint_printf("warning: some bisection steps failed!\n");
	  }

	if(arb_calc_verbose)
	  {
	    flint_printf("after bisection 2: ");
	    arf_interval_printd(blocks + i, 15);
	    flint_printf("\n");
	  }

	arf_interval_get_arb(v, t, high_prec);
	arb_calc_newton_conv_factor(C, function, params, v, low_prec);

	arf_interval_get_arb(w, blocks + i, high_prec);
	if (arb_calc_refine_root_newton(z, function, params, w, v, C, 10, high_prec) != ARB_CALC_SUCCESS) {
	    flint_printf("warning: some newton steps failed!\n");
	  }

	flint_printf("refined root (%wd/%wd):\n", i, num);
	arb_printn(z, digits + 2, 0);
	flint_printf("\n\n");
      }

    flint_printf("---------------------------------------------------------------\n");
    flint_printf("Found roots: %wd\n", found_roots);
    flint_printf("Subintervals possibly containing undetected roots: %wd\n", found_unknown);
    flint_printf("Function evaluations: %wd\n", eval_count);

    for (i = 0; i < num; i++)
      arf_interval_clear(blocks + i);
    flint_free(blocks);
    flint_free(info);

    arb_clear(interactions.J);                                         //^^^^^^
    arb_clear(interactions.K);                                         //^^^^^^

    
    arf_interval_clear(t);
    arf_interval_clear(interval);
    arf_clear(C);
    arb_clear(v);
    arb_clear(w);
    arb_clear(z);
    flint_cleanup();
  }

int main()
{
  run();
}

[saraedum]

To differentiate the problems, I installed the current version of arb (and flint,mpfr,gmp) on a different machine so that "if(arb_calc_verbose)" no longer fails.

I would have to debug into this to see what's going on. The code you posted shows that behavior if I perfom the following replacement?
[…]

I don't see the difference.

Sorry, I meant replacing b=5 with b=10.

On the new machine, b can be arbitrary, but a needs to be bigger than 0.3 and no roots are found. So I presume that the code doesn't really do what it is supposed to do thanks to lack of understanding of the underlying types.

I see. I guess one would have to separate this into more digestible chunks to verify that the code is doing what it's supposed to do.

While that would probably not be too hard to do, maybe a question is what you are actually trying to achieve. (Were you just playing with arb to learn about its capabalities or is computing roots of these sinh functions something that's relevant to you as such.)

[saraedum]

Also, note that the documentation for arb_calc.h says:

This code should be considered experimental.

So you might well have uncovered a bug here somehow.