Frank wrote:
function Log1P (x: Real): Real; asmname 'log1p';
According to the manpage:
: CONFORMING TO : BSD
Oops. My manpage says the same. OTOH, info libc doesn't say anything (it usually does for non-ISO C89 functions), the header says
#if defined __USE_MISC || defined __USE_XOPEN_EXTENDED || defined __USE_ISOC9X
("MISC" means the common intersection of BSD and System V), and I've succesfully compiled a simple test program even with -D__STRICT_ANSI__. This seems like quite a mess to me.
So this function might not be available everywhere, and we need a replacement for other systems. Can you please provide one (written in Pascal)? I can deal with autoconf then (i.e., checking if the function is available and using the replacement if not).
Here it is. The C library function is certainly preferable when present -- it is likely to have a hardware support, so it is faster (although I tried to optimize the Pascal implementation for speed, not code clarity).
--------------------------------------- var ExpEpsReal: Integer; Log1PCoefs: array [0 .. BitSizeOf (Real) div 3] of Real; DummyMantissa: LongReal;
function Log1P (X: Real) = Y: Real; var I, N: Integer; Z, U: Real; begin if (X <= -0.5) or (X >= 1) then Y := Ln (1 + X) else begin U := X / (2 + X); Z := Sqr (U); if Z <= MinReal then Y := X else begin SplitReal (Z, N, DummyMantissa); N := ExpEpsReal div N; Y := 0; for I := N downto 0 do Y := Y * Z + Log1PCoefs[I]; Y := Y * U end end end;
var I: Integer;
to begin do begin SplitReal (EpsReal, ExpEpsReal, DummyMantissa); Dec (ExpEpsReal); for I := 0 to BitSizeOf (Real) div 3 do Log1PCoefs[I] := 2 / (I * 2 + 1) end; -----------------------------------------------
Before I look closer into the routine (I'm quite short of time at the moment), one remark:
To save your time, here is some mathematics behind my Complex_Arctan reimplementation:
Let z = x+iy. Arctan (z) is defined as -i/2 Ln((1+iz)/(1-iz)). Thus
Re(arctan z) = 1/2 arg ((1+iz)/(1-iz)) = = 1/2 arg ( (1+iz)(1+iz') / |1-iz|^2 ) = = 1/2 arg ((1+iz)(1+iz')) = = 1/2 arg (1+i(z+z')-|z|^2) = = 1/2 arg (1-x^2-y^2+2ix) = = 1/2 Arctan2 (2x, (1+y)(1-y)-x^2)
where z' denotes the conjugate of z, and
Im(arctan z) = -1/2 ln |(1+iz)/(1-iz)| = = -1/2 ln |1+iz|/|1-iz| = = 1/2 ln |1-iz|/|1+iz| = = 1/2 ln \sqrt{((1+y)^2 + x^2) / ((1-y)^2 + x^2)} = = 1/4 ln (((1+y)^2 + x^2) / ((1-y)^2 + x^2)) = = 1/4 ln (1 + 4y / ((1-y)^2 + x^2))
The function deals with |y| instead of y to save some "if"s; it's clear from the equations above that arctan(z') = (arctan z)'.
The "if c <= MinReal" branch (needed to avoid some overflows) is executed when |y|=1 and x^2 <= MinReal. In this case,
1/4 ln (((1+|y|)^2 + x^2) / ((1-|y|)^2 + x^2) = 1/4 ln ((4+x^2) / x^2) = [not exactly, but recall that x^2 << 1] 1/4 ln (4 / x^2) = 1/2 ln (2/x)
The Pascal version of log1p() uses the power series
ln ((1+u) / (1-u)) = 2 \sum_{n=0}^{\infty} u^{2n+1} / (2n+1) (|u|<1)
together with ln (1+x) = ln ((1+u) / (1-u)), where u = x/(2+x).
Emil Jerabek