// 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 numberator 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, fallback 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; }