conversion of a denorm float -> double is busted

[Uncle George]

> conversion of a denorm float -> double is busted on the alpha.
>
> from the attached program , the iteration at I=-126 fails to properly
> convert the first denormal float number into a corresponding double
> ( the double is not denormal ). The exponent for the double should have
> been 896 ( 1 less than the previous iteration ( at i=-125), and frac=0.
> /gat
>
> does anyone know if the sparc, or ppc fails also ?
>
> alpha results:
> I=-126  prev=2.3509887016446e-38, d=1.1754943508223e-38
> Value:  f=800000, s=(0,0), exp=(1,1), frac=(0,0)
> V\=2.0  f=400000, s=(0,0), exp=(0,0), frac=(4194304,400000)
> Conversion: f=8000000000000, s=(0,0), exp=(0,0),
> frac=(2251799813685248,8000000000000)
> reverted: f=0, s=(0,0), exp=(0,0), frac=(0,0)
>  Conversion failed
>
> i386 results:
> I=-126  prev=2.3509887016446e-38, d=1.1754943508223e-38
> Value:  f=800000, s=(0,0), exp=(1,1), frac=(0,0)
> V\=2.0  f=400000, s=(0,0), exp=(0,0), frac=(4194304,400000)
> Conversion: f=0, s=(0,0), exp=(896,380), frac=(0,0)
> reverted: f=400000, s=(0,0), exp=(0,0), frac=(4194304,400000)
>
>

#include <stdio.h>

union _FP_UNION_D
{
  double flt;
  long   l;
  struct {
    unsigned long long frac : 64-12;
    unsigned exp  : 11 ;
    unsigned sign : 1;
  } bits __attribute__((packed));
} ud;
union _FP_UNION_S
{
  float flt;
  int    i;
  struct {
    unsigned int frac : 32-9;
    unsigned exp  : 8 ;
    unsigned sign : 1;
  } bits __attribute__((packed));
} uf;
int main()
{
        int errout;
#if defined(__alpha__) && 0
        unsigned long fpcr;
        asm( "excb; mf_fpcr %0" : "=f"(fpcr));

printf(" fpcr = %p\n", fpcr );
#if 1
        fpcr |= (1L << 47);     /* DNOD  Not avail in 21264(a) */
        fpcr &= ~(1L << 48 );   /* DNZ idenormal Operands to Zero */
        fpcr |= (1L << 49);     /* DNOD  Not avail in 21264(a) */
        fpcr &= ~(1L << 60 );   /* UNDZ underflow to zero OFF */
        fpcr |= (1L << 61);     /* UNFD underflow disable ON */
/*      fpcr |= (1L << 62);      UNFD */

        asm volatile ( "mt_fpcr %0 ; excb" : : "f"(fpcr));
#endif
        asm( "excb; mf_fpcr %0" : "=f"(fpcr));
printf(" fpcr = %p\n", fpcr );
#endif
        underflow();
} 
void printF( float, char *);
void printD( double, char *);
 
int underflow()
{
        int dsb=32-9, dme=-126;

        int expLimit = dme - dsb + 1;
        float divisor = 2.0;
        float  d = 1.0, e;
        float  prev=0;
        int i;
        double conversion;

        for (i=0; i>=expLimit; i--) {
                printf("I=%d\t prev=%16.14g, d=%16.14g\n",i, prev, d);
                printF( d, "Value:\t");
                prev = d;
                if ((d /= divisor) >= prev) {
                        printf("Oops. growing larger ? \n");
                        return;
                }
                else if (d == 0)
                        break;
                conversion = d;         /* double <= float */
                e = conversion;         /* float <= double */
                printF( d, "V\\=2.0\t");
                printD( conversion, "Conversion:");
                printF( e, "reverted:");
                if ( d != e ){          /* Still same ?    */
                   printf(" Conversion failed \n");
                   break;
                }
        };
printf("i = %d, expLimit= %d\n",i,expLimit);

        if (i > expLimit)
                printf("Oops, didnt do all denormal numbers\n");

        return;
}
void
printF( float d, char *msg )
{
                uf.flt = d;
printf("%s\tf=%x, s=(%d,%x), exp=(%d,%x), frac=(%d,%x)\n", 
                msg,
                uf.i,
                uf.bits.sign, uf.bits.sign,
                uf.bits.exp, uf.bits.exp,
                uf.bits.frac, uf.bits.frac );   
}

void
printD( double d, char *msg)
{
                ud.flt  = d;
printf("%s\tf=%lx, s=(%d,%x), exp=(%d,%x), frac=(%ld,%lx)\n", 
                msg,
                ud.l,
                ud.bits.sign, ud.bits.sign,
                ud.bits.exp, ud.bits.exp,
                ud.bits.frac, ud.bits.frac );           
}