Changes between Version 2 and Version 7 of ReplacingGMPNotes/TheCurrentGMPImplementation


Ignore:
Timestamp:
Jan 6, 2007 6:41:21 PM (8 years ago)
Author:
p_tanski
Comment:

add notes for Storage Manager, Special Functions sections

Legend:

Unmodified
Added
Removed
Modified
  • ReplacingGMPNotes/TheCurrentGMPImplementation

    v2 v7  
    66 
    77[[Image(GMP_interface.jpg)]] 
     8 
     9=== GMP and the Storage Manager === 
     10 
     11The most insidious and unique feature of the GHC implementation with GMP is memory management.  GHC uses the special memory functions in GMP ( `mp_set_memory_functions( (*alloc)(), (*realloc)(), (*dealloc)()` ) to let GMP use GHC's own garbage collector.  The memory functions are implemented in [[GhcFile(rts/sm/Storage.c)]] as: 
     12{{{ 
     13          static void* stgAllocForGMP (size_t size_in_bytes); 
     14          static void  stgDeallocForGMP (void *ptr STG_UNUSED,  
     15                                         size_t size STG_UNUSED); 
     16          static void* stgReallocForGMP (void *ptr,  
     17                                         size_t old_size,  
     18                                         size_t new_size); 
     19}}} 
     20These special allocation functions bring most GMP memory usage into the GHC heap but do not seem that efficient otherwise.  (I could be wrong --PDT.)  Allocation uses the internal `allocate()` interface, so no garbage collection is performed during a GMP operation.  Note that GMP operations may allocate more memory themselves.  The memory allocated is a simple array of words (`W_`), rounded up to a whole number. 
     21 
     22=== Special Functions === 
     23 
     24GHC adds its own functions for string conversion, least common multiple (lcm) and conversion to and from floating point (floats and doubles).  In particular, the GHC designers decided to perform their own operations for encoding a GMP number (really, the array of mp_limb_t) to floating point and doubles.  GMP provides functions for some of these operations: 
     25 
     26|| '''Operation''' || '''GMP function''' || '''Notes''' || 
     27|| convert to double || mpz_get_d() || rounds toward zero, avoids overflow || 
     28|| convert to double[[BR]]return exponent separately || mpz_get_d_2exp() || rounds toward zero, avoids overflow || 
     29|| convert to float || mpz_get_f() 
     30|| convert to string[[BR]]supports bases from 2 to 36 || mpz_get_str() || more efficient than Haskell version,[[BR]]with !PackedString, no conversion to Haskell string necessary || 
     31|| conversion from string || mpz_inp_str() || more efficient than Haskell version,[[BR]]with !PackedString, input string efficiently convertable to a c_str[[BR]](it is already an array of chars) || 
     32 
     33The real problem (no pun intended) is conversion to doubles and floats.  The `__encodeDouble` function, in [[GhcFile(rts/StgPrimFloat.c)]], may overflow: 
     34{{{ 
     35StgDouble 
     36__encodeDouble (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e */ 
     37{ 
     38    StgDouble r; 
     39    const mp_limb_t *const arr = (const mp_limb_t *)ba; 
     40    I_ i; 
     41 
     42    /* Convert MP_INT to a double; knows a lot about internal rep! */ 
     43    for(r = 0.0, i = __abs(size)-1; i >= 0; i--) 
     44        r = (r * GMP_BASE) + arr[i]; 
     45 
     46    /* Now raise to the exponent */ 
     47    if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */ 
     48        r = ldexp(r, e); 
     49 
     50    /* sign is encoded in the size */ 
     51    if (size < 0) 
     52        r = -r; 
     53 
     54    return r; 
     55} 
     56}}} 
     57If the problem section hurt your eyes, here is some more pain to drive the point home: 
     58{{{ 
     59    /* __abs(size) is the size of the array --PDT */ 
     60    for(r = 0.0, i = __abs(size)-1; i >= 0; i--) 
     61        r = (r * GMP_BASE) + arr[i]; 
     62}}}  
     63Note the problems: 
     64 * GMP numbers may be very large, much greater than the maximum value for a float (addition and multiplication operations may overflow);  
     65 * there is no specification of rounding mode--this operation relies on whatever hardware rounding may take effect; and, 
     66 * certain floating point hardware exceptions (traps) for overflow or underflow may be triggered (though this is not generally the case; note: on x86 machines floating point exceptions may not be triggered but the resulting float may denormalized). 
     67Nowhere in this function is there a check for whether the GMP number (the `const mp_limb_t *arr`) may be greater than the size of a double and either warn the user (possibly with an `ArithException` or round the resulting double toward zero.  Compare the problem section in `__encodeDouble` to the exponent check in the internal GMP function, `mpn_get_d`, that avoids overflow by comparing the relative magnitude of the GMP number with the maximum magnitude of the double: 
     68{{{ 
     69   /* note: exp is a long int, given as an argument to min_get_d as 0L --PDT */ 
     70  /* Adjust exp to a radix point just above {ptr,size}, guarding against 
     71     overflow.  After this exp can of course be reduced to anywhere within 
     72     the {ptr,size} region without underflow.  */ 
     73  if (UNLIKELY ((unsigned long) (GMP_NUMB_BITS * size) 
     74                > (unsigned long) (LONG_MAX - exp))) 
     75    { 
     76      if (_GMP_IEEE_FLOATS) 
     77        goto ieee_infinity; 
     78 
     79      /* generic */ 
     80      exp = LONG_MAX; 
     81    } 
     82  else 
     83    { 
     84      exp += GMP_NUMB_BITS * size; 
     85    } 
     86}}} 
     87(If you want to see the full `mpn_get_d` function, it is in the file `[toplevel gmp source]/mpn/generic/get_d.c` .)  A replacement library for GMP might use the GMP strategies of including a special bitwise conversion (with appropriate masks) and a hardware-based version.  An unconventional solution might perform the rounding manually (but with relatively portable predictability) using interval arithmetic. 
     88 
     89As for printing GMP numbers, [[GhcFile(libraries/base/GHC/Num.lhs)]] defines a special internal function `jtos` to handle `showsPrec` and `showList` for `Integer` as an instance of the class `Show`.  The `jtos` function uses the primitive GMP-based function `quotRemInteger` and performs many conversions from Integer to Int.  This is not as efficient as the internal GMP functions, especially for large numbers, because each conversion allocates extra storage for GMP as noted in [wiki:ReplacingGMPNotes#OptimisationOpportunities Replacing GMP -- Optimisation Opportunities].