Date:	Mon, 17 Feb 1992 16:08:16 +0100
From:	Dirk Zabel <dzabel AT cs DOT tu-berlin DOT de>
To:	djgpp AT sun DOT soe DOT clarkson DOT edu
Subject:	emu387 patch
Status:	O

Hello,

I recently read a few messages from this mailing-list indicating that
some people use djgcc for doing numerical computations.  If they use
the 80387 emulation package, they could get into trouble due to a
problem with the log-function.

You can see the problem using the sample calculator (samples/hexcalc)
and entering "log10(10)" or "log(e)". The answer does not make sense.

After some digging around I finally found the reason(s):
Logarithms are computed using the FYL2X instruction, which takes two
arguments (st(0) and st(1)) and replaces them with log2(st(0))*st(1).

The corresponding function (fyl2x() in e16.cc) uses the well-known
power series 
   log((1+x)/(1-x))/2 = x + x**3 / 3 + x**5 /5 + ...              (1)
    
In order to exploit this equation, first x is computed from
   arg = (1+x)/(1-x)  ==> x = (arg-1)/(arg+1)
The next step SHOULD have been to compute (a finite part of) the above
series (1), multiplying the result with 2*log2(e) (getting log2(arg))
and then with st(1) in order to follow the specification of the FYL2X
instruction. Instead, st(1) is multiplied with x BEFORE building the
sum. This is the main bug.

But, after I corrected this bug the emulator still did not give good
results - especially not for "large" arguments. And this is why:
(1) assumes |x| < 1 to converge. This is always true if arg has
a legal value (arg > 0 ), but large values of arg result in |x|
being nearly 1 (very small value have the same consequence; arg and 1/arg
result in the same |x|). So for arg being too large or too small
the series converges quite slow. Fortunately, there is a simple solution:
exploit the floating-point format of 
arg = 2**(exponent-EXP_BIAS) * significand.
Therefore, log2(arg) = 
  log2(2**(exponent-EXP_BIAS) *significand) = 
                   exponent-EXP_BIAS + log2(significand)
This can be done defining the new variable 
   save_exponent = arg.exponent-EXP_BIAS
and setting arg.exponent=EXP_BIAS.  I call the new value of arg arg'.   
So I get log2(arg) = save_exponent+log2(arg') .
                   
If one further exploits that |x| is the same for arg and 1/arg,
one can minimize |x| by dividing arg' by 2 if arg' > sqrt(2) 
(that is, by incrementing save_exponent and decrementing arg.exponent). 
This ensures that sqrt(2) >= arg' > 1/sqrt(2) and therefore 
|x| <= 1/(3+2*sqrt(x)) < 0.172 . 
This in turn helps to determine the  number of terms of (1) which 
has to be taken into account: if (0.172)**n / n < 2**-64,
no more terms can influence the result. The resulting limit is n=23, i.e.
we have to sum up the terms for n=1, 3, ... 23.

I have applied the corrections and got a working log-Funktion. (I also
added the "invalid operation" - exception for arguments <= 0). 

	Dirk Zabel


Dirk Zabel	<dzabel AT cs DOT tu-berlin DOT de>
Technische Universitaet Berlin


PS:  This patch applies to djgpp 1.05 (I think this is still the latest version)
PPS: djgpp is a great piece of software!!!

---------cut here -------------- cut here -----------------------------------
*** e16.cc	Wed Sep 11 14:48:40 1991
--- e16.new	Thu Dec 12 15:51:58 1991
***************
*** 32,54 ****
  {
    if (empty())
      return;
!   reg frac2, sum, div, term, pow, temp;
!   r_sub(st(), CONST_1, frac2);
!   r_add(st(), CONST_1, div);
!   r_div(frac2, div, sum);
!   r_mul(sum, sum, frac2);
!   r_mul(sum, st(1), pow);
!   for (long i=3; i<15; i+=2)
    {
!     r_mul(pow, frac2, temp);
      r_mov(temp, pow);
!     r_mov(&i, div);
!     r_div(temp, div, term);
      r_add(term, sum, temp);
      r_mov(temp, sum);
    }
    r_div(sum, CONST_LN2, temp);
    temp.exp++;
    r_mov(temp, st(1));
    st().tag = TW_E;
    top++;
--- 32,80 ----
  {
    if (empty())
      return;
!   reg z, x, nom, denom, xsquare, term, temp, sum, pow;
!   long exponent;
!   reg CONST_SQRT2 = { SIGN_POS, TW_V, EXP_BIAS, 0xf9de6000, 0xb504f333 };
! 
!   r_mov(st(), z);
!   if ((z.tag != TW_V) || (z.sign != SIGN_POS)) {
!     return exception(EX_I); // not valid, zero or negative
!   }  
!   exponent = (long)(z.exp - EXP_BIAS);
!   z.exp=EXP_BIAS;
!   if (compare(z, CONST_SQRT2) == COMP_A_GT_B) {
!     (z.exp)--;
!     exponent++;
!   }
! 
!   r_sub(z, CONST_1, nom);
!   r_add(z, CONST_1, denom);
!   r_div(nom, denom, x);
!   r_mov(x, pow);
!   r_mov(x, sum);
!   r_mul(x, x, xsquare);
!   
!   for (long i=3; i<25; i+=2)
    {
!     r_mul(pow, xsquare, temp);
      r_mov(temp, pow);
! 
!     r_mov(&i, denom);
!     r_div(pow, denom, term);
! 
      r_add(term, sum, temp);
      r_mov(temp, sum);
    }
    r_div(sum, CONST_LN2, temp);
    temp.exp++;
+   if (exponent) {
+     r_mov(&exponent, term);
+     r_add(term, temp, sum);
+   } else {
+     r_mov(temp, sum);
+   }
+ 
+   r_mul(sum, st(1), temp);
    r_mov(temp, st(1));
    st().tag = TW_E;
    top++;
---------cut here -------------- cut here -----------------------------------

- Raw text -