From Newsgroup: comp.lang.c
Hans Bezemer <
the.beez.speaks@gmail.com> wrote or quoted:
/* increase mantissa */
while ((zen -> mantissa < mant) && (zen -> exponent > exp))
{ /* while between limits */
zen -> mantissa *= 10; --(zen -> exponent);
} /* multiply by 10 */
Above, you mix pre-comments (referring to what follows) with
post-comments (referring to what precedes) giving them /exactly
the same format/, which makes them hard to read.
I'd prefer something like,
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <limits.h>
typedef struct {
long long mantissa;
long exponent;
} me;
/* integer version of pow() */
long power (long base, long exp)
{
long long a = 1LL;
do {
if (exp & 1) a *= base;
exp /= 2; base *= base;
} while (exp != 0);
return a; }
/* normalize mantissa/exponent */
void normalize10 (me* zen, long long mant, long exp)
{
/* save sign of mantissa */
char msign = (zen -> mantissa) < 0LL;
/* get absolute value of mantissa */
zen -> mantissa = llabs (zen -> mantissa);
/* increase mantissa */
while ((zen -> mantissa < mant) && (zen -> exponent > exp))
{ /* while between limits */
zen -> mantissa *= 10; --(zen -> exponent); /* multiply by 10 */ }
/* correct mantissa sign if required */
if (msign) zen -> mantissa = -(zen -> mantissa); }
/* convert from radix 10 to radix 2 */
void me10tome2 (me* zen)
{ /* save signs of mantissa and exponent */
char msign = (zen -> mantissa) < 0LL;
char esign = (zen -> exponent) < 0L;
long exp2;
long me10; /* radix 10 multiplication factor */
long me2; /* radix 2 multiplication factor */
if (zen -> exponent) { /* if exponent is non-zero */
zen -> mantissa = llabs (zen -> mantissa);
zen -> exponent = labs (zen -> exponent); /* get absolute value of mantissa and exponent */
exp2 = (((zen -> exponent) * 10) + 2) / 3; /* convert exponent to radix 2 */
me10 = power (10, zen -> exponent);
me2 = 1 << exp2; /* calculate multiplication factors */
if (esign) { /* if negative exponent */
zen -> exponent = -exp2;
zen -> mantissa = (zen -> mantissa * me2) / me10;
} else { /* if positive exponent */
zen -> exponent = exp2;
zen -> mantissa = (zen -> mantissa * me10) / me2;
}
/* correct mantissa sign if required */
if (msign) zen -> mantissa = -(zen -> mantissa); }}
int main (int argc, char** argv)
{
me z = { 1000LL, -4L }; /* 1000E-4 */
normalize10 (&z, LONG_MAX, -9); /* normalize mantissa and exponent */
me10tome2 (&z); /* convert to radix 2 */
/* check result */
printf("Verification:\t%lld * 2^(%ld) = %f\n",
z.mantissa, z.exponent, z.mantissa * pow(2.0, z.exponent)); }
. Here's what that former web search-engine said,
| This algorithm is a variant of a known, standard approach for base
| conversion in floating-point software libraries rather than a
| brand-new mathematical discovery.
|
| It implements a fixed-precision, scaling-factor conversion method
| using log approximations.
|
| Why It Is a Known Approach
|
| The core logic relies on two standard computer science techniques:
|
| Logarithmic Scaling Estimation The line exp2 = (((zen -> exponent) *
| 10) + 2) / 3; approximates log2(10)rCa~rCa3.3333. Multiplying by 10 and
| dividing by 3 yields 3.3333, which is a classic integer approximation
| used to guess the required binary exponent size for a given decimal
| exponent.
|
| Exponentiation by Squaring Your power function is the standard,
| well-known exponentiation-by-squaring algorithm used to compute powers
| in O(logrCan) time.
|
| Limitations and Risks of This Implementation
|
| While the general concept is standard, this specific implementation
| has severe edge-case flaws that prevent it from being used in
| production environments:
|
| Integer Overflow The power(10, zen->exponent) call will rapidly
| overflow standard 64-bit integer limits. 10^19 is the absolute maximum
| a 64-bit unsigned integer can hold. If the input exponent is larger
| than 18, me10 will overflow and corrupt the calculation.
|
| Bit-Shift Overflow The line me2 = 1 << exp2; uses standard integer
| shifting. In C, shifting a 32-bit signed integer by more than 31 bits
| (or a 64-bit integer by more than 63 bits) causes undefined behavior
| and severe overflow.
|
| Precision Loss Standard library conversion functions (like strtod or
| printf internals) use arbitrary-precision arithmetic (bignums) or
| specialized algorithms like GayrCOs algorithm, Dragon4, or Grisu to
| maintain exact bit-level precision without risking integer overflow.
.
--- Synchronet 3.22a-Linux NewsLink 1.2