fork download 
  1. #include <assert.h>         // assert
  2. #include <float.h>          // FLT_RADIX, DBL_MANT_DIG
  3. #include <iomanip>
  4. #include <iostream>
  5. #include <limits.h>         // CHAR_BIT
  6. #include <math.h>           // floor
  7. using namespace std;
  8.  
  9. #define STATIC_ASSERT( e )   \
  10.     typedef char StaticAssert[(e)? 1 : -1];
  11.  
  12. int const   bitsPerByte     = CHAR_BIT;
  13. int const   logHighestIntN  = sizeof( int )*bitsPerByte - 2;
  14. int const   highestIntN     = 1 << logHighestIntN;
  15.  
  16. bool isInteger( double const x )    { return (x - floor( x ) == 0); }
  17. int pow4( int const x )             { return x*x*x*x; }
  18. int pow8( int const x )             { return pow4( x*x ); }
  19.  
  20. double root( double const n, double x )
  21. {
  22.     assert( n > 0 );
  23.     assert( isInteger( n ) );
  24.  
  25.     return pow( x, 1.0/n );       // TODO
  26. }
  27.  
  28. // Returns x^(2^logN)
  29. double expPower( double const x, int const logN )
  30. {
  31.     assert( logN >= 0 );
  32.  
  33.     double  result  = x;
  34.     for( int i = 1;  i <= logN;  ++i )
  35.     {
  36.         result *= result;
  37.     }
  38.     return result;
  39. }
  40.  
  41. #if 0
  42. double intPower( double const x, int const n )
  43. {
  44.     assert( n >= 0 );
  45.  
  46.     if( n == 0 )
  47.     {
  48.         return 1;
  49.     }
  50.     else if( n % 2 == 1 )
  51.     {
  52.         return x*intPower( x, n - 1 );
  53.     }
  54.     else
  55.     {
  56.         return intPower( x*x, n/2 );
  57.     }
  58. }
  59. #else
  60. double intPower( double const x, int const n )
  61. {
  62.     assert( n >= 0 );
  63.  
  64.     int bitmask     = 1;
  65.  
  66.     for( int i = n;  i != 0;  i >>= 1 )
  67.     {
  68.         bitmask <<= 1;
  69.     }
  70.     // bitmask now has a 1 bit in the most significant bit pos of n.
  71.  
  72.     double  result  = 1.0;
  73.     for( ;; )
  74.     {
  75.         int const bit = (n & bitmask);
  76.         if( bit != 0 ) { result *= x; }
  77.         bitmask = (bitmask >> 1);
  78.         if( bitmask == 0 ) { break; }
  79.         result *= result;
  80.     }
  81.  
  82.     return result;
  83. }
  84. #endif
  85.  
  86. double fractionalPower( double const x, double const n )
  87. {
  88.     wclog << "fractionalPower(" << x << ", " << n << ")" << endl;
  89.     assert( 0 < n && n < 1 );
  90.     STATIC_ASSERT( FLT_RADIX == 2 );
  91.  
  92.     int const radix     = pow4( FLT_RADIX );
  93.     int const nDigits   = DBL_MANT_DIG/4;
  94.  
  95.     double  result          = 1.0;
  96.     double  digitWeight     = 1.0;
  97.     double  digits          = n;
  98.     for( int i = 1;  i <= nDigits && digits != 0;  ++i )
  99.     {
  100.         double const    leftShiftedDigits   =
  101.             radix*digits;
  102.         int const       digit               =
  103.             static_cast<int>( floor( leftShiftedDigits ) );
  104.  
  105.         digitWeight *= radix;
  106.  
  107.         double const    multiplier          =
  108.             root( digitWeight, intPower( x, digit ) );
  109.  
  110.         result *= multiplier;
  111.         wclog << "!" << setprecision( 17 ) << result << endl;
  112.         digits = leftShiftedDigits - digit;
  113.     }
  114.     return result;
  115. }
  116.  
  117. double integralPower( double const x, double const n )
  118. {
  119.     wclog << "integralPower(" << x << ", " << n << ")" << endl;
  120.     assert( n >= 0 );
  121.     assert( isInteger( n ) );
  122.  
  123.     if( n > highestIntN )
  124.     {
  125.         double const    xToHiN      = expPower( x, logHighestIntN );
  126.         double const    mAsDouble   = floor( n/highestIntN );
  127.         int const       m           = static_cast<int>( mAsDouble );
  128.         int const       r           = static_cast<int>( n - mAsDouble*highestIntN );
  129.  
  130.         return intPower( xToHiN, m )*intPower( x, r );
  131.     }
  132.     else
  133.     {
  134.         return intPower( x, static_cast<int>( n ) );
  135.     }
  136. }
  137.  
  138. double positivePower( double const x, double const n )
  139. {
  140.     assert( n > 0 );
  141.     double const    intPart         = floor( n );
  142.     double const    fractionalPart  = n - intPart;
  143.  
  144.     if( fractionalPart == 0 )
  145.     {
  146.         return integralPower( x, intPart );
  147.     }
  148.     else
  149.     {
  150.         return integralPower( x, intPart )*fractionalPower( x, fractionalPart );
  151.     }
  152. }
  153.  
  154. double power( double const x, double const n )
  155. {
  156.     if( n == 0 )
  157.     {
  158.         return 1;
  159.     }
  160.     else if( n < 0 )
  161.     {
  162.         return 1.0/positivePower( x, -n );
  163.     }
  164.     else
  165.     {
  166.         return positivePower( x, n );
  167.     }
  168. }
  169.  
  170. int main()
  171. {
  172.     double const x = 2.7; double const n = 3.14;
  173.  
  174.     wcout << setprecision( 16 );
  175.     wcout << "Stdlib function result: " << pow( x, n ) << endl;
  176.     wcout << "Custom function result: " << power( x, n ) << endl;
  177. }
  178.  
Success #stdin #stdout 0.01s 2896KB
comments (?)
stdin 
Standard input is empty
stdout 
Stdlib function result: 22.61945931074701
Custom function result: 22.61945931074701
https://ideone.com/Erdee
language:
C++ (gcc 8.3)
created:
14 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!