[BACK]Return to extract-dbl.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp

Diff for /OpenXM_contrib/gmp/Attic/extract-dbl.c between version 1.1.1.1 and 1.1.1.2

version 1.1.1.1, 2000/09/09 14:12:16 version 1.1.1.2, 2003/08/25 16:06:00
Line 1 
Line 1 
 /* __gmp_extract_double -- convert from double to array of mp_limb_t.  /* __gmp_extract_double -- convert from double to array of mp_limb_t.
   
 Copyright (C) 1996, 1999, 2000 Free Software Foundation, Inc.  Copyright 1996, 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
   
 This file is part of the GNU MP Library.  This file is part of the GNU MP Library.
   
Line 30  MA 02111-1307, USA. */
Line 30  MA 02111-1307, USA. */
 #define _GMP_IEEE_FLOATS 0  #define _GMP_IEEE_FLOATS 0
 #endif  #endif
   
   #define BITS_IN_MANTISSA 53
   
 /* Extract a non-negative double in d.  */  /* Extract a non-negative double in d.  */
   
 int  int
 #if __STDC__  
 __gmp_extract_double (mp_ptr rp, double d)  __gmp_extract_double (mp_ptr rp, double d)
 #else  
 __gmp_extract_double (rp, d)  
      mp_ptr rp;  
      double d;  
 #endif  
 {  {
   long exp;    long exp;
   unsigned sc;    unsigned sc;
   mp_limb_t manh, manl;  #ifdef _LONG_LONG_LIMB
   #define BITS_PER_PART 64        /* somewhat bogus */
     unsigned long long int manl;
   #else
   #define BITS_PER_PART GMP_LIMB_BITS
     unsigned long int manh, manl;
   #endif
   
   /* BUGS    /* BUGS
   
      1. Should handle Inf and NaN in IEEE specific code.       1. Should handle Inf and NaN in IEEE specific code.
      2. Handle Inf and NaN also in default code, to avoid hangs.       2. Handle Inf and NaN also in default code, to avoid hangs.
      3. Generalize to handle all BITS_PER_MP_LIMB >= 32.       3. Generalize to handle all GMP_LIMB_BITS >= 32.
      4. This lits is incomplete and misspelled.       4. This lits is incomplete and misspelled.
    */     */
   
     ASSERT (d >= 0.0);
   
   if (d == 0.0)    if (d == 0.0)
     {      {
       rp[0] = 0;        MPN_ZERO (rp, LIMBS_PER_DOUBLE);
       rp[1] = 0;  
 #if BITS_PER_MP_LIMB == 32  
       rp[2] = 0;  
 #endif  
       return 0;        return 0;
     }      }
   
Line 72  __gmp_extract_double (rp, d)
Line 72  __gmp_extract_double (rp, d)
     union ieee_double_extract x;      union ieee_double_extract x;
     x.d = d;      x.d = d;
     exp = x.s.exp;      exp = x.s.exp;
 #if BITS_PER_MP_LIMB == 64  #if BITS_PER_PART == 64         /* generalize this to BITS_PER_PART > BITS_IN_MANTISSA */
     manl = (((mp_limb_t) 1 << 63)      manl = (((mp_limb_t) 1 << 63)
             | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));              | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
     if (exp == 0)      if (exp == 0)
Line 87  __gmp_extract_double (rp, d)
Line 87  __gmp_extract_double (rp, d)
           }            }
         while ((mp_limb_signed_t) manl >= 0);          while ((mp_limb_signed_t) manl >= 0);
       }        }
 #else  #endif
   #if BITS_PER_PART == 32
     manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);      manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
     manl = x.s.manl << 11;      manl = x.s.manl << 11;
     if (exp == 0)      if (exp == 0)
Line 104  __gmp_extract_double (rp, d)
Line 105  __gmp_extract_double (rp, d)
         while ((mp_limb_signed_t) manh >= 0);          while ((mp_limb_signed_t) manh >= 0);
       }        }
 #endif  #endif
   #if BITS_PER_PART != 32 && BITS_PER_PART != 64
     You need to generalize the code above to handle this.
   #endif
     exp -= 1022;                /* Remove IEEE bias.  */      exp -= 1022;                /* Remove IEEE bias.  */
   }    }
 #else  #else
Line 140  __gmp_extract_double (rp, d)
Line 144  __gmp_extract_double (rp, d)
           }            }
       }        }
   
     d *= MP_BASE_AS_DOUBLE;      d *= (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2)));
 #if BITS_PER_MP_LIMB == 64  #if BITS_PER_PART == 64
     manl = d;      manl = d;
 #else  #else
     manh = d;      manh = d;
     manl = (d - manh) * MP_BASE_AS_DOUBLE;      manl = (d - manh) * (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2)));
 #endif  #endif
   }    }
 #endif  #endif /* IEEE */
   
   sc = (unsigned) exp % BITS_PER_MP_LIMB;    /* Up until here, we have ignored the actual limb size.  Remains
        to split manh,,manl into an array of LIMBS_PER_DOUBLE limbs.
     */
   
     sc = (unsigned) (exp + 64 * GMP_NUMB_BITS) % GMP_NUMB_BITS;
   
   /* We add something here to get rounding right.  */    /* We add something here to get rounding right.  */
   exp = (exp + 2048) / BITS_PER_MP_LIMB - 2048 / BITS_PER_MP_LIMB + 1;    exp = (exp + 64 * GMP_NUMB_BITS) / GMP_NUMB_BITS - 64 * GMP_NUMB_BITS / GMP_NUMB_BITS + 1;
   
 #if BITS_PER_MP_LIMB == 64  #if LIMBS_PER_DOUBLE == 2
   #if GMP_NAIL_BITS == 0
   if (sc != 0)    if (sc != 0)
     {      {
       rp[1] = manl >> (BITS_PER_MP_LIMB - sc);        rp[1] = manl >> (GMP_LIMB_BITS - sc);
       rp[0] = manl << sc;        rp[0] = manl << sc;
     }      }
   else    else
Line 168  __gmp_extract_double (rp, d)
Line 177  __gmp_extract_double (rp, d)
       exp--;        exp--;
     }      }
 #else  #else
     if (sc > GMP_NAIL_BITS)
       {
         rp[1] = manl >> (GMP_LIMB_BITS - sc);
         rp[0] = (manl << (sc - GMP_NAIL_BITS)) & GMP_NUMB_MASK;
       }
     else
       {
         if (sc == 0)
           {
             rp[1] = manl >> GMP_NAIL_BITS;
             rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS) & GMP_NUMB_MASK;
             exp--;
           }
         else
           {
             rp[1] = manl >> (GMP_LIMB_BITS - sc);
             rp[0] = (manl >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK;
           }
       }
   #endif
   #endif
   
   #if LIMBS_PER_DOUBLE == 3
   #if GMP_NAIL_BITS == 0
   if (sc != 0)    if (sc != 0)
     {      {
       rp[2] = manh >> (BITS_PER_MP_LIMB - sc);        rp[2] = manh >> (GMP_LIMB_BITS - sc);
       rp[1] = (manl >> (BITS_PER_MP_LIMB - sc)) | (manh << sc);        rp[1] = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc));
       rp[0] = manl << sc;        rp[0] = manl << sc;
     }      }
   else    else
Line 181  __gmp_extract_double (rp, d)
Line 214  __gmp_extract_double (rp, d)
       rp[0] = 0;        rp[0] = 0;
       exp--;        exp--;
     }      }
   #else
     if (sc > GMP_NAIL_BITS)
       {
         rp[2] = (manh >> (GMP_LIMB_BITS - sc));
         rp[1] = ((manh << (sc - GMP_NAIL_BITS)) |
                  (manl >> (GMP_LIMB_BITS - sc + GMP_NAIL_BITS))) & GMP_NUMB_MASK;
         if (sc >= 2 * GMP_NAIL_BITS)
           rp[0] = (manl << sc - 2 * GMP_NAIL_BITS) & GMP_NUMB_MASK;
         else
           rp[0] = manl >> (2 * GMP_NAIL_BITS - sc) & GMP_NUMB_MASK;
       }
     else
       {
         if (sc == 0)
           {
             rp[2] = manh >> GMP_NAIL_BITS;
             rp[1] = ((manh << GMP_NUMB_BITS - GMP_NAIL_BITS) | (manl >> 2 * GMP_NAIL_BITS)) & GMP_NUMB_MASK;
             rp[0] = 0;
             exp--;
           }
         else
           {
             rp[2] = (manh >> (GMP_LIMB_BITS - sc));
             rp[1] = (manh >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK;
             rp[0] = ((manh << (GMP_NUMB_BITS - GMP_NAIL_BITS + sc))
                      | (manl >> (GMP_LIMB_BITS - (GMP_NUMB_BITS - GMP_NAIL_BITS + sc)))) & GMP_NUMB_MASK;
           }
       }
   #endif
   #endif
   
   #if LIMBS_PER_DOUBLE > 3
     /* Insert code for splitting manh,,manl into LIMBS_PER_DOUBLE
        mp_limb_t's at rp.  */
     if (sc != 0)
       {
         /* This is not perfect, and would fail for GMP_LIMB_BITS == 16.
            The ASSERT_ALWAYS should catch the problematic cases.  */
         ASSERT_ALWAYS ((manl << sc) == 0);
         manl = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc));
         manh = manh >> (GMP_LIMB_BITS - sc);
       }
     {
       int i;
       for (i = LIMBS_PER_DOUBLE - 1; i >= 0; i--)
         {
           rp[i] = manh >> (BITS_PER_LONGINT - GMP_LIMB_BITS);
           manh = ((manh << GMP_LIMB_BITS)
                   | (manl >> (BITS_PER_LONGINT - GMP_LIMB_BITS)));
           manl = manl << GMP_LIMB_BITS;
         }
     }
 #endif  #endif
   
   return exp;    return exp;
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.1.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>