A year ago there was a post from Jonas Maebe on this list
(http://www.mail-archive.com/[email protected]/msg28483.html)
asking about the possibility to reuse parts of AMath/DAMath in FPC.

Today I downloaded the development source files from
(ftp://ftp.freepascal.org/pub/fpc/snapshot/trunk/source/fpc.zip)
and found very close source code similarity between the changed
(compared with FPC264) math.pp routines and AMath routines.
See the list below.

Although, as Jonas wrote, zlib is compatible with FPC, there is the
need to give references and/or mark the changes. Unfortunately this
is not done.

I am quite sure that these routines are 'inspired' by AMath.
If this is the case, I would appreciate reference to 'AMath by W.Ehrhardt'
or so, as it is done in the FPC264 source for Arjan van Dijk.

Best wishes

Wolfgang Ehrhardt


hypot
---------------------------------------------------------------------------
AMATH
function hypot(x,y: extended): extended;
  {-Return sqrt(x*x + y*y)}
begin
  x := abs(x);
  y := abs(y);
  if x>y then hypot := x*sqrt(1.0+sqr(y/x))
  else if x>0.0 then hypot := y*sqrt(1.0+sqr(x/y)) {here y >= x > 0}
  else hypot := y;                                 {here x=0}
end;


FPC-Dev
function hypot(x,y : float) : float;
  begin
    x:=abs(x);
    y:=abs(y);
    if (x>y) then
      hypot:=x*sqrt(1.0+sqr(y/x))
    else if (x>0.0) then
      hypot:=y*sqrt(1.0+sqr(x/y))
    else
      hypot:=y;
  end;

FPC 264
function hypot(x,y : float) : float;

  begin
     hypot:=Sqrt(x*x+y*y)
  end;


arcsin and arccos
---------------------------------------------------------------------------
AMATH
function arcsin(x: extended): extended;
  {-Return the inverse circular sine of x, |x| <= 1}
begin
  {basic formula arcsin(x) = arctan(x/sqrt(1-x^2))}
  arcsin := arctan2(x, sqrt((1.0-x)*(1.0+x)))
end;

function arccos(x: extended): extended;
  {-Return the inverse circular cosine of x, |x| <= 1}
begin
  {basic formula arccos(x) = arctan(sqrt(1-x^2)/x))}
  if abs(x)=1.0 then begin
    if x<0.0 then arccos := Pi else arccos := 0.0;
  end
  else arccos := arctan2(sqrt((1.0-x)*(1.0+x)),x)
end;

FPC-Dev
function arcsin(x : float) : float;
begin
  arcsin:=arctan2(x,sqrt((1.0-x)*(1.0+x)));
end;

function Arccos(x : Float) : Float;
begin
  if abs(x)=1.0 then
    if x<0.0 then
      arccos:=Pi
    else
      arccos:=0
  else
    arccos:=arctan2(sqrt((1.0-x)*(1.0+x)),x);
end;

FPC 264
{ ArcSin and ArcCos from Arjan van Dijk ([email protected]) }
function arcsin(x : float) : float;
begin
  if abs(x) > 1 then InvalidArgument
  else if abs(x) < 0.5 then
    arcsin := arctan(x/sqrt(1-sqr(x)))
  else
    arcsin := sign(x) * (pi*0.5 - arctan(sqrt(1 / sqr(x) - 1)));
end;


function Arccos(x : Float) : Float;
begin
  arccos := pi*0.5 - arcsin(x);
end;


lnxp1 = ln1p
---------------------------------------------------------------------------
AMath (up to V1.74 of 04.Oct.2013, V1.75 introduces an additional branch)

function ln1p(x: extended): extended;
  {-Return ln(1+x), accurate even for x near 0}
var
  y: extended;
const
  x0 = 4.0;
begin
  if x >= x0 then ln1p := ln(1.0 + x)
  else begin
    y := 1.0 + x;
    {The following formula is more accurate than Goldberg [3], Theorem 4. The}
    {Taylor series f(x) = f(x0) + (x-x0)*f'(x0) + .. for f(x) = ln(1+x) gives}
    {ln1p(x) = ln(1+x0) + (x-x0)/(1+x0) = ln(y) + (x-(y-1))/y, with y = 1+x0.}
    if y=1.0 then ln1p := x
    else ln1p := ln(y) + (x-(y-1.0))/y;
  end
end;

FPC-Dev
function lnxp1(x : float) : float;
  var
    y: float;
  begin
    if (x>=4.0) then
      lnxp1:=ln(1.0+x)
    else
      begin
        y:=1.0+x;
        if (y=1.0) then
          lnxp1:=x
        else
          lnxp1:=ln(y)+(x-(y-1.0))/y;
      end;
  end;

FPC 264
function lnxp1(x : float) : float;
  begin
     if x<-1 then
       InvalidArgument;
     lnxp1:=ln(1+x);
  end;






_______________________________________________
fpc-devel maillist  -  [email protected]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-devel

Reply via email to