/*
 * 5799-WZQ (C) COPYRIGHT IBM CORPORATION 1987
 * LICENSED MATERIALS - PROPERTY OF IBM
 * REFER TO COPYRIGHT INSTRUCTIONS FORM NUMBER G120-2083
 */
/* $Header:l10.c 12.0$ */
/* $ACIS:l10.c 12.0$ */
/* $Source: /ibm/acis/usr/src/lib/libc/ca/gen/RCS/l10.c,v $ */

#ifndef lint
static char *rcsid = "$Header:l10.c 12.0$";
#endif

/* APPROXIMATE BASE-10 LOGARITHM OF DOUBLE FLOATING-POINT NUMBER

This function returns the approximate base-10 logarithm of its 
IEEE P754 double floating-point argument, disguised as an unsigned 
long two-element vector. The result is an integer ranging between 
-324 and 308, corresponding to arguments 0x0000000000000001 and 
0x7fefffffffffffff, respectively. This function does not handle 
zero, infinite, or NAN arguments, and does not test to exclude 
them. Negative numbers are treated as positive. 

The algorithm used is modified from one given in "Contributions 
to a Proposed Standard for Binary Floating-Point Arithmetic", 
chapter 7, section 2.3, "A Poor Man's Logarithm", by Jerome 
Coonen, Ph. D. dissertation, University of California, Berkeley, 
California, 1984. It differs by providing for subnormal numbers.
Subnormal numbers are unnormalized numbers with exponent 0x000.

Coonen's analysis shows that the result is the accurate floor of 
the common logarithm of the argument, with an average error of 
about 0.017. The result may be perhaps too low by 1. Rapidity is 
gained by casting the algorithm completely in terms of integer 
operations. */

extern _mul211();
int _l10 (w)
unsigned long w[2];

{long int b, t5, t6, x, t11, p[2], z;
 union {unsigned long ul[2]; long l;} u;
 unsigned long pp[2];

/* Two cases must be handled, normalized and subnormalized 
numbers. Subnormalized numbers are detected by masking out all 
but the exponent field; if what remains is 0, the argument is 
subnormal. */
         
 if (w[0] & 0x7ff00000)

/* [Normal number.] The only thing needed is to unbias the 
exponent, after insuring that the sign bit is 0. */ 

 {x = (w[0] & 0x7fffffff) - 0x3ff00000;}

 else

/* [Subnormal number.] We convert this to a normalized number by 
counting how many leading zeros it has, and incidentally 
left-adjusting it, then losing the leading bit, and positioning the 
remaining bits as if they were a normalized mantissa, and lastly 
by converting the exponent 0x000 to an unbiassed one, taking into 
account the argument's number of leading zeros, less the combined 
width of the sign and exponent fields. We copy w into pp in order 
to avoid affecting the global variable w.*/

  {pp[0] = w[0] & 0x7fffffff; pp[1] = w[1];
   b = _shlun(pp); /* b counts leading zeros, and as a side-effect,
                    left-adjusts pp. */
   u.ul[0] = (pp[0] & 0x7fffffff) >> 11;
   t5 = 0xc01 - (b - 12);
   t6 = t5 << 20;
   x = t6 + u.l;}

/* The multiplier, an integer form of the common logarithm of 2, 
is increased by one if the argument is negative. */

 t11 = 19728 + (x < 0);

/* Multiply the unbiassed number's leading 32 bits by the 
approximation to the logarithm of 2. */

 _mul211(x, t11, p);

/* Shift the product right to take its floor, giving the desired 
 result. */

 z = p[0] >> 4; 

/* [Exit.] */

 return (z);}

