Re: compile-time conversion of floating-point expressions to long longs

Re: compile-time conversion of floating-point expressions to long longs

[Nelson H. F. Beebe]

Jim Cromie <jim.cromie@gmail.com> asks on Sun, 07 May 2006 13:41:41 -0400
why gcc won't evaluate an initialization involving frexp() at compile time.

In general, compilers only permit initializations using constant
expressions that involve no external references.  This is mandated
by the ISO C Standards.  From ISO/IEC 9899:1999 (E):

>> ...
>>      6.6 Constant expressions
>>
>>      Syntax
>>
>> 1	       constant-expression:
>> 		       conditional-expression
>>
>>      Description
>>
>> 2    A constant expression can be evaluated during translation rather
>>      than runtime, and accordingly may be used in any place that a
>>      constant may be.
>>
>>      Constraints
>>
>> 3    Constant expressions shall not contain assignment, increment,
>>      decrement, function-call, or comma operators, except when they
                   ^^^^^^^^^^^^^ NOTE
>>      are contained within a subexpression that is not evaluated.
>>
>> 4    Each constant expression shall evaluate to a constant that is in
>>      the range of representable values for its type.
>> ...

Recall that it is common practice today to use shared libraries for
widely-used code.  Thus, an initialization like "const float root13 =
sqrt(13);" at outer level could have different meanings on different
runs or different systems of the same architecture, depending on which
version of the math library were used, and the accuracy of the
implementation of the sqrt() function.  Also, if the compiler is
cross-compiling for another architecture, it would not have access to
that other platform's run-time library, and thus could not compute the
expression.  Thus, the value is x is not really a constant as far as
machine implementations are concerned.

Historically, across CPU architectures, the meaning of "float", and
its precision, could change, although all modern machines use the IEEE
754 system, so that is no longer the case.  For example, I have a DEC
PDP-10 on my desktop, with a 36-bit word, and three bits more
precision in float than my other systems.

If you want to have initializations with mathematical expressions, the
correct way is to do one of two things:

(1) For C99 compilation, start with an accurate value of the
    expression evaluated in higher precision, perhaps via a symbolic
    algebra system with user-specifiable precision, and then express
    that as a hexadecimal constant.  For the above example, I might
    then write:

	const float root13 = +0x1.cd82b446159f360fedeccf37f9e5p+1F; /*	sqrt(13) */

    Some compilers will complain about the extra digits, but that
    doesn't matter: you'll get a correct-to-last-bit value for root13.

(2) For either C99 or pre-C99, express the constant as the sum of an
    exactly-representable value and a small correction.  For example,

	const float root13 = 945173.0F / 262144.0F + 2.4168214111681192212674704961299214e-06;

    Both 945173.0f and 262144.0F are exactly representable in the IEEE
    754 32-bit floating-point format, which has 24 bits of precision.
    Older floating-point systems always had at least that much
    precision for the float data type.  Importantly, 262144.0F is a
    power of two (2 to the 18), so the division is EXACT, since in the
    absence of underflow or overflow, multiplication or division by a
    power of the base simply adjusts the exponent, with changing the
    significand.   Thus, the first term is computed exactly at compile
    time.  The second term is a small correction, so even if the
    decimal-to-binary conversion is subject to small errors, they will be
    well beyond the precision of the result, and the stored constant
    computed at compile time will be correct to the last bit.

Technique (2) is common in carefully-written mathematical software.
As C99 compilers become more common, technique (1) will be used more
often.

It certainly is NOT sufficient to simply express the number in decimal
and assign it with, e.g.,

	const float root13 = 3.60555127546398929311922126747049613F;

The accuracy of decimal-to-binary conversion varies across systems,
and the result is that the constant could get values that differ by
one or two units in the last place from the exact value obtained by
either of the techniques shown above.

-------------------------------------------------------------------------------
- Nelson H. F. Beebe                    Tel: +1 801 581 5254                  -
- University of Utah                    FAX: +1 801 581 4148                  -
- Department of Mathematics, 110 LCB    Internet e-mail: beebe@math.utah.edu  -
- 155 S 1400 E RM 233                       beebe@acm.org  beebe@computer.org -
- Salt Lake City, UT 84112-0090, USA    URL: http://www.math.utah.edu/~beebe/ -
-------------------------------------------------------------------------------

Re: compile-time conversion of floating-point expressions to long longs

[Jim Cromie]

Nelson H. F. Beebe wrote:
> Jim Cromie <jim.cromie@gmail.com> asks on Wed, 03 May 2006 12:17:38 -0400:
>
>   
>>> ...
>>> Are there any macros that can tear into a floating point number and
>>> pull out the exponent and mantissa ?
>>> ...
>>>       
>
> Yes, the C89 and C99 ISO C Standards include ldexp(x,n), which forms x
> * 2**n, and "f = frexp(x,&n);", which returns the fraction as a value
> in [1/2,1) (or 0 if x is zero) and the exponent of 2 in n.  Both are
> exact, and can be implemented reasonably efficiently.
>   

yes, frexp() is almost what I want.
heres basically what I tried to do with it:

#include <stdio.h>
#include <math.h>
#include <float.h>

/* evaluate these constants at compile-time ?!? */
float f = 3.14159 * 2 * 100;
int myexp;
double mant = frexp(f,&myexp);

int main(int c, char** v)
{
    printf("float %g mantissa %f exp %d \n", f, mant, myexp);
}


as Id feared, it wouldnt compile;
$ make fp-comp
cc -g    fp-comp.c   -o fp-comp
fp-comp.c:37: error: initializer element is not constant
make: *** [fp-comp] Error 1

Im guessing its referring to myexp, the addy of which is passed.
This seems to foul up the constant-folder / compile-time evaluator,
but maybe this just wont work on more fundamental grounds.

If the problem is the 2nd arg being modified during evaluation,
then Id consider trying to split the function into 2 separate ones,
(since its compile-time, I dont care about evaluation cost of 2 calls),
but again, I think im oversimplifying.

IOW, I suspect that the compiler just wont call the function for me at 
compile-time,
and fold the constant into the code for me.  Id love to hear differently.
It also sounds like a rather unique requirement that would be tossed 
with scorn.

tia
jimc

Re: compile-time conversion of floating-point expressions to long longs

[Jim Cromie]

Nelson H. F. Beebe wrote:
> Jim Cromie <jim.cromie@gmail.com> asks on Wed, 03 May 2006 12:17:38 -0400:
>
>   
>>> ...
>>> Are there any macros that can tear into a floating point number and
>>> pull out the exponent and mantissa ?
>>> ...
>>>       
>
> Yes, the C89 and C99 ISO C Standards include ldexp(x,n), which forms x
> * 2**n, and "f = frexp(x,&n);", which returns the fraction as a value
> in [1/2,1) (or 0 if x is zero) and the exponent of 2 in n.  Both are
> exact, and can be implemented reasonably efficiently.
>   
thank you for the thorough answer.

yes, frexp() is almost what I want.
heres basically what I tried to do with it:

#include <stdio.h>
#include <math.h>
#include <float.h>

/* evaluate these constants at compile-time ?!? */
float f = 3.14159 * 2 * 100;
int myexp;
double mant = frexp(f,&myexp);

/* this confirms gcc's unwillingness to run libmath routines at 
compile-time */
float ff = ldexp(2,8);

int main(int c, char** v)
{
    printf("float %g mantissa %f exp %d \n", f, mant, myexp);
}


as Id feared, it wouldnt compile;
$ make fp-comp
cc -g    fp-comp.c   -o fp-comp
fp-comp.c:37: error: initializer element is not constant
make: *** [fp-comp] Error 1

(the code above is slightly edited for clarity - err line numbers dont 
match, but result does ;-(

I also tried these additions to compile cmd:
   -fkeep-static-consts  -fmerge-constants  -fmerge-all-constants
    -fsingle-precision-constant

So it seems that the compiler just wont call the function for me at 
compile-time,
and fold the constant into the code for me.  Id love to hear differently.

tia
jimc

Re: compile-time conversion of floating-point expressions to long longs

[Nelson H. F. Beebe]

Jim Cromie <jim.cromie@gmail.com> asks on Wed, 03 May 2006 12:17:38 -0400:

>> ...
>> Are there any macros that can tear into a floating point number and
>> pull out the exponent and mantissa ?
>> ...

Yes, the C89 and C99 ISO C Standards include ldexp(x,n), which forms x
* 2**n, and "f = frexp(x,&n);", which returns the fraction as a value
in [1/2,1) (or 0 if x is zero) and the exponent of 2 in n.  Both are
exact, and can be implemented reasonably efficiently.

Two related, and useful, functions are

(1) modf(x,*y), which splits x into an integral part stored in y as an
    exactly-representable floating-point number, and returns the
    fractional part as the function value.  Both have the same sign as x.

(2) fmod(x,y), which returns the remainder of x divided by y.

All four of these functions are EXACT; no rounding errors whatsoever
are possible.  Note that for fmod(), the computation involves finding
the integer close to x/y, and that can be very large: in IEEE 754
128-bit arithmetic, it is a number of nearly 9900 decimal digits.
Fortunately, it is possible to implement fmod() iteratively, without
access to high-precision floating-point.

ldexp() and frexp() were in C in Unix Version 7 in the late 1970s; the
other two are more recent.  You can safely use all four of them on
modern systems.  Pre-C99 implementations may lack the float and long
double counterparts.  However, glibc on GNU/Linux offers all twelve
variants.

-------------------------------------------------------------------------------
- Nelson H. F. Beebe                    Tel: +1 801 581 5254                  -
- University of Utah                    FAX: +1 801 581 4148                  -
- Department of Mathematics, 110 LCB    Internet e-mail: beebe@math.utah.edu  -
- 155 S 1400 E RM 233                       beebe@acm.org  beebe@computer.org -
- Salt Lake City, UT 84112-0090, USA    URL: http://www.math.utah.edu/~beebe/ -
-------------------------------------------------------------------------------

compile-time conversion of floating-point expressions to long longs

[Jim Cromie]

Im working on a program that cannot use floating-point, so it must
do shift and integer multiply, with the 2 constants being determined by
the need to maintain maximum precision.

Id like to bury the math in macros, and to have it all computed
at compile-time, and reduced to constants.  GCC only solutions are fine.

The obvious way apparently involves floating point.. 
(all those f* instructions - pun not intended)

int main(int c, char** v)
{
    float f = 3.14159 * 2 * 100;
    long long run = (long long) f;
    printf("float %g cast-to-long %lld \n", f, run);
}
$ fp-comp
float 628.318 cast-to-long 628

int main(int c, char** v)
{
 80483b4:       55                      push   %ebp
 80483b5:       89 e5                   mov    %esp,%ebp
 80483b7:       83 ec 18                sub    $0x18,%esp
 80483ba:       83 e4 f0                and    $0xfffffff0,%esp
 80483bd:       b8 00 00 00 00          mov    $0x0,%eax
 80483c2:       83 c0 0f                add    $0xf,%eax
 80483c5:       83 c0 0f                add    $0xf,%eax
 80483c8:       c1 e8 04                shr    $0x4,%eax
 80483cb:       c1 e0 04                shl    $0x4,%eax
 80483ce:       29 c4                   sub    %eax,%esp
        float f = 3.14159 * 2 * 100;
 80483d0:       b8 5a 14 1d 44          mov    $0x441d145a,%eax
 80483d5:       89 45 f4                mov    %eax,0xfffffff4(%ebp)
        long long run = (long long) f;
 80483d8:       d9 45 f4                flds   0xfffffff4(%ebp)
 80483db:       d9 7d ee                fnstcw 0xffffffee(%ebp)
 80483de:       66 8b 45 ee             mov    0xffffffee(%ebp),%ax
 80483e2:       b4 0c                   mov    $0xc,%ah
 80483e4:       66 89 45 ec             mov    %ax,0xffffffec(%ebp)
 80483e8:       d9 6d ec                fldcw  0xffffffec(%ebp)
 80483eb:       df 7d f8                fistpll 0xfffffff8(%ebp)
 80483ee:       d9 6d ee                fldcw  0xffffffee(%ebp)


The actual conversion Im after is cycles -> nanoseconds,
where cycles is pulled from a Time Stamp Counter, or similar.


Are there any macros that can tear into a floating point number
and pull out the exponent and mantissa ?  Arch-specific is ok, as long 
as they exist.