#include <cmath>
#include <iostream>

// In a previous exercise, you were asked to calculate
// integer powers of a number where the base could be
// either an integer or a real number.
//
// Write two overloaded functions that implement these
// calculations. The compiler will determine which
// function to call from the types of the arguments.
//
// If the base is an integer and the exponent is negative,
// the result generally cannot be represented by an
// integer. For example,
//
//       -3
//      2   = 1/8.
//
// There are two exceptions:
//
//       m
//      1  = 1
//
// and
//
//        m
//    (-1)  = 1   if m is even,
//           = -1  if m is odd.
//
// For all other integer bases with negative exponents,
// return 0.
//
// Finally, compare our calculation with std::pow() from
// the C++ standard library.

// Function declarations
int main();
int my_pow( int n, int m );
double my_pow( double x, int m );

// Function definitions
int main() {
    std::cout << my_pow( 3, 5 ) << std::endl;
    std::cout << my_pow( 1.0032, 105 ) << std::endl;

    std::cout << my_pow( 3, 0 ) << std::endl;
    std::cout << my_pow( 1.0032, 0 ) << std::endl;

    std::cout << my_pow( 3, -5 ) << std::endl;
    std::cout << my_pow( 1.0032, -105 ) << std::endl;

    std::cout << my_pow( -1, -1 ) << std::endl;
    std::cout << my_pow( -1, -13 ) << std::endl;
    std::cout << my_pow( -1, -14 ) << std::endl;
    std::cout << my_pow(  1, -13 ) << std::endl;
    std::cout << my_pow(  1, -14 ) << std::endl;

    // Specify that std::cout should print floating-point
    // numbers with sixteen digits of precision.
    std::cout.precision( 16 );

    // Compare our calculation with std::pow().
    //
    // Are either of these answers exactly equal to
    //
    //             1530513531
    //      (1.0000000001)             ?
    //
    // Why might they differ?
    std::cout << my_pow( 1.0000000001, 1530513531 )
              << std::endl;
    std::cout << std::pow( 1.0000000001, 1530513531 )
              << std::endl;

    return 0;
}

int my_pow( int n, int m ) {
    if ( m < 0 ) {
        if ( n == 1 ) {
            return 1;
        } else if ( n == -1 ) {
            if ( m%2 == 0 ) {
                return 1;
            } else {
                return -1;
            }
        } else {
            return 0;
        }
    }

    int result{ 1 };

    for ( int k{ 1 }; k <= m; ++k ) {
        result *= n;
    }

    return result;
}

double my_pow( double x, int m ) {
    // If m < 0, then
    //
    //       m    / 1 \ -m
    //      x  = | --- |   .
    //            \ x /
    //
    // Thus, replace x by its reciprocal and negate m.
    //
    // For simplicity, we assume that m is not the
    // smallest possible int, because negating that
    // value cannot be represented by an int.
    if ( m < 0 ) {
        x = 1.0/x;
        m = -m;
    }

    double result{ 1.0 };

    for ( int k{ 1 }; k <= m; ++k ) {
        result *= x;
    }

    return result;
}

// Answer to the question:
//
// Neither calculation gives the exact mathematical value
//
//                 1530513531
//      (1.0000000001)             .
//
// Using the exact decimal value 1.0000000001, the result
// is approximately
//
//      1.1653848235205852168631164303491760776145...
//
// Our repeated-multiplication calculation and std::pow()
// produce slightly different results.
//
// Why?
//
// The first source of error occurs before either power
// calculation begins. The decimal number
//
//      1.0000000001
//
// cannot be represented exactly as a binary
// floating-point number of type double. The stored value
// is therefore slightly different from the exact decimal
// value.
//
// Raising this slightly different value to the enormous
// power 1530513531 magnifies that initial difference.
//
// Additional rounding occurs during the calculations.
// Our function performs more than 1.5 billion
// multiplications, with each multiplication rounded to a
// representable double. The implementation of std::pow()
// uses a different algorithm and therefore accumulates
// rounding errors differently.
//
// To find the value corresponding to the exact decimal
// number, we can use the binomial theorem. Let
//
//      e = 0.0000000001.
//
// Then
//
//                 1530513531
//      (1 + e)
//
// can be expanded as
//
//                         / 1530513531 \  2
//      1 + 1530513531*e + |            | e
//                         \      2     /
//
//          / 1530513531 \  3
//        + |            | e
//          \      3     /
//
//          / 1530513531 \  4
//        + |            | e  + ...
//          \      4     /
//
// where the quantities in parentheses are binomial
// coefficients.
//
// The first terms are approximately
//
//      1
//      0.1530513531
//      0.01171235833521787215
//      0.000597530762951562349605577245
//      0.0000228632229023381846201700359430367590
//      0.0000006998494384568958297958719410178792...
//      0.0000000178521505286963939245929304115290...
//      0.0000000003903279690643675931473901444890...
//      0.0000000000074675279431058528477426292690...
//      0.0000000000001269905833367098044879794902...
//      0.0000000000000019436080496350223475882890...
//      0.0000000000000000270428945408258256192822...
//    + 0.0000000000000000003449126309555676506289...
//      -----------------------------------------------
//      1.1653848235205852168631164303491760776145...
//
// Consequently, this is an important lesson about
// floating-point arithmetic: a decimal value may already
// have been approximated when it is stored as a double,
// and subsequent calculations can magnify that initial
// approximation. Different algorithms may also produce
// slightly different floating-point results.