[fpc-devel] Changed Math functions in FPC development trunk 2.7

[info at wolfgang-ehrhardt.de]

A year ago there was a post from Jonas Maebe on this list
(http://www.mail-archive.com/fpc-devel@lists.freepascal.org/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 (arjan.vanDijk at User.METAIR.WAU.NL) }
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;