Arb v. 2.23.0: doubt about arb_calc_isolate_roots

[Rudi Gaelzer]

Dear Fredrik.

Using Arb v. 2.23.0, I'm trying to implement the function arb_calc_isolate_root to find intervals containing the roots of a given continuous real function.

I was under the impression that the algorithm uses basically bisection to isolate the intervals, but in my specific implementation I noticed that the function func (as a arb_calc_func_t struct) is called with the argument order (as a slong) equal to 1 or even 2.

Is this correct?  the algorithm can employ the first and second derivatives of func to isolate the individual roots?

But in this case, I also noticed in the example code real_roots.c that individual functions take specific action depending on order.  For example, the function sin_x contains the declaration/assignment 

int xlen = FLINT_MIN(2, order);

Could you please explain a bit about the role played by the argument order?

Thanks.

[Stephen Crowley]

if order = 0 then the function isnt required to be continuous, if order = 1, then the result at position 0 should be continuous, and of order>1, for instance 2, means the value of the function goes in spot 0 and its derivative goes in index 1, if order=3 then the function is required to calculate the fuinctions value and 1st two derivatives. it should thrown an error if its asked to compute a derivative it cant compute

[Rudi Gaelzer]

Dear Stephen.  Thanks for the reply.

Indeed, by trial-and-error I figured all that out.  I just don't know why a bisection method would need to evaluate the derivative(s) of func, but I confess that I did not take time to study the source code.

Anyway, I wrote the function that evaluates all the derivatives up to order - 1 and managed to implement the method to find all roots within a given interval.

In my tests, I verified that the refinement with arb_calc_refine_root_newton frequently fails, and so I'm doing all refinement for the individual roots using arb_calc_refine_root_bisect.  Of course,  this is slow, but right now I value accuracy more than speed.  

However, my final aim is the evaluation of the complex roots of an analytic function, and thus I will next use acb_calc_integrate to bisect the complex roots in a given region using  the argument principle.

[Stephen Crowley]

No problem Rudi, glad you figured it out. I did study the source extensively and ported it to Java, you might find this useful

https://github.com/crowlogic/arb4j/blob/7da4e80504a7fc410a16ac251a719ce9797ea276/src/main/java/arb/functions/real/RealFunction.java#L334

I also encountered many situations where the Newton iteration fails, but the code  I linked to here has just as robust abilities as the C

implementation and I think I fixed a small glitch or two here or there and perhaps introduced some more. I do have quite a bit of unit tests,

it's just for my purposes I'm not so focused on root finding any more as I am on implementing integral transforms

Best regards

--Stephen

--

---
You received this message because you are subscribed to a topic in the Google Groups "flint-devel" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/flint-devel/jfhvk38OLmE/unsubscribe.
To unsubscribe from this group and all its topics, send an email to flint-devel...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/7a811d3d-4b83-4681-9855-a62a50519197n%40googlegroups.com.

[Stephen Crowley]

More pertinent to your question, this code shows how the bisection is first did and then a Newton step is attempted

package arb;

import static arb.utensils.Utensils.println;

import static java.lang.String.format;

import static java.lang.System.out;

import arb.functions.real.RealFunction;

public class RealRootInterval extends

FloatInterval

{

@Override

public RealRootInterval set(FloatInterval interval)

{

return (RealRootInterval) super.set(interval);

}

@Override

public String toString()

{

return String.format("[%s, %s]=%s", getA(), getB(), status);

}

public RealRootInterval(double left,

double right)

{

super(

left,

right);

assert left < right;

}

public RealRootInterval()

{

}

public RealRootInterval(Float left,

Float right)

{

super(

left,

right);

assert left.compareTo(right) <= 0;

}

public RealRootInterval(FloatInterval interval)

{

this(

interval.getA(),

interval.getB());

}

public RealRootInterval swap(RealRootInterval u)

{

RealRootInterval r = (RealRootInterval) super.swap(u);

RootStatus status = u.status;

u.status = r.status;

r.status = status;

return r;

}

public RootStatus status = RootStatus.RootUnknown;

/**

* Splits the interval. Chiefly called by

* {@link RealFunction#recursivelyLocateRoots(Roots, RealRootInterval, int, int, int, RootLocatorOptions)}

*

* @param func

* @param found

* @param asign

* @param bsign

* @param depth

* @param config

* @return

*/

public boolean split(RealFunction func, Roots found, int asign, int bsign, int depth, RootLocatorOptions config)

{

if (RealFunction.verbose)

{

println(String.format("%s.split(func=%s, found=%s, asign=%s, bsign=%s, depth=%s)",

this,

func,

found,

asign,

bsign,

depth));

}

try ( RealRootInterval L = new RealRootInterval(); RealRootInterval R = new RealRootInterval();)

{

int msign = func.calculatePartition(L, R, this, config.prec);

found.evals++;

if (msign == 0)

{

return false;

}

func.recursivelyLocateRoots(found, L, asign, msign, depth + 1, config);

func.recursivelyLocateRoots(found, R, msign, bsign, depth + 1, config);

}

return true;

}

public RootStatus determineRootStatus(Roots found, RealFunction func, int asign, int bsign, int prec)

{

status = RootStatus.RootUnknown;

try ( Real t = Real.newVector(2); Real x = new Real())

{

func.evaluate(getReal(x, prec), 1, prec, t.get(0));

found.evals++;

if (t.isPositive() || t.isNegative())

{

status = RootStatus.NoRoot;

}

else

{

if ((asign < 0 && bsign > 0) || (asign > 0 && bsign < 0))

{

Real firstDerivative = func.evaluate(x, 2, prec, t).get(1);

found.evals++;

if (firstDerivative.isFinite() && !firstDerivative.containsZero())

{

status = RootStatus.RootLocated;

}

}

}

}

return status;

}

public static enum RefinementResult

{

ImpreciseInput,

Success,

NoConvergence

}

public static final int FLINT_BITS = 64;

/**

* Increases the precision of the root

*

* @param f

* @param highPrec

* @param w 3-vector of Taylor jet

* @param v input/output result

* @param u

* @param convergenceRegion output

* @param convergenceFactor output

* @param rootInterval

* @return

*/

public Real refine(RealFunction func,

int lowPrec,

int digits,

Real w,

Real v,

RealRootInterval u,

RealRootInterval convergenceRegion,

Float convergenceFactor,

boolean verbose)

{

if (RealFunction.verbose)

{

println(format("%s.refine(func=%s, lowPrec=%d, digits=%d))", this, func, lowPrec, digits));

}

int highPrec = (int) (digits * 3.32192809488736 + 10);

if (status == RootStatus.NoRoot)

{

return null;

}

if (bisectAndRefine(func, v, convergenceRegion, u, 5, lowPrec) != RefinementResult.Success)

{

if (RealFunction.verbose)

{

System.out.println("Bisection failed");

}

}

else

{

if (convergenceRegion.bisectAndRefine(func, v, this, u, 5, lowPrec) != RefinementResult.Success)

{

if (RealFunction.verbose)

{

System.out.println("second Bisection failed");

}

}

}

convergenceRegion.getReal(v, highPrec);

func.getNewtonConvergenceFactor(v, w, highPrec, convergenceFactor);

RefinementResult refinementResult = func.refineRootNewton(v,

getReal(w, highPrec),

convergenceFactor,

10,

highPrec,

verbose);

if (refinementResult != RefinementResult.Success)

{

if (RealFunction.verbose)

{

out.format("Warning: some newton steps failed: refinementResult=%s\n", refinementResult);

}

}

if (RealFunction.verbose)

{

System.out.println("Refined root: " + v);

}

return v;

}

public RefinementResult

bisectAndRefine(RealFunction func, Real v, RealRootInterval t, RealRootInterval u, int iters, int prec)

{

if (RealFunction.verbose)

{

println(String.format("%s.bisectAndRefine(func=%s, t=%s, iters=%s, prec=%s)\n", this, func, t, iters, prec));

}

try ( Real m = new Real();)

{

m.setMid(getA());

int asign = func.evaluate(m, 1, prec, v).sign();

m.setMid(getB());

int bsign = func.evaluate(m, 1, prec, v).sign();

/* must have proper sign changes */

if (asign == 0 || bsign == 0 || asign == bsign)

{

return RefinementResult.ImpreciseInput;

}

else

{

for (int i = 0; i < iters; i++)

{

int msign = func.calculatePartition(t, u, this, prec);

/*

* TODO: handle the case where the value at the midpoint is actually zero even

* though it can't be distinguished via sign comparison zero, this should be

* just checking isZero before returning NoConvergence.. right?

*/

if (msign == 0)

{

return RefinementResult.NoConvergence;

}

swap(msign == asign ? u : t);

}

}

}

return RefinementResult.Success;

}

}

[Rudi Gaelzer]

Very interesting.  Unfortunately, I don't use Java, but Modern Fortran instead.  

In fact, I intend to employ the highly accurate quadrature evaluation provided by Arb and adapt it into an existing Fortran code that locates the roots in the complex plane using the argument principle.

I find Arb so useful that in the future I'd like to try and do an implementation in Fortran.  However, I'm still trying to figure out how to implement the Flint structs, particularly those that employ unsigned integers, since those are not supported by the Fortran standard.

[Stephen Crowley]

Sounds fun . Enjoy😁

--

---
You received this message because you are subscribed to a topic in the Google Groups "flint-devel" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/flint-devel/jfhvk38OLmE/unsubscribe.
To unsubscribe from this group and all its topics, send an email to flint-devel...@googlegroups.com.

To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/d168b5c1-f2af-4f0d-9dcc-8bf590b71e55n%40googlegroups.com.

[Akhil Akkapelli]

I'm intrested in integration of Fortran functions using functions in Arb library. I'm eager to follow the progress of your implementation of Arb in Fortran.