DevWiki/calculator
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
ef3f5e9cbb
calculator/src/CalcManager/Ratpack/exp.cpp
Will 1113ff4b86 Updating comments per the C++ core guidelines and removing trailing whitespace (#194) 
Fixed comments that were inconsistent with the style guidelines described in C++ core guidelines and the modern C++/WinRT language projections and removed trailing whitespace.

Inserted a space after the beginning of the comment so the text wasn't touching the // on all occurrences.

Removed all occurrences of trailing whitespace
2019-03-14 23:30:07 -07:00

541 lines
16 KiB
C++
Raw Blame History

// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
//-----------------------------------------------------------------------------
// Package Title ratpak
// File exp.c
// Copyright (C) 1995-96 Microsoft
// Date 01-16-95
//
//
// Description
//
// Contains exp, and log functions for rationals
//
//
//-----------------------------------------------------------------------------
#include "pch.h"
#include "ratpak.h"
//-----------------------------------------------------------------------------
//
// FUNCTION: exprat
//
// ARGUMENTS: x PRAT representation of number to exponentiate
//
// RETURN: exp of x in PRAT form.
//
// EXPLANATION: This uses Taylor series
//
// n
// ___
// \ ] X
// \ thisterm ; where thisterm = thisterm * ---------
// / j j+1 j j+1
// /__]
// j=0
//
// thisterm = X ; and stop when thisterm < precision used.
// 0 n
//
//-----------------------------------------------------------------------------
void _exprat( PRAT *px, int32_t precision)
{
CREATETAYLOR();
addnum(&(pret->pp),num_one, BASEX);
addnum(&(pret->pq),num_one, BASEX);
DUPRAT(thisterm,pret);
n2=longtonum(0L, BASEX);
do {
NEXTTERM(*px, INC(n2) DIVNUM(n2), precision);
} while ( !SMALL_ENOUGH_RAT( thisterm, precision) );
DESTROYTAYLOR();
}
void exprat( PRAT *px, uint32_t radix, int32_t precision)
{
PRAT pwr= nullptr;
PRAT pint= nullptr;
long intpwr;
if ( rat_gt( *px, rat_max_exp, precision) || rat_lt( *px, rat_min_exp, precision) )
{
// Don't attempt exp of anything large.
throw( CALC_E_DOMAIN );
}
DUPRAT(pwr,rat_exp);
DUPRAT(pint,*px);
intrat(&pint, radix, precision);
intpwr = rattolong(pint, radix, precision);
ratpowlong( &pwr, intpwr, precision);
subrat(px, pint, precision);
// It just so happens to be an integral power of e.
if ( rat_gt( *px, rat_negsmallest, precision) && rat_lt( *px, rat_smallest, precision) )
{
DUPRAT(*px,pwr);
}
else
{
_exprat(px, precision);
mulrat(px, pwr, precision);
}
destroyrat( pwr );
destroyrat( pint );
}
//-----------------------------------------------------------------------------
//
// FUNCTION: lograt, _lograt
//
// ARGUMENTS: x PRAT representation of number to logarithim
//
// RETURN: log of x in PRAT form.
//
// EXPLANATION: This uses Taylor series
//
// n
// ___
// \ ] j*(1-X)
// \ thisterm ; where thisterm = thisterm * ---------
// / j j+1 j j+1
// /__]
// j=0
//
// thisterm = X ; and stop when thisterm < precision used.
// 0 n
//
// Number is scaled between one and e_to_one_half prior to taking the
// log. This is to keep execution time from exploding.
//
//
//-----------------------------------------------------------------------------
void _lograt( PRAT *px, int32_t precision)
{
CREATETAYLOR();
createrat(thisterm);
// sub one from x
(*px)->pq->sign *= -1;
addnum(&((*px)->pp),(*px)->pq, BASEX);
(*px)->pq->sign *= -1;
DUPRAT(pret,*px);
DUPRAT(thisterm,*px);
n2=longtonum(1L, BASEX);
(*px)->pp->sign *= -1;
do {
NEXTTERM(*px, MULNUM(n2) INC(n2) DIVNUM(n2), precision);
TRIMTOP(*px, precision);
} while ( !SMALL_ENOUGH_RAT( thisterm, precision) );
DESTROYTAYLOR();
}
void lograt( PRAT *px, int32_t precision)
{
bool fneglog;
PRAT pwr = nullptr; // pwr is the large scaling factor.
PRAT offset = nullptr; // offset is the incremental scaling factor.
// Check for someone taking the log of zero or a negative number.
if ( rat_le( *px, rat_zero, precision) )
{
throw( CALC_E_DOMAIN );
}
// Get number > 1, for scaling
fneglog = rat_lt( *px, rat_one, precision);
if ( fneglog )
{
// WARNING: This is equivalent to doing *px = 1 / *px
PNUMBER pnumtemp= nullptr;
pnumtemp = (*px)->pp;
(*px)->pp = (*px)->pq;
(*px)->pq = pnumtemp;
}
// Scale the number within BASEX factor of 1, for the large scale.
// log(x*2^(BASEXPWR*k)) = BASEXPWR*k*log(2)+log(x)
if ( LOGRAT2(*px) > 1 )
{
// Take advantage of px's base BASEX to scale quickly down to
// a reasonable range.
long intpwr;
intpwr=LOGRAT2(*px)-1;
(*px)->pq->exp += intpwr;
pwr=longtorat(intpwr*BASEXPWR);
mulrat(&pwr, ln_two, precision);
// ln(x+e)-ln(x) looks close to e when x is close to one using some
// expansions. This means we can trim past precision digits+1.
TRIMTOP(*px, precision);
}
else
{
DUPRAT(pwr,rat_zero);
}
DUPRAT(offset,rat_zero);
// Scale the number between 1 and e_to_one_half, for the small scale.
while ( rat_gt( *px, e_to_one_half, precision) )
{
divrat( px, e_to_one_half, precision);
addrat( &offset, rat_one, precision);
}
_lograt(px, precision);
// Add the large and small scaling factors, take into account
// small scaling was done in e_to_one_half chunks.
divrat(&offset, rat_two, precision);
addrat(&pwr, offset, precision);
// And add the resulting scaling factor to the answer.
addrat(px, pwr, precision);
trimit(px, precision);
// If number started out < 1 rescale answer to negative.
if ( fneglog )
{
(*px)->pp->sign *= -1;
}
destroyrat(offset);
destroyrat(pwr);
}
void log10rat( PRAT *px, int32_t precision)
{
lograt(px, precision);
divrat(px, ln_ten, precision);
}
//
// return if the given x is even number. The assumption here is its denominator is 1 and we are testing the numerator is
// even or not
bool IsEven(PRAT x, uint32_t radix, int32_t precision)
{
PRAT tmp = nullptr;
bool bRet = false;
DUPRAT(tmp, x);
divrat(&tmp, rat_two, precision);
fracrat(&tmp, radix, precision);
addrat(&tmp, tmp, precision);
subrat(&tmp, rat_one, precision);
if ( rat_lt( tmp, rat_zero, precision))
{
bRet = true;
}
destroyrat(tmp);
return bRet;
}
//---------------------------------------------------------------------------
//
// FUNCTION: powrat
//
// ARGUMENTS: PRAT *px, PRAT y, uint32_t radix, int32_t precision
//
// RETURN: none, sets *px to *px to the y.
//
// EXPLANATION: Calculates the power of both px and
// handles special cases where px is a perfect root.
// Assumes, all checking has been done on validity of numbers.
//
//
//---------------------------------------------------------------------------
void powrat(PRAT *px, PRAT y, uint32_t radix, int32_t precision)
{
// Handle cases where px or y is 0 by calling powratcomp directly
if (zerrat(*px) || zerrat(y))
{
powratcomp(px, y, radix, precision);
return;
}
// When y is 1, return px
if (rat_equ(y, rat_one, precision))
{
return;
}
try
{
powratNumeratorDenominator(px, y, radix, precision);
}
catch (...)
{
// If calculating the power using numerator/denominator
// failed, fall back to the less accurate method of
// passing in the original y
powratcomp(px, y, radix, precision);
}
}
void powratNumeratorDenominator(PRAT *px, PRAT y, uint32_t radix, int32_t precision)
{
// Prepare rationals
PRAT yNumerator = nullptr;
PRAT yDenominator = nullptr;
DUPRAT(yNumerator, rat_zero); // yNumerator->pq is 1 one
DUPRAT(yDenominator, rat_zero); // yDenominator->pq is 1 one
DUPNUM(yNumerator->pp, y->pp);
DUPNUM(yDenominator->pp, y->pq);
// Calculate the following use the Powers of Powers rule:
// px ^ (yNum/yDenom) == px ^ yNum ^ (1/yDenom)
// 1. For px ^ yNum, we call powratcomp directly which will call ratpowlong
// and store the result in pxPowNum
// 2. For pxPowNum ^ (1/yDenom), we call powratcomp
// 3. Validate the result of 2 by adding/subtracting 0.5, flooring and call powratcomp with yDenom
// on the floored result.
// 1. Initialize result.
PRAT pxPow = nullptr;
DUPRAT(pxPow, *px);
// 2. Calculate pxPow = px ^ yNumerator
// if yNumerator is not 1
if (!rat_equ(yNumerator, rat_one, precision))
{
powratcomp(&pxPow, yNumerator, radix, precision);
}
// 2. Calculate pxPowNumDenom = pxPowNum ^ (1/yDenominator),
// if yDenominator is not 1
if (!rat_equ(yDenominator, rat_one, precision))
{
// Calculate 1 over y
PRAT oneoveryDenom = nullptr;
DUPRAT(oneoveryDenom, rat_one);
divrat(&oneoveryDenom, yDenominator, precision);
// ##################################
// Take the oneoveryDenom power
// ##################################
PRAT originalResult = nullptr;
DUPRAT(originalResult, pxPow);
powratcomp(&originalResult, oneoveryDenom, radix, precision);
// ##################################
// Round the originalResult to roundedResult
// ##################################
PRAT roundedResult = nullptr;
DUPRAT(roundedResult, originalResult);
if (roundedResult->pp->sign == -1)
{
subrat(&roundedResult, rat_half, precision);
}
else
{
addrat(&roundedResult, rat_half, precision);
}
intrat(&roundedResult, radix, precision);
// ##################################
// Take the yDenom power of the roundedResult.
// ##################################
PRAT roundedPower = nullptr;
DUPRAT(roundedPower, roundedResult);
powratcomp(&roundedPower, yDenominator, radix, precision);
// ##################################
// if roundedPower == px,
// we found an exact power in roundedResult
// ##################################
if (rat_equ(roundedPower, pxPow, precision))
{
DUPRAT(*px, roundedResult);
}
else
{
DUPRAT(*px, originalResult);
}
destroyrat(oneoveryDenom);
destroyrat(originalResult);
destroyrat(roundedResult);
destroyrat(roundedPower);
}
else
{
DUPRAT(*px, pxPow);
}
destroyrat(yNumerator);
destroyrat(yDenominator);
destroyrat(pxPow);
}
//---------------------------------------------------------------------------
//
// FUNCTION: powratcomp
//
// ARGUMENTS: PRAT *px, and PRAT y
//
// RETURN: none, sets *px to *px to the y.
//
// EXPLANATION: This uses x^y=e(y*ln(x)), or a more exact calculation where
// y is an integer.
// Assumes, all checking has been done on validity of numbers.
//
//
//---------------------------------------------------------------------------
void powratcomp(PRAT *px, PRAT y, uint32_t radix, int32_t precision)
{
long sign = ((*px)->pp->sign * (*px)->pq->sign);
// Take the absolute value
(*px)->pp->sign = 1;
(*px)->pq->sign = 1;
if ( zerrat( *px ) )
{
// *px is zero.
if ( rat_lt( y, rat_zero, precision) )
{
throw( CALC_E_DOMAIN );
}
else if ( zerrat( y ) )
{
// *px and y are both zero, special case a 1 return.
DUPRAT(*px,rat_one);
// Ensure sign is positive.
sign = 1;
}
}
else
{
PRAT pxint= nullptr;
DUPRAT(pxint,*px);
subrat(&pxint, rat_one, precision);
if ( rat_gt( pxint, rat_negsmallest, precision) &&
rat_lt( pxint, rat_smallest, precision) && ( sign == 1 ) )
{
// *px is one, special case a 1 return.
DUPRAT(*px,rat_one);
// Ensure sign is positive.
sign = 1;
}
else
{
// Only do the exp if the number isn't zero or one
PRAT podd = nullptr;
DUPRAT(podd,y);
fracrat(&podd, radix, precision);
if ( rat_gt( podd, rat_negsmallest, precision) && rat_lt( podd, rat_smallest, precision) )
{
// If power is an integer let ratpowlong deal with it.
PRAT iy = nullptr;
long inty;
DUPRAT(iy,y);
subrat(&iy, podd, precision);
inty = rattolong(iy, radix, precision);
PRAT plnx = nullptr;
DUPRAT(plnx,*px);
lograt(&plnx, precision);
mulrat(&plnx, iy, precision);
if ( rat_gt( plnx, rat_max_exp, precision) || rat_lt( plnx, rat_min_exp, precision) )
{
// Don't attempt exp of anything large or small.A
destroyrat(plnx);
destroyrat(iy);
destroyrat(pxint);
destroyrat(podd);
throw( CALC_E_DOMAIN );
}
destroyrat(plnx);
ratpowlong(px, inty, precision);
if ( ( inty & 1 ) == 0 )
{
sign=1;
}
destroyrat(iy);
}
else
{
// power is a fraction
if ( sign == -1 )
{
// Need to throw an error if the exponent has an even denominator.
// As a first step, the numerator and denominator must be divided by 2 as many times as
// possible, so that 2/6 is allowed.
// If the final numerator is still even, the end result should be positive.
PRAT pNumerator = nullptr;
PRAT pDenominator = nullptr;
bool fBadExponent = false;
// Get the numbers in arbitrary precision rational number format
DUPRAT(pNumerator, rat_zero); // pNumerator->pq is 1 one
DUPRAT(pDenominator, rat_zero); // pDenominator->pq is 1 one
DUPNUM(pNumerator->pp, y->pp);
pNumerator->pp->sign = 1;
DUPNUM(pDenominator->pp, y->pq);
pDenominator->pp->sign = 1;
while (IsEven(pNumerator, radix, precision) && IsEven(pDenominator, radix, precision)) // both Numerator & denominator is even
{
divrat(&pNumerator, rat_two, precision);
divrat(&pDenominator, rat_two, precision);
}
if (IsEven(pDenominator, radix, precision)) // denominator is still even
{
fBadExponent = true;
}
if (IsEven(pNumerator, radix, precision)) // numerator is still even
{
sign = 1;
}
destroyrat(pNumerator);
destroyrat(pDenominator);
if (fBadExponent)
{
throw( CALC_E_DOMAIN );
}
}
else
{
// If the exponent is not odd disregard the sign.
sign = 1;
}
lograt(px, precision);
mulrat(px, y, precision);
exprat(px, radix, precision);
}
destroyrat(podd);
}
destroyrat(pxint);
}
(*px)->pp->sign *= sign;
}
Reference in New Issue View Git Blame Copy Permalink