0

I made a function that uses Newton's method to compute integer roots for positive real numbers. It works for every case I have tried, except the square root of 2.

I have a repeat loop that iterates Newton's method as long as X[n] does not equal X[n-1]. For the square root of 2, I eventually get X[n] = 1.4142135624 and X[n-1] = 1.4142135624, yet the loop continues despite X[n] equaling X[n-1]. How do I fix this problem? Is it a bug within the Delphi language? (This was made in Delphi 7 - SE)

function IntPower(Const Base : Real; Const Power : Integer):Real; //Integer power function
var
  K : Integer;
  Output : Real;
begin
  Output := 1;
  If Power >= 0 then
    For K := 1 to Power do
    begin
      Output := Output * Base;
    end
  else
    For K := -1 downto Power do
    begin
      Output := Output / Base;
    end;
  Result := Output;
end;

function IntRoot(Base : Real; Root : Integer):Real;
var
  K : Integer;
  Output, Greater, Less, Previous : Real;
begin
  if (Root < 0) AND (Base = 0) then
  begin
    Result := 0;
    Exit;
  end;
  if (Root = 1) OR (Base = 1) then
  begin
    Result := Base;
    Exit;
  end;
  K := 2;
  Previous := 0;
  if Root < 0 then
  begin
    Root := Root * (-1);
    Base := 1/Base;
  end;
  //BEGINNING OF INITIAL GUESS//
  if Base < 1 then
  begin
    Greater := 1;
    Less := 0;
    While Less = 0 do
    begin
      if IntPower(1/K,Root) > Base then
        Greater := 1/K
      else if IntPower(1/K,Root) < Base then
        Less := 1/K
      else
      begin
        Result := 1/K;
        Exit;
      end;
      Inc(K);
    end;
  end
  else if Base > 0 then
  begin
    Greater := 0;
    Less := 1;
    While Greater = 0 do
    begin
      if IntPower(K,Root) > Base then
        Greater := K
      else if IntPower(K,Root) < Base then
        Less := K
      else
      begin
        Result := K;
        Exit;
      end;
      Inc(K);
    end;
  end;
  Output := (Greater+Less)/2;
  //END OF INITIAL GUESS//
  Repeat
    Previous := Output;
    Output := Previous - (IntPower(Previous,Root)-base)/((Root)*IntPower(Previous,Root-1));
  Until(Output = Previous);
  Result := Output;
end;

begin
  Writeln(FloatToStr(IntRoot(2,2));
  Readln;
end.

Yes, I know that I made the nth root of 0 = 0 for n < 0 and I know there is a problem with 0^n for n < 0.

Justin
  • 91
  • 1
  • 5
  • 6
    Use e.g. `...until SameValue(Output, Previous);`. `SameValue` function is from Math unit. – TLama Jul 28 '15 at 21:08
  • 3
    See [How to compare double in delphi?](http://stackoverflow.com/q/6106119/576719). – LU RD Jul 28 '15 at 21:17
  • @TLama You also have to decide on how to choose a tolerance. That requires analysis of the algorithm. There's a lot more to this than you let on. – David Heffernan Jul 29 '15 at 06:31
  • @David, true, I should have rather show a call with the `Epsilon` parameter. My point was essentially about the used equality operator. – TLama Jul 29 '15 at 07:15
  • @LURD FWIW I think that the canonical dupe is this: http://stackoverflow.com/questions/588004/is-floating-point-math-broken – David Heffernan Jul 29 '15 at 07:16
  • @DavidHeffernan, yes I saw that one too. But I find closing with a delphi tag more comfortable, maybe it is wrong thinking? – LU RD Jul 29 '15 at 08:04
  • @LURD Personally I don't regard this as a delphi issue – David Heffernan Jul 29 '15 at 08:07

0 Answers0