0

I ran into a problem with Ted Ersek's RootSearch (https://library.wolfram.com/infocenter/Demos/4482/) while solving the problem described here. RootSearch sometimes hangs after reporting that a number is less than the smallest number allowed in the Wolfram Mathematica. The reason turned out to be the internal function Ulp2. It entered into an infinite loop if its first argument x1 was equal to 0. (That's right: zero with a dot). Here is the corrected version of this function:

Ulp2=Function[{x1,x2},
      Module[{eps,sgn=Sign2[x2-x1],tol},
        (*Print["RootSearch Ulp2: sgn=",sgn,", x1=",x1,", x2=",x2];*)
        Which[
          (*test1*)sgn==0,
          (*Print["RootSearch Ulp2-0: sgn=",sgn,", x1=",x1,", x2=",x2];*)
          Message[Ulp2::fail];Throw[False],
          (*test2*)MachineNumberQ[x1],
          (*Print["***RootSearch Ulp2-test2: MachineNumberQ[x1]=",MachineNumberQ[x1]];*)
          (*tol=$MachineEpsilon*Max[Abs[x1],$MinMachineNumber]/2;*)
          tol=Max[$MachineEpsilon*Abs[x1],$MinMachineNumber];
          (*Print["***RootSearch Ulp2-test2: tol=",NumberForm[tol,{40,36}]];*)
          While[(eps=Plus@@{x1+sgn*tol,-x1})==0.0,tol*=2];
          (*Print["***RootSearch Ulp2-test2: eps=",eps];*)
          Abs[eps],
          (*test3*)True,
          (*Print["***RootSearch Ulp2-test3: sgn=",sgn];*)
          Block[{$MinPrecision=0},
            tol=SetPrecision[16*10^-Accuracy[x1],20];

        While[(eps=Plus@@{x1+sgn*tol,-x1};Precision[eps]<=0),
          tol*=2];
        SetPrecision[Abs[eps],20]
        ]
      ]
    ]];

Diagnostics came from the line

tol=$MachineEpsilon*Max[Abs[x1],$MinMachineNumber]/2;

The tol paramter was evaluated to 0 in this line so the next line enters to endless cycle:

While[(eps=Plus@@{x1+sgn*tol,-x1})==0.0,tol*=2];

I have changed definition of tol as follows:

tol=Max[$MachineEpsilon*Abs[x1],$MinMachineNumber];

which seems to me more adequate. Can someone propose a more robust solution?

Igor Kotelnikov
  • 749
  • 3
  • 12
  • 2
    If RootSearch is not part of Mathematica itself, then please include a link, otherwise everyone interested will have to use a search engine. Also, I understand the tag "bugs" to be for bugs in Mathematica itself, but I am just a random user. – user293787 Jul 03 '22 at 06:54
  • 2
    I think questions of the form "Is this correct" will rarely have a definite answer in floating point matters. The original code allows for $MachineEpsilon*$MinMachineNumber which produces a warning message, so I agree that that does not make much sense. But numerical code always has edge cases, even more so if it is by a third party. So if your solution works, and you can validate it, why not run with it. – user293787 Jul 03 '22 at 08:11
  • 2
    I don't think this is a good fit for this site. It appears to be not a question, but a bug report. It should go directly to the package's author. Have you contacted him before posting here? – Szabolcs Jul 03 '22 at 09:36
  • 3
    It looks to me that instead of insisting on using RootSearch, which has not been touched in a long while, you should be using any of the other alternatives presented in this question. – J. M.'s missing motivation Jul 03 '22 at 15:38
  • @J.M.: Solutions based on Plot and its associated documented and undocumented functions are slow for more or less complex equations. I have my own version of RootSearch with all known bugs fixed. – Igor Kotelnikov Jul 05 '22 at 02:06
  • @Szabolcs: I tried to contact the package's author in the past, but he did not react. – Igor Kotelnikov Jul 05 '22 at 02:09

0 Answers0