Hello GitHub

This commit is contained in:
Howard Wolosky
2019-01-28 16:24:37 -08:00
parent 456fe5e355
commit c13b8a099e
822 changed files with 276650 additions and 75 deletions

View File

@@ -0,0 +1,82 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
// CalcErr.h
//
// Defines the error codes thrown by ratpak and caught by Calculator
//
//
// Ratpak errors are 32 bit values layed out as follows:
//
// 3 3 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1
// 1 0 9 8 7 6 5 4 3 2 1 0 9 8 7 6 5 4 3 2 1 0 9 8 7 6 5 4 3 2 1 0
// +-+-------+---------------------+-------------------------------+
// |S| R | Facility | Code |
// +-+-------+---------------------+-------------------------------+
//
// where
//
// S - Severity - indicates success/fail
//
// 0 - Success
// 1 - Fail
//
// R - Reserved - not currently used for anything
//
// r - reserved portion of the facility code. Reserved for internal
// use. Used to indicate HRESULT values that are not status
// values, but are instead message ids for display strings.
//
// Facility - is the facility code
//
// Code - is the actual error code
//
// This format is based losely on an OLE HRESULT and is compatible with the
// SUCCEEDED and FAILED marcos as well as the HRESULT_CODE macro
// CALC_E_DIVIDEBYZERO
//
// The current operation would require a divide by zero to complete
#define CALC_E_DIVIDEBYZERO ((DWORD)0x80000000)
// CALC_E_DOMAIN
//
// The given input is not within the domain of this function
#define CALC_E_DOMAIN ((DWORD)0x80000001)
// CALC_E_INDEFINITE
//
// The result of this function is undefined
#define CALC_E_INDEFINITE ((DWORD)0x80000002)
// CALC_E_POSINFINITY
//
// The result of this function is Positive Infinity.
#define CALC_E_POSINFINITY ((DWORD)0x80000003)
// CALC_E_NEGINFINITY
//
// The result of this function is Negative Infinity
#define CALC_E_NEGINFINITY ((DWORD)0x80000004)
// CALC_E_INVALIDRANGE
//
// The given input is within the domain of the function but is beyond
// the range for which calc can successfully compute the answer
#define CALC_E_INVALIDRANGE ((DWORD)0x80000006)
// CALC_E_OUTOFMEMORY
//
// There is not enough free memory to complete the requested function
#define CALC_E_OUTOFMEMORY ((DWORD)0x80000007)
// CALC_E_OVERFLOW
//
// The result of this operation is an overflow
#define CALC_E_OVERFLOW ((DWORD)0x80000008)
// CALC_E_NORESULT
//
// The result of this operation is undefined
#define CALC_E_NORESULT ((DWORD)0x80000009)

View File

@@ -0,0 +1,367 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
//-----------------------------------------------------------------------------
// Package Title ratpak
// File basex.c
// Copyright (C) 1995-97 Microsoft
// Date 03-14-97
//
//
// Description
//
// Contains number routines for internal base computations, these assume
// internal base is a power of 2.
//
//-----------------------------------------------------------------------------
#include "pch.h"
#include "ratpak.h"
void _mulnumx( PNUMBER *pa, PNUMBER b );
//----------------------------------------------------------------------------
//
// FUNCTION: mulnumx
//
// ARGUMENTS: pointer to a number and a second number, the
// base is always BASEX.
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the number equivalent of *pa *= b.
// This is a stub which prevents multiplication by 1, this is a big speed
// improvement.
//
//----------------------------------------------------------------------------
void __inline mulnumx( PNUMBER *pa, PNUMBER b )
{
if ( b->cdigit > 1 || b->mant[0] != 1 || b->exp != 0 )
{
// If b is not one we multiply
if ( (*pa)->cdigit > 1 || (*pa)->mant[0] != 1 || (*pa)->exp != 0 )
{
// pa and b are both nonone.
_mulnumx( pa, b );
}
else
{
// if pa is one and b isn't just copy b. and adjust the sign.
long sign = (*pa)->sign;
DUPNUM(*pa,b);
(*pa)->sign *= sign;
}
}
else
{
// B is +/- 1, But we do have to set the sign.
(*pa)->sign *= b->sign;
}
}
//----------------------------------------------------------------------------
//
// FUNCTION: _mulnumx
//
// ARGUMENTS: pointer to a number and a second number, the
// base is always BASEX.
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the number equivalent of *pa *= b.
// Assumes the base is BASEX of both numbers. This algorithm is the
// same one you learned in gradeschool, except the base isn't 10 it's
// BASEX.
//
//----------------------------------------------------------------------------
void _mulnumx( PNUMBER *pa, PNUMBER b )
{
PNUMBER c= nullptr; // c will contain the result.
PNUMBER a= nullptr; // a is the dereferenced number pointer from *pa
MANTTYPE *ptra; // ptra is a pointer to the mantissa of a.
MANTTYPE *ptrb; // ptrb is a pointer to the mantissa of b.
MANTTYPE *ptrc; // ptrc is a pointer to the mantissa of c.
MANTTYPE *ptrcoffset; // ptrcoffset, is the anchor location of the next
// single digit multiply partial result.
long iadigit=0; // Index of digit being used in the first number.
long ibdigit=0; // Index of digit being used in the second number.
MANTTYPE da=0; // da is the digit from the fist number.
TWO_MANTTYPE cy=0; // cy is the carry resulting from the addition of
// a multiplied row into the result.
TWO_MANTTYPE mcy=0; // mcy is the resultant from a single
// multiply, AND the carry of that multiply.
long icdigit=0; // Index of digit being calculated in final result.
a=*pa;
ibdigit = a->cdigit + b->cdigit - 1;
createnum( c, ibdigit + 1 );
c->cdigit = ibdigit;
c->sign = a->sign * b->sign;
c->exp = a->exp + b->exp;
ptra = a->mant;
ptrcoffset = c->mant;
for ( iadigit = a->cdigit; iadigit > 0; iadigit-- )
{
da = *ptra++;
ptrb = b->mant;
// Shift ptrc, and ptrcoffset, one for each digit
ptrc = ptrcoffset++;
for ( ibdigit = b->cdigit; ibdigit > 0; ibdigit-- )
{
cy = 0;
mcy = (DWORDLONG)da * (*ptrb);
if ( mcy )
{
icdigit = 0;
if ( ibdigit == 1 && iadigit == 1 )
{
c->cdigit++;
}
}
// If result is nonzero, or while result of carry is nonzero...
while ( mcy || cy )
{
// update carry from addition(s) and multiply.
cy += (TWO_MANTTYPE)ptrc[icdigit]+((DWORD)mcy&((DWORD)~BASEX));
// update result digit from
ptrc[icdigit++]=(MANTTYPE)((DWORD)cy&((DWORD)~BASEX));
// update carries from
mcy >>= BASEXPWR;
cy >>= BASEXPWR;
}
ptrb++;
ptrc++;
}
}
// prevent different kinds of zeros, by stripping leading duplicate zeroes.
// digits are in order of increasing significance.
while ( c->cdigit > 1 && c->mant[c->cdigit-1] == 0 )
{
c->cdigit--;
}
destroynum( *pa );
*pa=c;
}
//-----------------------------------------------------------------------------
//
// FUNCTION: numpowlongx
//
// ARGUMENTS: root as number power as long
// number.
//
// RETURN: None root is changed.
//
// DESCRIPTION: changes numeric representation of root to
// root ** power. Assumes base BASEX
// decomposes the exponent into it's sums of powers of 2, so on average
// it will take n+n/2 multiplies where n is the highest on bit.
//
//-----------------------------------------------------------------------------
void numpowlongx( _Inout_ PNUMBER *proot, _In_ long power )
{
PNUMBER lret = longtonum( 1, BASEX );
// Once the power remaining is zero we are done.
while ( power > 0 )
{
// If this bit in the power decomposition is on, multiply the result
// by the root number.
if ( power & 1 )
{
mulnumx( &lret, *proot );
}
// multiply the root number by itself to scale for the next bit (i.e.
// square it.
mulnumx( proot, *proot );
// move the next bit of the power into place.
power >>= 1;
}
destroynum( *proot );
*proot=lret;
}
void _divnumx( PNUMBER *pa, PNUMBER b, int32_t precision);
//----------------------------------------------------------------------------
//
// FUNCTION: divnumx
//
// ARGUMENTS: pointer to a number, a second number and precision.
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the number equivalent of *pa /= b.
// Assumes radix is the internal radix representation.
// This is a stub which prevents division by 1, this is a big speed
// improvement.
//
//----------------------------------------------------------------------------
void __inline divnumx( PNUMBER *pa, PNUMBER b, int32_t precision)
{
if ( b->cdigit > 1 || b->mant[0] != 1 || b->exp != 0 )
{
// b is not one.
if ( (*pa)->cdigit > 1 || (*pa)->mant[0] != 1 || (*pa)->exp != 0 )
{
// pa and b are both not one.
_divnumx( pa, b, precision);
}
else
{
// if pa is one and b is not one, just copy b, and adjust the sign.
long sign = (*pa)->sign;
DUPNUM(*pa,b);
(*pa)->sign *= sign;
}
}
else
{
// b is one so don't divide, but set the sign.
(*pa)->sign *= b->sign;
}
}
//----------------------------------------------------------------------------
//
// FUNCTION: _divnumx
//
// ARGUMENTS: pointer to a number, a second number and precision.
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the number equivalent of *pa /= b.
// Assumes radix is the internal radix representation.
//
//----------------------------------------------------------------------------
void _divnumx( PNUMBER *pa, PNUMBER b, int32_t precision)
{
PNUMBER a= nullptr; // a is the dereferenced number pointer from *pa
PNUMBER c= nullptr; // c will contain the result.
PNUMBER lasttmp = nullptr; // lasttmp allows a backup when the algorithm
// guesses one bit too far.
PNUMBER tmp = nullptr; // current guess being worked on for divide.
PNUMBER rem = nullptr; // remainder after applying guess.
long cdigits; // count of digits for answer.
MANTTYPE *ptrc; // ptrc is a pointer to the mantissa of c.
long thismax = precision + g_ratio; // set a maximum number of internal digits
// to shoot for in the divide.
a=*pa;
if ( thismax < a->cdigit )
{
// a has more digits than precision specified, bump up digits to shoot
// for.
thismax = a->cdigit;
}
if ( thismax < b->cdigit )
{
// b has more digits than precision specified, bump up digits to shoot
// for.
thismax = b->cdigit;
}
// Create c (the divide answer) and set up exponent and sign.
createnum( c, thismax + 1 );
c->exp = (a->cdigit+a->exp) - (b->cdigit+b->exp) + 1;
c->sign = a->sign * b->sign;
ptrc = c->mant + thismax;
cdigits = 0;
DUPNUM( rem, a );
rem->sign = b->sign;
rem->exp = b->cdigit + b->exp - rem->cdigit;
while ( cdigits++ < thismax && !zernum(rem) )
{
long digit = 0;
*ptrc = 0;
while ( !lessnum( rem, b ) )
{
digit = 1;
DUPNUM( tmp, b );
destroynum( lasttmp );
lasttmp=longtonum( 0, BASEX );
while ( lessnum( tmp, rem ) )
{
destroynum( lasttmp );
DUPNUM(lasttmp,tmp);
addnum( &tmp, tmp, BASEX );
digit *= 2;
}
if ( lessnum( rem, tmp ) )
{
// too far, back up...
destroynum( tmp );
digit /= 2;
tmp=lasttmp;
lasttmp= nullptr;
}
tmp->sign *= -1;
addnum( &rem, tmp, BASEX );
destroynum( tmp );
destroynum( lasttmp );
*ptrc |= digit;
}
rem->exp++;
ptrc--;
}
cdigits--;
if ( c->mant != ++ptrc )
{
memmove( c->mant, ptrc, (int)(cdigits*sizeof(MANTTYPE)) );
}
if ( !cdigits )
{
// A zero, make sure no wierd exponents creep in
c->exp = 0;
c->cdigit = 1;
}
else
{
c->cdigit = cdigits;
c->exp -= cdigits;
// prevent different kinds of zeros, by stripping leading duplicate
// zeroes. digits are in order of increasing significance.
while ( c->cdigit > 1 && c->mant[c->cdigit-1] == 0 )
{
c->cdigit--;
}
}
destroynum( rem );
destroynum( *pa );
*pa=c;
}

File diff suppressed because it is too large Load Diff

View File

@@ -0,0 +1,540 @@
// 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;
}

View File

@@ -0,0 +1,258 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
//-----------------------------------------------------------------------------
// Package Title ratpak
// File fact.c
// Copyright (C) 1995-96 Microsoft
// Date 01-16-95
//
//
// Description
//
// Contains fact(orial) and supporting _gamma functions.
//
//-----------------------------------------------------------------------------
#include "pch.h"
#include "ratpak.h"
#define ABSRAT(x) (((x)->pp->sign=1),((x)->pq->sign=1))
#define NEGATE(x) ((x)->pp->sign *= -1)
//-----------------------------------------------------------------------------
//
// FUNCTION: factrat, _gamma, gamma
//
// ARGUMENTS: x PRAT representation of number to take the sine of
//
// RETURN: factorial of x in PRAT form.
//
// EXPLANATION: This uses Taylor series
//
// n
// ___ 2j
// n \ ] A 1 A
// A \ -----[ ---- - ---------------]
// / (2j)! n+2j (n+2j+1)(2j+1)
// /__]
// j=0
//
// / oo
// | n-1 -x __
// This was derived from | x e dx = |
// | | (n) { = (n-1)! for +integers}
// / 0
//
// It can be shown that the above series is within precision if A is chosen
// big enough.
// A n precision
// Based on the relation ne = A 10 A was chosen as
//
// precision
// A = ln(Base /n)+1
// A += n*ln(A) This is close enough for precision > base and n < 1.5
//
//
//-----------------------------------------------------------------------------
void _gamma( PRAT *pn, uint32_t radix, int32_t precision)
{
PRAT factorial= nullptr;
PNUMBER count= nullptr;
PRAT tmp= nullptr;
PRAT one_pt_five= nullptr;
PRAT a= nullptr;
PRAT a2= nullptr;
PRAT term= nullptr;
PRAT sum= nullptr;
PRAT err= nullptr;
PRAT mpy= nullptr;
PRAT ratprec = nullptr;
PRAT ratRadix = nullptr;
long oldprec;
// Set up constants and initial conditions
oldprec = precision;
ratprec = longtorat( oldprec );
// Find the best 'A' for convergence to the required precision.
a=longtorat( radix );
lograt(&a, precision);
mulrat(&a, ratprec, precision);
// Really is -ln(n)+1, but -ln(n) will be < 1
// if we scale n between 0.5 and 1.5
addrat(&a, rat_two, precision);
DUPRAT(tmp,a);
lograt(&tmp, precision);
mulrat(&tmp, *pn, precision);
addrat(&a, tmp, precision);
addrat(&a, rat_one, precision);
// Calculate the necessary bump in precision and up the precision.
// The following code is equivalent to
// precision += ln(exp(a)*pow(a,n+1.5))-ln(radix));
DUPRAT(tmp,*pn);
one_pt_five=longtorat( 3L );
divrat( &one_pt_five, rat_two, precision);
addrat( &tmp, one_pt_five, precision);
DUPRAT(term,a);
powratcomp( &term, tmp, radix, precision);
DUPRAT( tmp, a );
exprat( &tmp, radix, precision);
mulrat( &term, tmp, precision);
lograt( &term, precision);
ratRadix = longtorat(radix);
DUPRAT(tmp,ratRadix);
lograt( &tmp, precision);
subrat( &term, tmp, precision);
precision += rattolong( term, radix, precision);
// Set up initial terms for series, refer to series in above comment block.
DUPRAT(factorial,rat_one); // Start factorial out with one
count = longtonum( 0L, BASEX );
DUPRAT(mpy,a);
powratcomp(&mpy,*pn, radix, precision);
// a2=a^2
DUPRAT(a2,a);
mulrat(&a2, a, precision);
// sum=(1/n)-(a/(n+1))
DUPRAT(sum,rat_one);
divrat(&sum, *pn, precision);
DUPRAT(tmp,*pn);
addrat(&tmp, rat_one, precision);
DUPRAT(term,a);
divrat(&term, tmp, precision);
subrat(&sum, term, precision);
DUPRAT(err,ratRadix);
NEGATE(ratprec);
powratcomp(&err,ratprec, radix, precision);
divrat(&err, ratRadix, precision);
// Just get something not tiny in term
DUPRAT(term, rat_two );
// Loop until precision is reached, or asked to halt.
while ( !zerrat( term ) && rat_gt( term, err, precision) )
{
addrat(pn, rat_two, precision);
// WARNING: mixing numbers and rationals here.
// for speed and efficiency.
INC(count);
mulnumx(&(factorial->pp),count);
INC(count)
mulnumx(&(factorial->pp),count);
divrat(&factorial, a2, precision);
DUPRAT(tmp,*pn);
addrat( &tmp, rat_one, precision);
destroyrat(term);
createrat(term);
DUPNUM(term->pp,count);
DUPNUM(term->pq,num_one);
addrat( &term, rat_one, precision);
mulrat( &term, tmp, precision);
DUPRAT(tmp,a);
divrat( &tmp, term, precision);
DUPRAT(term,rat_one);
divrat( &term, *pn, precision);
subrat( &term, tmp, precision);
divrat (&term, factorial, precision);
addrat( &sum, term, precision);
ABSRAT(term);
}
// Multiply by factor.
mulrat( &sum, mpy, precision);
// And cleanup
precision = oldprec;
destroyrat(ratprec);
destroyrat(err);
destroyrat(term);
destroyrat(a);
destroyrat(a2);
destroyrat(tmp);
destroyrat(one_pt_five);
destroynum(count);
destroyrat(factorial);
destroyrat(*pn);
DUPRAT(*pn,sum);
destroyrat(sum);
}
void factrat( PRAT *px, uint32_t radix, int32_t precision)
{
PRAT fact = nullptr;
PRAT frac = nullptr;
PRAT neg_rat_one = nullptr;
if ( rat_gt( *px, rat_max_fact, precision) || rat_lt( *px, rat_min_fact, precision) )
{
// Don't attempt factorial of anything too large or small.
throw CALC_E_OVERFLOW;
}
DUPRAT(fact,rat_one);
DUPRAT(neg_rat_one,rat_one);
neg_rat_one->pp->sign *= -1;
DUPRAT( frac, *px );
fracrat( &frac, radix, precision);
// Check for negative integers and throw an error.
if ( ( zerrat(frac) || ( LOGRATRADIX(frac) <= -precision) ) &&
( (*px)->pp->sign * (*px)->pq->sign == -1 ) )
{
throw CALC_E_DOMAIN;
}
while ( rat_gt( *px, rat_zero, precision) &&
( LOGRATRADIX(*px) > -precision) )
{
mulrat( &fact, *px, precision);
subrat( px, rat_one, precision);
}
// Added to make numbers 'close enough' to integers use integer factorial.
if ( LOGRATRADIX(*px) <= -precision)
{
DUPRAT((*px),rat_zero);
intrat(&fact, radix, precision);
}
while ( rat_lt( *px, neg_rat_one, precision) )
{
addrat( px, rat_one, precision);
divrat( &fact, *px, precision);
}
if ( rat_neq( *px, rat_zero, precision) )
{
addrat( px, rat_one, precision);
_gamma( px, radix, precision);
mulrat( px, fact, precision);
}
else
{
DUPRAT(*px,fact);
}
destroyrat(fact);
destroyrat(frac);
destroyrat(neg_rat_one);
}

View File

@@ -0,0 +1,342 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
//-----------------------------------------------------------------------------
// Package Title ratpak
// File itrans.c
// Copyright (C) 1995-96 Microsoft
// Date 01-16-95
//
//
// Description
//
// Contains inverse sin, cos, tan functions for rationals
//
// Special Information
//
//-----------------------------------------------------------------------------
#include "pch.h"
#include "ratpak.h"
void ascalerat( _Inout_ PRAT *pa, ANGLE_TYPE angletype, int32_t precision)
{
switch ( angletype )
{
case ANGLE_RAD:
break;
case ANGLE_DEG:
divrat( pa, two_pi, precision);
mulrat( pa, rat_360, precision);
break;
case ANGLE_GRAD:
divrat( pa, two_pi, precision);
mulrat( pa, rat_400, precision);
break;
}
}
//-----------------------------------------------------------------------------
//
// FUNCTION: asinrat, _asinrat
//
// ARGUMENTS: x PRAT representation of number to take the inverse
// sine of
// RETURN: asin of x in PRAT form.
//
// EXPLANATION: This uses Taylor series
//
// n
// ___ 2 2
// \ ] (2j+1) X
// \ thisterm ; where thisterm = thisterm * ---------
// / j j+1 j (2j+2)*(2j+3)
// /__]
// j=0
//
// thisterm = X ; and stop when thisterm < precision used.
// 0 n
//
// If abs(x) > 0.85 then an alternate form is used
// pi/2-sgn(x)*asin(sqrt(1-x^2)
//
//
//-----------------------------------------------------------------------------
void _asinrat( PRAT *px, int32_t precision)
{
CREATETAYLOR();
DUPRAT(pret,*px);
DUPRAT(thisterm,*px);
DUPNUM(n2,num_one);
do
{
NEXTTERM(xx,MULNUM(n2) MULNUM(n2)
INC(n2) DIVNUM(n2) INC(n2) DIVNUM(n2), precision);
}
while ( !SMALL_ENOUGH_RAT( thisterm, precision) );
DESTROYTAYLOR();
}
void asinanglerat( _Inout_ PRAT *pa, ANGLE_TYPE angletype, uint32_t radix, int32_t precision)
{
asinrat( pa, radix, precision);
ascalerat( pa, angletype, precision);
}
void asinrat( PRAT *px, uint32_t radix, int32_t precision)
{
long sgn;
PRAT pret= nullptr;
PRAT phack= nullptr;
sgn = (*px)->pp->sign* (*px)->pq->sign;
(*px)->pp->sign = 1;
(*px)->pq->sign = 1;
// Avoid the really bad part of the asin curve near +/-1.
DUPRAT(phack,*px);
subrat(&phack, rat_one, precision);
// Since *px might be epsilon near zero we must set it to zero.
if ( rat_le(phack, rat_smallest, precision) && rat_ge(phack, rat_negsmallest, precision) )
{
destroyrat(phack);
DUPRAT( *px, pi_over_two );
}
else
{
destroyrat(phack);
if ( rat_gt( *px, pt_eight_five, precision) )
{
if ( rat_gt( *px, rat_one, precision) )
{
subrat( px, rat_one, precision);
if ( rat_gt( *px, rat_smallest, precision) )
{
throw( CALC_E_DOMAIN );
}
else
{
DUPRAT(*px,rat_one);
}
}
DUPRAT(pret,*px);
mulrat( px, pret, precision);
(*px)->pp->sign *= -1;
addrat( px, rat_one, precision);
rootrat( px, rat_two, radix, precision);
_asinrat( px, precision);
(*px)->pp->sign *= -1;
addrat( px, pi_over_two, precision);
destroyrat(pret);
}
else
{
_asinrat( px, precision);
}
}
(*px)->pp->sign = sgn;
(*px)->pq->sign = 1;
}
//-----------------------------------------------------------------------------
//
// FUNCTION: acosrat, _acosrat
//
// ARGUMENTS: x PRAT representation of number to take the inverse
// cosine of
// RETURN: acos of x in PRAT form.
//
// EXPLANATION: This uses Taylor series
//
// n
// ___ 2 2
// \ ] (2j+1) X
// \ thisterm ; where thisterm = thisterm * ---------
// / j j+1 j (2j+2)*(2j+3)
// /__]
// j=0
//
// thisterm = 1 ; and stop when thisterm < precision used.
// 0 n
//
// In this case pi/2-asin(x) is used. At least for now _acosrat isn't
// called.
//
//-----------------------------------------------------------------------------
void acosanglerat( _Inout_ PRAT *pa, ANGLE_TYPE angletype, uint32_t radix, int32_t precision)
{
acosrat( pa, radix, precision);
ascalerat( pa, angletype, precision);
}
void _acosrat( PRAT *px, int32_t precision)
{
CREATETAYLOR();
createrat(thisterm);
thisterm->pp=longtonum( 1L, BASEX );
thisterm->pq=longtonum( 1L, BASEX );
DUPNUM(n2,num_one);
do
{
NEXTTERM(xx,MULNUM(n2) MULNUM(n2)
INC(n2) DIVNUM(n2) INC(n2) DIVNUM(n2), precision);
}
while ( !SMALL_ENOUGH_RAT( thisterm, precision) );
DESTROYTAYLOR();
}
void acosrat( PRAT *px, uint32_t radix, int32_t precision)
{
long sgn;
sgn = (*px)->pp->sign*(*px)->pq->sign;
(*px)->pp->sign = 1;
(*px)->pq->sign = 1;
if ( rat_equ( *px, rat_one, precision) )
{
if ( sgn == -1 )
{
DUPRAT(*px,pi);
}
else
{
DUPRAT( *px, rat_zero );
}
}
else
{
(*px)->pp->sign = sgn;
asinrat( px, radix, precision);
(*px)->pp->sign *= -1;
addrat(px, pi_over_two, precision);
}
}
//-----------------------------------------------------------------------------
//
// FUNCTION: atanrat, _atanrat
//
// ARGUMENTS: x PRAT representation of number to take the inverse
// hyperbolic tangent of
//
// RETURN: atanh of x in PRAT form.
//
// EXPLANATION: This uses Taylor series
//
// n
// ___ 2
// \ ] (2j)*X (-1^j)
// \ thisterm ; where thisterm = thisterm * ---------
// / j j+1 j (2j+2)
// /__]
// j=0
//
// thisterm = X ; and stop when thisterm < precision used.
// 0 n
//
// If abs(x) > 0.85 then an alternate form is used
// asin(x/sqrt(q+x^2))
//
// And if abs(x) > 2.0 then this form is used.
//
// pi/2 - atan(1/x)
//
//-----------------------------------------------------------------------------
void atananglerat( _Inout_ PRAT *pa, ANGLE_TYPE angletype, uint32_t radix, int32_t precision)
{
atanrat( pa, radix, precision);
ascalerat( pa, angletype, precision);
}
void _atanrat( PRAT *px, int32_t precision)
{
CREATETAYLOR();
DUPRAT(pret,*px);
DUPRAT(thisterm,*px);
DUPNUM(n2,num_one);
xx->pp->sign *= -1;
do {
NEXTTERM(xx,MULNUM(n2) INC(n2) INC(n2) DIVNUM(n2), precision);
} while ( !SMALL_ENOUGH_RAT( thisterm, precision) );
DESTROYTAYLOR();
}
void atanrat( PRAT *px, uint32_t radix, int32_t precision)
{
long sgn;
PRAT tmpx= nullptr;
sgn = (*px)->pp->sign * (*px)->pq->sign;
(*px)->pp->sign = 1;
(*px)->pq->sign = 1;
if ( rat_gt( (*px), pt_eight_five, precision) )
{
if ( rat_gt( (*px), rat_two, precision) )
{
(*px)->pp->sign = sgn;
(*px)->pq->sign = 1;
DUPRAT(tmpx,rat_one);
divrat(&tmpx, (*px), precision);
_atanrat(&tmpx, precision);
tmpx->pp->sign = sgn;
tmpx->pq->sign = 1;
DUPRAT(*px,pi_over_two);
subrat(px, tmpx, precision);
destroyrat( tmpx );
}
else
{
(*px)->pp->sign = sgn;
DUPRAT(tmpx,*px);
mulrat( &tmpx, *px, precision);
addrat( &tmpx, rat_one, precision);
rootrat( &tmpx, rat_two, radix, precision);
divrat( px, tmpx, precision);
destroyrat( tmpx );
asinrat( px, radix, precision);
(*px)->pp->sign = sgn;
(*px)->pq->sign = 1;
}
}
else
{
(*px)->pp->sign = sgn;
(*px)->pq->sign = 1;
_atanrat( px, precision);
}
if ( rat_gt( *px, pi_over_two, precision) )
{
subrat( px, pi, precision);
}
}

View File

@@ -0,0 +1,160 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
//-----------------------------------------------------------------------------
// Package Title ratpak
// File itransh.c
// Copyright (C) 1995-97 Microsoft
// Date 01-16-95
//
//
// Description
//
// Contains inverse hyperbolic sin, cos, and tan functions.
//
// Special Information
//
//
//-----------------------------------------------------------------------------
#include "pch.h"
#include "ratpak.h"
//-----------------------------------------------------------------------------
//
// FUNCTION: asinhrat
//
// ARGUMENTS: x PRAT representation of number to take the inverse
// hyperbolic sine of
// RETURN: asinh of x in PRAT form.
//
// EXPLANATION: This uses Taylor series
//
// n
// ___ 2 2
// \ ] -(2j+1) X
// \ thisterm ; where thisterm = thisterm * ---------
// / j j+1 j (2j+2)*(2j+3)
// /__]
// j=0
//
// thisterm = X ; and stop when thisterm < precision used.
// 0 n
//
// For abs(x) < .85, and
//
// asinh(x) = log(x+sqrt(x^2+1))
//
// For abs(x) >= .85
//
//-----------------------------------------------------------------------------
void asinhrat( PRAT *px, uint32_t radix, int32_t precision)
{
PRAT neg_pt_eight_five = nullptr;
DUPRAT(neg_pt_eight_five,pt_eight_five);
neg_pt_eight_five->pp->sign *= -1;
if ( rat_gt( *px, pt_eight_five, precision) || rat_lt( *px, neg_pt_eight_five, precision) )
{
PRAT ptmp = nullptr;
DUPRAT(ptmp,(*px));
mulrat(&ptmp, *px, precision);
addrat(&ptmp, rat_one, precision);
rootrat(&ptmp, rat_two, radix, precision);
addrat(px, ptmp, precision);
lograt(px, precision);
destroyrat(ptmp);
}
else
{
CREATETAYLOR();
xx->pp->sign *= -1;
DUPRAT(pret,(*px));
DUPRAT(thisterm,(*px));
DUPNUM(n2,num_one);
do
{
NEXTTERM(xx,MULNUM(n2) MULNUM(n2)
INC(n2) DIVNUM(n2) INC(n2) DIVNUM(n2), precision);
}
while ( !SMALL_ENOUGH_RAT( thisterm, precision) );
DESTROYTAYLOR();
}
destroyrat(neg_pt_eight_five);
}
//-----------------------------------------------------------------------------
//
// FUNCTION: acoshrat
//
// ARGUMENTS: x PRAT representation of number to take the inverse
// hyperbolic cose of
// RETURN: acosh of x in PRAT form.
//
// EXPLANATION: This uses
//
// acosh(x)=ln(x+sqrt(x^2-1))
//
// For x >= 1
//
//-----------------------------------------------------------------------------
void acoshrat( PRAT *px, uint32_t radix, int32_t precision)
{
if ( rat_lt( *px, rat_one, precision) )
{
throw CALC_E_DOMAIN;
}
else
{
PRAT ptmp = nullptr;
DUPRAT(ptmp,(*px));
mulrat(&ptmp, *px, precision);
subrat(&ptmp, rat_one, precision);
rootrat(&ptmp,rat_two, radix, precision);
addrat(px, ptmp, precision);
lograt(px, precision);
destroyrat(ptmp);
}
}
//-----------------------------------------------------------------------------
//
// FUNCTION: atanhrat
//
// ARGUMENTS: x PRAT representation of number to take the inverse
// hyperbolic tangent of
//
// RETURN: atanh of x in PRAT form.
//
// EXPLANATION: This uses
//
// 1 x+1
// atanh(x) = -*ln(----)
// 2 x-1
//
//-----------------------------------------------------------------------------
void atanhrat( PRAT *px, int32_t precision)
{
PRAT ptmp = nullptr;
DUPRAT(ptmp,(*px));
subrat(&ptmp, rat_one, precision);
addrat(px, rat_one, precision);
divrat(px, ptmp, precision);
(*px)->pp->sign *= -1;
lograt(px, precision);
divrat(px, rat_two, precision);
destroyrat(ptmp);
}

View File

@@ -0,0 +1,219 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
//---------------------------------------------------------------------------
// Package Title ratpak
// File num.c
// Copyright (C) 1995-99 Microsoft
// Date 01-16-95
//
//
// Description
//
// Contains routines for and, or, xor, not and other support
//
//---------------------------------------------------------------------------
#include "pch.h"
#include "ratpak.h"
void lshrat( PRAT *pa, PRAT b, uint32_t radix, int32_t precision)
{
PRAT pwr= nullptr;
long intb;
intrat(pa, radix, precision);
if ( !zernum( (*pa)->pp ) )
{
// If input is zero we're done.
if ( rat_gt( b, rat_max_exp, precision) )
{
// Don't attempt lsh of anything big
throw( CALC_E_DOMAIN );
}
intb = rattolong(b, radix, precision);
DUPRAT(pwr,rat_two);
ratpowlong(&pwr, intb, precision);
mulrat(pa, pwr, precision);
destroyrat(pwr);
}
}
void rshrat( PRAT *pa, PRAT b, uint32_t radix, int32_t precision)
{
PRAT pwr= nullptr;
long intb;
intrat(pa, radix, precision);
if ( !zernum( (*pa)->pp ) )
{
// If input is zero we're done.
if ( rat_lt( b, rat_min_exp, precision) )
{
// Don't attempt rsh of anything big and negative.
throw( CALC_E_DOMAIN );
}
intb = rattolong(b, radix, precision);
DUPRAT(pwr,rat_two);
ratpowlong(&pwr, intb, precision);
divrat(pa, pwr, precision);
destroyrat(pwr);
}
}
void boolrat( PRAT *pa, PRAT b, int func, uint32_t radix, int32_t precision);
void boolnum( PNUMBER *pa, PNUMBER b, int func );
enum {
FUNC_AND,
FUNC_OR,
FUNC_XOR
} BOOL_FUNCS;
void andrat( PRAT *pa, PRAT b, uint32_t radix, int32_t precision)
{
boolrat( pa, b, FUNC_AND, radix, precision);
}
void orrat( PRAT *pa, PRAT b, uint32_t radix, int32_t precision)
{
boolrat( pa, b, FUNC_OR, radix, precision);
}
void xorrat( PRAT *pa, PRAT b, uint32_t radix, int32_t precision)
{
boolrat( pa, b, FUNC_XOR, radix, precision);
}
//---------------------------------------------------------------------------
//
// FUNCTION: boolrat
//
// ARGUMENTS: pointer to a rational a second rational.
//
// RETURN: None, changes pointer.
//
// DESCRIPTION: Does the rational equivalent of *pa op= b;
//
//---------------------------------------------------------------------------
void boolrat( PRAT *pa, PRAT b, int func, uint32_t radix, int32_t precision)
{
PRAT tmp= nullptr;
intrat( pa, radix, precision);
DUPRAT(tmp,b);
intrat( &tmp, radix, precision);
boolnum( &((*pa)->pp), tmp->pp, func );
destroyrat(tmp);
}
//---------------------------------------------------------------------------
//
// FUNCTION: boolnum
//
// ARGUMENTS: pointer to a number a second number
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the number equivalent of *pa &= b.
// radix doesn't matter for logicals.
// WARNING: Assumes numbers are unsigned.
//
//---------------------------------------------------------------------------
void boolnum( PNUMBER *pa, PNUMBER b, int func )
{
PNUMBER c= nullptr;
PNUMBER a= nullptr;
MANTTYPE *pcha;
MANTTYPE *pchb;
MANTTYPE *pchc;
long cdigits;
long mexp;
MANTTYPE da;
MANTTYPE db;
a=*pa;
cdigits = max( a->cdigit+a->exp, b->cdigit+b->exp ) -
min( a->exp, b->exp );
createnum( c, cdigits );
c->exp = min( a->exp, b->exp );
mexp = c->exp;
c->cdigit = cdigits;
pcha = a->mant;
pchb = b->mant;
pchc = c->mant;
for ( ;cdigits > 0; cdigits--, mexp++ )
{
da = ( ( ( mexp >= a->exp ) && ( cdigits + a->exp - c->exp >
(c->cdigit - a->cdigit) ) ) ?
*pcha++ : 0 );
db = ( ( ( mexp >= b->exp ) && ( cdigits + b->exp - c->exp >
(c->cdigit - b->cdigit) ) ) ?
*pchb++ : 0 );
switch ( func )
{
case FUNC_AND:
*pchc++ = da & db;
break;
case FUNC_OR:
*pchc++ = da | db;
break;
case FUNC_XOR:
*pchc++ = da ^ db;
break;
}
}
c->sign = a->sign;
while ( c->cdigit > 1 && *(--pchc) == 0 )
{
c->cdigit--;
}
destroynum( *pa );
*pa=c;
}
//-----------------------------------------------------------------------------
//
// FUNCTION: modrat
//
// ARGUMENTS: pointer to a rational a second rational.
//
// RETURN: None, changes pointer.
//
// DESCRIPTION: Does the rational equivalent of frac(*pa);
//
//-----------------------------------------------------------------------------
void modrat( PRAT *pa, PRAT b )
{
PRAT tmp = nullptr;
if ( zerrat( b ) )
{
throw CALC_E_INDEFINITE;
}
DUPRAT(tmp,b);
mulnumx( &((*pa)->pp), tmp->pq );
mulnumx( &(tmp->pp), (*pa)->pq );
remnum( &((*pa)->pp), tmp->pp, BASEX );
mulnumx( &((*pa)->pq), tmp->pq );
//Get *pa back in the integer over integer form.
RENORMALIZE(*pa);
destroyrat( tmp );
}

View File

@@ -0,0 +1,655 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
//-----------------------------------------------------------------------------
// Package Title ratpak
// File num.c
// Copyright (C) 1995-97 Microsoft
// Date 01-16-95
//
//
// Description
//
// Contains number routines for add, mul, div, rem and other support
// and longs.
//
// Special Information
//
//
//-----------------------------------------------------------------------------
#include "pch.h"
#include "ratpak.h"
using namespace std;
//----------------------------------------------------------------------------
//
// FUNCTION: addnum
//
// ARGUMENTS: pointer to a number a second number, and the
// radix.
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the number equivalent of *pa += b.
// Assumes radix is the base of both numbers.
//
// ALGORITHM: Adds each digit from least significant to most
// significant.
//
//
//----------------------------------------------------------------------------
void _addnum( PNUMBER *pa, PNUMBER b, uint32_t radix);
void __inline addnum( PNUMBER *pa, PNUMBER b, uint32_t radix)
{
if ( b->cdigit > 1 || b->mant[0] != 0 )
{ // If b is zero we are done.
if ( (*pa)->cdigit > 1 || (*pa)->mant[0] != 0 )
{ // pa and b are both nonzero.
_addnum( pa, b, radix);
}
else
{ // if pa is zero and b isn't just copy b.
DUPNUM(*pa,b);
}
}
}
void _addnum( PNUMBER *pa, PNUMBER b, uint32_t radix)
{
PNUMBER c= nullptr; // c will contain the result.
PNUMBER a= nullptr; // a is the dereferenced number pointer from *pa
MANTTYPE *pcha; // pcha is a pointer to the mantissa of a.
MANTTYPE *pchb; // pchb is a pointer to the mantissa of b.
MANTTYPE *pchc; // pchc is a pointer to the mantissa of c.
long cdigits; // cdigits is the max count of the digits results
// used as a counter.
long mexp; // mexp is the exponent of the result.
MANTTYPE da; // da is a single 'digit' after possible padding.
MANTTYPE db; // db is a single 'digit' after possible padding.
MANTTYPE cy=0; // cy is the value of a carry after adding two 'digits'
long fcompla = 0; // fcompla is a flag to signal a is negative.
long fcomplb = 0; // fcomplb is a flag to signal b is negative.
a=*pa;
// Calculate the overlap of the numbers after alignment, this includes
// necessary padding 0's
cdigits = max( a->cdigit+a->exp, b->cdigit+b->exp ) -
min( a->exp, b->exp );
createnum( c, cdigits + 1 );
c->exp = min( a->exp, b->exp );
mexp = c->exp;
c->cdigit = cdigits;
pcha = a->mant;
pchb = b->mant;
pchc = c->mant;
// Figure out the sign of the numbers
if ( a->sign != b->sign )
{
cy = 1;
fcompla = ( a->sign == -1 );
fcomplb = ( b->sign == -1 );
}
// Loop over all the digits, real and 0 padded. Here we know a and b are
// aligned
for ( ;cdigits > 0; cdigits--, mexp++ )
{
// Get digit from a, taking padding into account.
da = ( ( ( mexp >= a->exp ) && ( cdigits + a->exp - c->exp >
(c->cdigit - a->cdigit) ) ) ?
*pcha++ : 0 );
// Get digit from b, taking padding into account.
db = ( ( ( mexp >= b->exp ) && ( cdigits + b->exp - c->exp >
(c->cdigit - b->cdigit) ) ) ?
*pchb++ : 0 );
// Handle complementing for a and b digit. Might be a better way, but
// haven't found it yet.
if ( fcompla )
{
da = (MANTTYPE)(radix) - 1 - da;
}
if ( fcomplb )
{
db = (MANTTYPE)(radix) - 1 - db;
}
// Update carry as necessary
cy = da + db + cy;
*pchc++ = (MANTTYPE)(cy % (MANTTYPE)radix);
cy /= (MANTTYPE)radix;
}
// Handle carry from last sum as extra digit
if ( cy && !(fcompla || fcomplb) )
{
*pchc++ = cy;
c->cdigit++;
}
// Compute sign of result
if ( !(fcompla || fcomplb) )
{
c->sign = a->sign;
}
else
{
if ( cy )
{
c->sign = 1;
}
else
{
// In this particular case an overflow or underflow has occoured
// and all the digits need to be complemented, at one time an
// attempt to handle this above was made, it turned out to be much
// slower on average.
c->sign = -1;
cy = 1;
for ( ( cdigits = c->cdigit ), (pchc = c->mant);
cdigits > 0;
cdigits-- )
{
cy = (MANTTYPE)radix - (MANTTYPE)1 - *pchc + cy;
*pchc++ = (MANTTYPE)( cy % (MANTTYPE)radix);
cy /= (MANTTYPE)radix;
}
}
}
// Remove leading zeroes, remember digits are in order of
// increasing significance. i.e. 100 would be 0,0,1
while ( c->cdigit > 1 && *(--pchc) == 0 )
{
c->cdigit--;
}
destroynum( *pa );
*pa=c;
}
//----------------------------------------------------------------------------
//
// FUNCTION: mulnum
//
// ARGUMENTS: pointer to a number a second number, and the
// radix.
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the number equivalent of *pa *= b.
// Assumes radix is the radix of both numbers. This algorithm is the
// same one you learned in gradeschool.
//
//----------------------------------------------------------------------------
void _mulnum( PNUMBER *pa, PNUMBER b, uint32_t radix);
void __inline mulnum( PNUMBER *pa, PNUMBER b, uint32_t radix)
{
if ( b->cdigit > 1 || b->mant[0] != 1 || b->exp != 0 )
{ // If b is one we don't multiply exactly.
if ( (*pa)->cdigit > 1 || (*pa)->mant[0] != 1 || (*pa)->exp != 0 )
{ // pa and b are both nonone.
_mulnum( pa, b, radix);
}
else
{ // if pa is one and b isn't just copy b, and adjust the sign.
long sign = (*pa)->sign;
DUPNUM(*pa,b);
(*pa)->sign *= sign;
}
}
else
{ // But we do have to set the sign.
(*pa)->sign *= b->sign;
}
}
void _mulnum( PNUMBER *pa, PNUMBER b, uint32_t radix)
{
PNUMBER c= nullptr; // c will contain the result.
PNUMBER a= nullptr; // a is the dereferenced number pointer from *pa
MANTTYPE *pcha; // pcha is a pointer to the mantissa of a.
MANTTYPE *pchb; // pchb is a pointer to the mantissa of b.
MANTTYPE *pchc; // pchc is a pointer to the mantissa of c.
MANTTYPE *pchcoffset; // pchcoffset, is the anchor location of the next
// single digit multiply partial result.
long iadigit = 0; // Index of digit being used in the first number.
long ibdigit = 0; // Index of digit being used in the second number.
MANTTYPE da = 0; // da is the digit from the fist number.
TWO_MANTTYPE cy = 0; // cy is the carry resulting from the addition of
// a multiplied row into the result.
TWO_MANTTYPE mcy = 0; // mcy is the resultant from a single
// multiply, AND the carry of that multiply.
long icdigit = 0; // Index of digit being calculated in final result.
a=*pa;
ibdigit = a->cdigit + b->cdigit - 1;
createnum( c, ibdigit + 1 );
c->cdigit = ibdigit;
c->sign = a->sign * b->sign;
c->exp = a->exp + b->exp;
pcha = a->mant;
pchcoffset = c->mant;
for ( iadigit = a->cdigit; iadigit > 0; iadigit-- )
{
da = *pcha++;
pchb = b->mant;
// Shift pchc, and pchcoffset, one for each digit
pchc = pchcoffset++;
for ( ibdigit = b->cdigit; ibdigit > 0; ibdigit-- )
{
cy = 0;
mcy = (TWO_MANTTYPE)da * *pchb;
if ( mcy )
{
icdigit = 0;
if ( ibdigit == 1 && iadigit == 1 )
{
c->cdigit++;
}
}
// If result is nonzero, or while result of carry is nonzero...
while ( mcy || cy )
{
// update carry from addition(s) and multiply.
cy += (TWO_MANTTYPE)pchc[icdigit]+(mcy%(TWO_MANTTYPE)radix);
// update result digit from
pchc[icdigit++]=(MANTTYPE)(cy%(TWO_MANTTYPE)radix);
// update carries from
mcy /= (TWO_MANTTYPE)radix;
cy /= (TWO_MANTTYPE)radix;
}
pchb++;
pchc++;
}
}
// prevent different kinds of zeros, by stripping leading duplicate zeroes.
// digits are in order of increasing significance.
while ( c->cdigit > 1 && c->mant[c->cdigit-1] == 0 )
{
c->cdigit--;
}
destroynum( *pa );
*pa=c;
}
//----------------------------------------------------------------------------
//
// FUNCTION: remnum
//
// ARGUMENTS: pointer to a number a second number, and the
// radix.
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the number equivalent of *pa %= b.
// Repeatedly subtracts off powers of 2 of b until *pa < b.
//
//
//----------------------------------------------------------------------------
void remnum( PNUMBER *pa, PNUMBER b, uint32_t radix)
{
PNUMBER tmp = nullptr; // tmp is the working remainder.
PNUMBER lasttmp = nullptr; // lasttmp is the last remainder which worked.
// Once *pa is less than b, *pa is the remainder.
while ( !lessnum( *pa, b ) )
{
DUPNUM( tmp, b );
if ( lessnum( tmp, *pa ) )
{
// Start off close to the right answer for subtraction.
tmp->exp = (*pa)->cdigit+(*pa)->exp - tmp->cdigit;
if ( MSD(*pa) <= MSD(tmp) )
{
// Don't take the chance that the numbers are equal.
tmp->exp--;
}
}
destroynum( lasttmp );
lasttmp=longtonum( 0, radix);
while ( lessnum( tmp, *pa ) )
{
DUPNUM( lasttmp, tmp );
addnum( &tmp, tmp, radix);
}
if ( lessnum( *pa, tmp ) )
{
// too far, back up...
destroynum( tmp );
tmp=lasttmp;
lasttmp= nullptr;
}
// Subtract the working remainder from the remainder holder.
tmp->sign = -1*(*pa)->sign;
addnum( pa, tmp, radix);
destroynum( tmp );
destroynum( lasttmp );
}
}
//---------------------------------------------------------------------------
//
// FUNCTION: divnum
//
// ARGUMENTS: pointer to a number a second number, and the
// radix.
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the number equivalent of *pa /= b.
// Assumes radix is the radix of both numbers.
//
//---------------------------------------------------------------------------
void _divnum( PNUMBER *pa, PNUMBER b, uint32_t radix, int32_t precision);
void __inline divnum( PNUMBER *pa, PNUMBER b, uint32_t radix, int32_t precision)
{
if ( b->cdigit > 1 || b->mant[0] != 1 || b->exp != 0 )
{
// b is not one
_divnum( pa, b, radix, precision);
}
else
{ // But we do have to set the sign.
(*pa)->sign *= b->sign;
}
}
void _divnum( PNUMBER *pa, PNUMBER b, uint32_t radix, int32_t precision)
{
PNUMBER a = *pa;
long thismax = precision + 2;
if (thismax < a->cdigit)
{
thismax = a->cdigit;
}
if (thismax < b->cdigit)
{
thismax = b->cdigit;
}
PNUMBER c = nullptr;
createnum(c, thismax + 1);
c->exp = (a->cdigit + a->exp) - (b->cdigit + b->exp) + 1;
c->sign = a->sign * b->sign;
MANTTYPE *ptrc = c->mant + thismax;
PNUMBER rem = nullptr;
PNUMBER tmp = nullptr;
DUPNUM(rem, a);
DUPNUM(tmp, b);
tmp->sign = a->sign;
rem->exp = b->cdigit + b->exp - rem->cdigit;
// Build a table of multiplications of the divisor, this is quicker for
// more than radix 'digits'
list<PNUMBER> numberList{ longtonum(0L, radix) };
for (unsigned long i = 1; i < radix; i++)
{
PNUMBER newValue = nullptr;
DUPNUM(newValue, numberList.front());
addnum(&newValue, tmp, radix);
numberList.emplace_front(newValue);
}
destroynum(tmp);
long digit;
long cdigits = 0;
while (cdigits++ < thismax && !zernum(rem))
{
digit = radix - 1;
PNUMBER multiple = nullptr;
for (const auto& num : numberList)
{
if (!lessnum(rem, num) || !--digit)
{
multiple = num;
break;
}
}
if (digit)
{
multiple->sign *= -1;
addnum(&rem, multiple, radix);
multiple->sign *= -1;
}
rem->exp++;
*ptrc-- = (MANTTYPE)digit;
}
cdigits--;
if (c->mant != ++ptrc)
{
memmove(c->mant, ptrc, (int)(cdigits * sizeof(MANTTYPE)));
}
// Cleanup table structure
for (auto& num : numberList)
{
destroynum(num);
}
if (!cdigits)
{
c->cdigit = 1;
c->exp = 0;
}
else
{
c->cdigit = cdigits;
c->exp -= cdigits;
while (c->cdigit > 1 && c->mant[c->cdigit - 1] == 0)
{
c->cdigit--;
}
}
destroynum(rem);
destroynum(*pa);
*pa = c;
}
//---------------------------------------------------------------------------
//
// FUNCTION: equnum
//
// ARGUMENTS: two numbers.
//
// RETURN: Boolean
//
// DESCRIPTION: Does the number equivalent of ( a == b )
// Only assumes that a and b are the same radix.
//
//---------------------------------------------------------------------------
bool equnum( PNUMBER a, PNUMBER b )
{
long diff;
MANTTYPE *pa;
MANTTYPE *pb;
long cdigits;
long ccdigits;
MANTTYPE da;
MANTTYPE db;
diff = ( a->cdigit + a->exp ) - ( b->cdigit + b->exp );
if ( diff < 0 )
{
// If the exponents are different, these are different numbers.
return false;
}
else
{
if ( diff > 0 )
{
// If the exponents are different, these are different numbers.
return false;
}
else
{
// OK the exponents match.
pa = a->mant;
pb = b->mant;
pa += a->cdigit - 1;
pb += b->cdigit - 1;
cdigits = max( a->cdigit, b->cdigit );
ccdigits = cdigits;
// Loop over all digits until we run out of digits or there is a
// difference in the digits.
for ( ;cdigits > 0; cdigits-- )
{
da = ( (cdigits > (ccdigits - a->cdigit) ) ?
*pa-- : 0 );
db = ( (cdigits > (ccdigits - b->cdigit) ) ?
*pb-- : 0 );
if ( da != db )
{
return false;
}
}
// In this case, they are equal.
return true;
}
}
}
//---------------------------------------------------------------------------
//
// FUNCTION: lessnum
//
// ARGUMENTS: two numbers.
//
// RETURN: Boolean
//
// DESCRIPTION: Does the number equivalent of ( abs(a) < abs(b) )
// Only assumes that a and b are the same radix, WARNING THIS IS AN.
// UNSIGNED COMPARE!
//
//---------------------------------------------------------------------------
bool lessnum( PNUMBER a, PNUMBER b )
{
long diff;
MANTTYPE *pa;
MANTTYPE *pb;
long cdigits;
long ccdigits;
MANTTYPE da;
MANTTYPE db;
diff = ( a->cdigit + a->exp ) - ( b->cdigit + b->exp );
if ( diff < 0 )
{
// The exponent of a is less than b
return true;
}
else
{
if ( diff > 0 )
{
return false;
}
else
{
pa = a->mant;
pb = b->mant;
pa += a->cdigit - 1;
pb += b->cdigit - 1;
cdigits = max( a->cdigit, b->cdigit );
ccdigits = cdigits;
for ( ;cdigits > 0; cdigits-- )
{
da = ( (cdigits > (ccdigits - a->cdigit) ) ?
*pa-- : 0 );
db = ( (cdigits > (ccdigits - b->cdigit) ) ?
*pb-- : 0 );
diff = da-db;
if ( diff )
{
return( diff < 0 );
}
}
// In this case, they are equal.
return false;
}
}
}
//----------------------------------------------------------------------------
//
// FUNCTION: zernum
//
// ARGUMENTS: number
//
// RETURN: Boolean
//
// DESCRIPTION: Does the number equivalent of ( !a )
//
//----------------------------------------------------------------------------
bool zernum( PNUMBER a )
{
long length;
MANTTYPE *pcha;
length = a->cdigit;
pcha = a->mant;
// loop over all the digits until you find a nonzero or until you run
// out of digits
while ( length-- > 0 )
{
if ( *pcha++ )
{
// One of the digits isn't zero, therefore the number isn't zero
return false;
}
}
// All of the digits are zero, therefore the number is zero
return true;
}

View File

@@ -0,0 +1,303 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
//-----------------------------------------------------------------------------
// Package Title ratpak
// File rat.c
// Copyright (C) 1995-96 Microsoft
// Date 01-16-95
//
//
// Description
//
// Contains mul, div, add, and other support functions for rationals.
//
//
//
//-----------------------------------------------------------------------------
#include "pch.h"
#include "ratpak.h"
using namespace std;
//-----------------------------------------------------------------------------
//
// FUNCTION: gcdrat
//
// ARGUMENTS: pointer to a rational.
//
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Divides p and q in rational by the G.C.D.
// of both. It was hoped this would speed up some
// calculations, and until the above trimming was done it
// did, but after trimming gcdratting, only slows things
// down.
//
//-----------------------------------------------------------------------------
void gcdrat( PRAT *pa, uint32_t radix, int32_t precision)
{
PNUMBER pgcd= nullptr;
PRAT a= nullptr;
a=*pa;
pgcd = gcd( a->pp, a->pq );
if ( !zernum( pgcd ) )
{
divnumx( &(a->pp), pgcd, precision);
divnumx( &(a->pq), pgcd, precision);
}
destroynum( pgcd );
*pa=a;
RENORMALIZE(*pa);
}
//-----------------------------------------------------------------------------
//
// FUNCTION: fracrat
//
// ARGUMENTS: pointer to a rational a second rational.
//
// RETURN: None, changes pointer.
//
// DESCRIPTION: Does the rational equivalent of frac(*pa);
//
//-----------------------------------------------------------------------------
void fracrat( PRAT *pa , uint32_t radix, int32_t precision)
{
// Only do the intrat operation if number is nonzero.
// and only if the bottom part is not one.
if ( !zernum( (*pa)->pp ) && !equnum( (*pa)->pq, num_one ) )
{
wstring ratStr = RatToString(*pa, FMT_FLOAT, radix, precision);
PNUMBER pnum = StringToNumber(ratStr, radix, precision);
destroyrat( *pa );
*pa = numtorat( pnum, radix);
destroynum( pnum );
}
remnum( &((*pa)->pp), (*pa)->pq, BASEX );
//Get *pa back in the integer over integer form.
RENORMALIZE(*pa);
}
//-----------------------------------------------------------------------------
//
// FUNCTION: mulrat
//
// ARGUMENTS: pointer to a rational a second rational.
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the rational equivalent of *pa *= b.
// Assumes radix is the radix of both numbers.
//
//-----------------------------------------------------------------------------
void mulrat( PRAT *pa, PRAT b, int32_t precision)
{
// Only do the multiply if it isn't zero.
if ( !zernum( (*pa)->pp ) )
{
mulnumx( &((*pa)->pp), b->pp );
mulnumx( &((*pa)->pq), b->pq );
trimit(pa, precision);
}
else
{
// If it is zero, blast a one in the denominator.
DUPNUM( ((*pa)->pq), num_one );
}
#ifdef MULGCD
gcdrat( pa );
#endif
}
//-----------------------------------------------------------------------------
//
// FUNCTION: divrat
//
// ARGUMENTS: pointer to a rational a second rational.
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the rational equivalent of *pa /= b.
// Assumes radix is the radix of both numbers.
//
//-----------------------------------------------------------------------------
void divrat( PRAT *pa, PRAT b, int32_t precision)
{
if ( !zernum( (*pa)->pp ) )
{
// Only do the divide if the top isn't zero.
mulnumx( &((*pa)->pp), b->pq );
mulnumx( &((*pa)->pq), b->pp );
if ( zernum( (*pa)->pq ) )
{
// raise an exception if the bottom is 0.
throw( CALC_E_DIVIDEBYZERO );
}
trimit(pa, precision);
}
else
{
// Top is zero.
if ( zerrat( b ) )
{
// If bottom is zero
// 0 / 0 is indefinite, raise an exception.
throw( CALC_E_INDEFINITE );
}
else
{
// 0/x make a unique 0.
DUPNUM( ((*pa)->pq), num_one );
}
}
#ifdef DIVGCD
gcdrat( pa );
#endif
}
//-----------------------------------------------------------------------------
//
// FUNCTION: subrat
//
// ARGUMENTS: pointer to a rational a second rational.
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the rational equivalent of *pa += b.
// Assumes base is internal througought.
//
//-----------------------------------------------------------------------------
void subrat( PRAT *pa, PRAT b, int32_t precision)
{
b->pp->sign *= -1;
addrat( pa, b, precision);
b->pp->sign *= -1;
}
//-----------------------------------------------------------------------------
//
// FUNCTION: addrat
//
// ARGUMENTS: pointer to a rational a second rational.
//
// RETURN: None, changes first pointer.
//
// DESCRIPTION: Does the rational equivalent of *pa += b.
// Assumes base is internal througought.
//
//-----------------------------------------------------------------------------
void addrat( PRAT *pa, PRAT b, int32_t precision)
{
PNUMBER bot= nullptr;
if ( equnum( (*pa)->pq, b->pq ) )
{
// Very special case, q's match.,
// make sure signs are involved in the calculation
// we have to do this since the optimization here is only
// working with the top half of the rationals.
(*pa)->pp->sign *= (*pa)->pq->sign;
(*pa)->pq->sign = 1;
b->pp->sign *= b->pq->sign;
b->pq->sign = 1;
addnum( &((*pa)->pp), b->pp, BASEX );
}
else
{
// Usual case q's aren't the same.
DUPNUM( bot, (*pa)->pq );
mulnumx( &bot, b->pq );
mulnumx( &((*pa)->pp), b->pq );
mulnumx( &((*pa)->pq), b->pp );
addnum( &((*pa)->pp), (*pa)->pq, BASEX );
destroynum( (*pa)->pq );
(*pa)->pq = bot;
trimit(pa, precision);
// Get rid of negative zeroes here.
(*pa)->pp->sign *= (*pa)->pq->sign;
(*pa)->pq->sign = 1;
}
#ifdef ADDGCD
gcdrat( pa );
#endif
}
//-----------------------------------------------------------------------------
//
// FUNCTION: rootrat
//
// PARAMETERS: y prat representation of number to take the root of
// n prat representation of the root to take.
//
// RETURN: bth root of a in rat form.
//
// EXPLANATION: This is now a stub function to powrat().
//
//-----------------------------------------------------------------------------
void rootrat( PRAT *py, PRAT n, uint32_t radix, int32_t precision)
{
// Initialize 1/n
PRAT oneovern= nullptr;
DUPRAT(oneovern,rat_one);
divrat(&oneovern, n, precision);
powrat(py, oneovern, radix, precision);
destroyrat(oneovern);
}
//-----------------------------------------------------------------------------
//
// FUNCTION: zerrat
//
// ARGUMENTS: Rational number.
//
// RETURN: Boolean
//
// DESCRIPTION: Returns true if input is zero.
// False otherwise.
//
//-----------------------------------------------------------------------------
bool zerrat( PRAT a )
{
return( zernum(a->pp) );
}

View File

@@ -0,0 +1,482 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
// Autogenerated by _dumprawrat in support.c
NUMBER init_num_one= {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_num_two= {
1,
1,
0,
{ 2,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_num_five= {
1,
1,
0,
{ 5,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_num_six= {
1,
1,
0,
{ 6,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_num_ten= {
1,
1,
0,
{ 10,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_smallest = {
1,
1,
0,
{ 1,}
};
NUMBER init_q_rat_smallest = {
1,
4,
0,
{ 0, 190439170, 901055854, 10097,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_negsmallest = {
-1,
1,
0,
{ 1,}
};
NUMBER init_q_rat_negsmallest = {
1,
4,
0,
{ 0, 190439170, 901055854, 10097,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_pt_eight_five = {
1,
1,
0,
{ 85,}
};
NUMBER init_q_pt_eight_five = {
1,
1,
0,
{ 100,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_six = {
1,
1,
0,
{ 6,}
};
NUMBER init_q_rat_six = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_two = {
1,
1,
0,
{ 2,}
};
NUMBER init_q_rat_two = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_zero = {
1,
1,
0,
{ 0,}
};
NUMBER init_q_rat_zero = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_one = {
1,
1,
0,
{ 1,}
};
NUMBER init_q_rat_one = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_neg_one = {
-1,
1,
0,
{ 1,}
};
NUMBER init_q_rat_neg_one = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_half = {
1,
1,
0,
{ 1,}
};
NUMBER init_q_rat_half = {
1,
1,
0,
{ 2,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_ten = {
1,
1,
0,
{ 10,}
};
NUMBER init_q_rat_ten = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_pi = {
1,
6,
0,
{ 125527896, 283898350, 1960493936, 1672850762, 1288168272, 8,}
};
NUMBER init_q_pi = {
1,
6,
0,
{ 1288380402, 1120116153, 1860424692, 1944118326, 1583591604, 2,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_two_pi = {
1,
6,
0,
{ 251055792, 567796700, 1773504224, 1198217877, 428852897, 17,}
};
NUMBER init_q_two_pi = {
1,
6,
0,
{ 1288380402, 1120116153, 1860424692, 1944118326, 1583591604, 2,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_pi_over_two = {
1,
6,
0,
{ 125527896, 283898350, 1960493936, 1672850762, 1288168272, 8,}
};
NUMBER init_q_pi_over_two = {
1,
6,
0,
{ 429277156, 92748659, 1573365737, 1740753005, 1019699561, 5,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_one_pt_five_pi = {
1,
6,
0,
{ 1241201312, 270061909, 1051574664, 1924965045, 1340320627, 70,}
};
NUMBER init_q_one_pt_five_pi = {
1,
6,
0,
{ 1579671539, 1837970263, 1067644340, 523549916, 2119366659, 14,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_e_to_one_half = {
1,
6,
0,
{ 256945612, 216219427, 223516738, 477442596, 581063757, 23,}
};
NUMBER init_q_e_to_one_half = {
1,
6,
0,
{ 1536828363, 698484484, 1127331835, 224219346, 245499408, 14,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_exp = {
1,
6,
0,
{ 943665199, 1606559160, 1094967530, 1759391384, 1671799163, 1123581,}
};
NUMBER init_q_rat_exp = {
1,
6,
0,
{ 879242208, 2022880100, 617392930, 1374929092, 1367479163, 413342,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_ln_ten = {
1,
6,
0,
{ 2086268922, 165794492, 1416063951, 1851428830, 1893239400, 65366841,}
};
NUMBER init_q_ln_ten = {
1,
6,
0,
{ 26790652, 564532679, 783998273, 216030448, 1564709968, 28388458,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_ln_two = {
1,
6,
0,
{ 1789230241, 1057927868, 715399197, 908801241, 1411265331, 3,}
};
NUMBER init_q_ln_two = {
1,
6,
0,
{ 1559869847, 1930657510, 1228561531, 219003871, 593099283, 5,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rad_to_deg = {
1,
6,
0,
{ 2127722024, 1904928383, 2016479213, 2048947859, 1578647346, 492,}
};
NUMBER init_q_rad_to_deg = {
1,
6,
0,
{ 125527896, 283898350, 1960493936, 1672850762, 1288168272, 8,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rad_to_grad = {
1,
6,
0,
{ 2125526288, 684931327, 570267400, 129125085, 1038224725, 547,}
};
NUMBER init_q_rad_to_grad = {
1,
6,
0,
{ 125527896, 283898350, 1960493936, 1672850762, 1288168272, 8,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_qword = {
1,
3,
0,
{ 2147483647, 2147483647, 3,}
};
NUMBER init_q_rat_qword = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_dword = {
1,
2,
0,
{ 2147483647, 1,}
};
NUMBER init_q_rat_dword = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_max_long = {
1,
1,
0,
{ 2147483647,}
};
NUMBER init_q_rat_max_long = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_min_long = {
-1,
2,
0,
{ 0, 1,}
};
NUMBER init_q_rat_min_long = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_word = {
1,
1,
0,
{ 65535,}
};
NUMBER init_q_rat_word = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_byte = {
1,
1,
0,
{ 255,}
};
NUMBER init_q_rat_byte = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_400 = {
1,
1,
0,
{ 400,}
};
NUMBER init_q_rat_400 = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_360 = {
1,
1,
0,
{ 360,}
};
NUMBER init_q_rat_360 = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_200 = {
1,
1,
0,
{ 200,}
};
NUMBER init_q_rat_200 = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_180 = {
1,
1,
0,
{ 180,}
};
NUMBER init_q_rat_180 = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_max_exp = {
1,
1,
0,
{ 100000,}
};
NUMBER init_q_rat_max_exp = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_min_exp = {
-1,
1,
0,
{ 100000,}
};
NUMBER init_q_rat_min_exp = {
1,
1,
0,
{ 1,}
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_max_fact = {
1,
1,
0,
{ 3249, }
};
NUMBER init_q_rat_max_fact = {
1,
1,
0,
{ 1, }
};
// Autogenerated by _dumprawrat in support.c
NUMBER init_p_rat_min_fact = {
-1,
1,
0,
{ 1000, }
};
NUMBER init_q_rat_min_fact = {
1,
1,
0,
{ 1, }
};

View File

@@ -0,0 +1,452 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
#pragma once
//-----------------------------------------------------------------------------
// Package Title ratpak
// File ratpak.h
// Copyright (C) 1995-99 Microsoft
// Date 01-16-95
//
//
// Description
//
// Infinite precision math package header file, if you use ratpak.lib you
// need to include this header.
//
//-----------------------------------------------------------------------------
#include "CalcErr.h"
static constexpr uint32_t BASEXPWR = 31L;// Internal log2(BASEX)
static constexpr uint32_t BASEX = 0x80000000; // Internal radix used in calculations, hope to raise
// this to 2^32 after solving scaling problems with
// overflow detection esp. in mul
typedef unsigned long MANTTYPE;
typedef unsigned __int64 TWO_MANTTYPE;
enum eNUMOBJ_FMT {
FMT_FLOAT, // returns floating point, or exponential if number is too big
FMT_SCIENTIFIC, // always returns scientific notation
FMT_ENGINEERING // always returns engineering notation such that exponent is a multiple of 3
};
enum eANGLE_TYPE {
ANGLE_DEG, // Calculate trig using 360 degrees per revolution
ANGLE_RAD, // Calculate trig using 2 pi radians per revolution
ANGLE_GRAD // Calculate trig using 400 gradients per revolution
};
typedef enum eNUMOBJ_FMT NUMOBJ_FMT;
typedef enum eANGLE_TYPE ANGLE_TYPE;
//-----------------------------------------------------------------------------
//
// NUMBER type is a representation of a generic sized generic radix number
//
//-----------------------------------------------------------------------------
#pragma warning(push)
#pragma warning(disable:4200) // nonstandard extension used : zero-sized array in struct/union
typedef struct _number
{
long sign; // The sign of the mantissa, +1, or -1
long cdigit; // The number of digits, or what passes for digits in the
// radix being used.
long exp; // The offset of digits from the radix point
// (decimal point in radix 10)
MANTTYPE mant[];
// This is actually allocated as a continuation of the
// NUMBER structure.
} NUMBER, *PNUMBER, **PPNUMBER;
#pragma warning(pop)
//-----------------------------------------------------------------------------
//
// RAT type is a representation radix on 2 NUMBER types.
// pp/pq, where pp and pq are pointers to integral NUMBER types.
//
//-----------------------------------------------------------------------------
typedef struct _rat
{
PNUMBER pp;
PNUMBER pq;
} RAT, *PRAT;
static constexpr uint32_t MAX_LONG_SIZE = 33; // Base 2 requires 32 'digits'
//-----------------------------------------------------------------------------
//
// List of useful constants for evaluation, note this list needs to be
// initialized.
//
//-----------------------------------------------------------------------------
extern PNUMBER num_one;
extern PNUMBER num_two;
extern PNUMBER num_five;
extern PNUMBER num_six;
extern PNUMBER num_ten;
extern PRAT ln_ten;
extern PRAT ln_two;
extern PRAT rat_zero;
extern PRAT rat_neg_one;
extern PRAT rat_one;
extern PRAT rat_two;
extern PRAT rat_six;
extern PRAT rat_half;
extern PRAT rat_ten;
extern PRAT pt_eight_five;
extern PRAT pi;
extern PRAT pi_over_two;
extern PRAT two_pi;
extern PRAT one_pt_five_pi;
extern PRAT e_to_one_half;
extern PRAT rat_exp;
extern PRAT rad_to_deg;
extern PRAT rad_to_grad;
extern PRAT rat_qword;
extern PRAT rat_dword;
extern PRAT rat_word;
extern PRAT rat_byte;
extern PRAT rat_360;
extern PRAT rat_400;
extern PRAT rat_180;
extern PRAT rat_200;
extern PRAT rat_nRadix;
extern PRAT rat_smallest;
extern PRAT rat_negsmallest;
extern PRAT rat_max_exp;
extern PRAT rat_min_exp;
extern PRAT rat_max_fact;
extern PRAT rat_min_fact;
extern PRAT rat_max_long;
extern PRAT rat_min_long;
// DUPNUM Duplicates a number taking care of allocation and internals
#define DUPNUM(a,b) destroynum(a);createnum( a, (b)->cdigit );_dupnum(a, b);
// DUPRAT Duplicates a rational taking care of allocation and internals
#define DUPRAT(a,b) destroyrat(a);createrat(a);DUPNUM((a)->pp,(b)->pp);DUPNUM((a)->pq,(b)->pq);
// LOG*RADIX calculates the integral portion of the log of a number in
// the base currently being used, only accurate to within g_ratio
#define LOGNUMRADIX(pnum) (((pnum)->cdigit+(pnum)->exp)*g_ratio)
#define LOGRATRADIX(prat) (LOGNUMRADIX((prat)->pp)-LOGNUMRADIX((prat)->pq))
// LOG*2 calculates the integral portion of the log of a number in
// the internal base being used, only accurate to within g_ratio
#define LOGNUM2(pnum) ((pnum)->cdigit+(pnum)->exp)
#define LOGRAT2(prat) (LOGNUM2((prat)->pp)-LOGNUM2((prat)->pq))
#if defined( DEBUG_RATPAK )
//-----------------------------------------------------------------------------
//
// Debug versions of rational number creation and destruction routines.
// used for debugging allocation errors.
//
//-----------------------------------------------------------------------------
#define createrat(y) (y)=_createrat(); \
{ \
std::wstringstream outputString; \
outputString << "createrat " << y << " " << # y << " file= " << __FILE__ << ", line= " << __LINE__ << "\n"; \
OutputDebugString(outputString.str().c_str()); \
}
#define destroyrat(x) \
{ \
std::wstringstream outputString; \
outputString << "destroyrat " << x << " file= " << __FILE__ << ", line= " << __LINE__ << "\n"; \
OutputDebugString(outputString.str().c_str()); \
} \
_destroyrat(x),(x)=nullptr
#define createnum(y,x) (y)=_createnum(x); \
{ \
std::wstringstream outputString; \
outputString << "createnum " << y << " " << # y << " file= " << __FILE__ << ", line= " << __LINE__ << "\n"; \
OutputDebugString(outputString.str().c_str()); \
}
#define destroynum(x) \
{ \
std::wstringstream outputString; \
outputString << "destroynum " << x << " file= " << __FILE__ << ", line= " << __LINE__ << "\n"; \
OutputDebugString(outputString.str().c_str()); \
} \
_destroynum(x),(x)=nullptr
#else
#define createrat(y) (y)=_createrat()
#define destroyrat(x) _destroyrat(x),(x)=nullptr
#define createnum(y,x) (y)=_createnum(x)
#define destroynum(x) _destroynum(x),(x)=nullptr
#endif
//-----------------------------------------------------------------------------
//
// Defines for checking when to stop taylor series expansions due to
// precision satisfaction.
//
//-----------------------------------------------------------------------------
// RENORMALIZE, gets the exponents non-negative.
#define RENORMALIZE(x) if ( (x)->pp->exp < 0 ) { \
(x)->pq->exp -= (x)->pp->exp; \
(x)->pp->exp = 0; \
} \
if ( (x)->pq->exp < 0 ) { \
(x)->pp->exp -= (x)->pq->exp; \
(x)->pq->exp = 0; \
}
// TRIMNUM ASSUMES the number is in radix form NOT INTERNAL BASEX!!!
#define TRIMNUM(x, precision) if ( !g_ftrueinfinite ) { \
long trim = (x)->cdigit - precision-g_ratio;\
if ( trim > 1 ) \
{ \
memmove( (x)->mant, &((x)->mant[trim]), sizeof(MANTTYPE)*((x)->cdigit-trim) ); \
(x)->cdigit -= trim; \
(x)->exp += trim; \
} \
}
// TRIMTOP ASSUMES the number is in INTERNAL BASEX!!!
#define TRIMTOP(x, precision) if ( !g_ftrueinfinite ) { \
long trim = (x)->pp->cdigit - (precision/g_ratio) - 2;\
if ( trim > 1 ) \
{ \
memmove( (x)->pp->mant, &((x)->pp->mant[trim]), sizeof(MANTTYPE)*((x)->pp->cdigit-trim) ); \
(x)->pp->cdigit -= trim; \
(x)->pp->exp += trim; \
} \
trim = min((x)->pp->exp,(x)->pq->exp);\
(x)->pp->exp -= trim;\
(x)->pq->exp -= trim;\
}
#define SMALL_ENOUGH_RAT(a, precision) (zernum((a)->pp) || ( ( ( (a)->pq->cdigit + (a)->pq->exp ) - ( (a)->pp->cdigit + (a)->pp->exp ) - 1 ) * g_ratio > precision ) )
//-----------------------------------------------------------------------------
//
// Defines for setting up taylor series expansions for infinite precision
// functions.
//
//-----------------------------------------------------------------------------
#define CREATETAYLOR() PRAT xx=nullptr;\
PNUMBER n2=nullptr; \
PRAT pret=nullptr; \
PRAT thisterm=nullptr; \
DUPRAT(xx,*px); \
mulrat(&xx,*px, precision); \
createrat(pret); \
pret->pp=longtonum( 0L, BASEX ); \
pret->pq=longtonum( 0L, BASEX );
#define DESTROYTAYLOR() destroynum( n2 ); \
destroyrat( xx );\
destroyrat( thisterm );\
destroyrat( *px );\
trimit(&pret, precision);\
*px=pret;
// INC(a) is the rational equivalent of a++
// Check to see if we can avoid doing this the hard way.
#define INC(a) if ( (a)->mant[0] < BASEX - 1 ) \
{ \
(a)->mant[0]++; \
} \
else \
{ \
addnum( &(a), num_one, BASEX); \
}
#define MSD(x) ((x)->mant[(x)->cdigit-1])
// MULNUM(b) is the rational equivalent of thisterm *= b where thisterm is
// a rational and b is a number, NOTE this is a mixed type operation for
// efficiency reasons.
#define MULNUM(b) mulnumx( &(thisterm->pp), b);
// DIVNUM(b) is the rational equivalent of thisterm /= b where thisterm is
// a rational and b is a number, NOTE this is a mixed type operation for
// efficiency reasons.
#define DIVNUM(b) mulnumx( &(thisterm->pq), b);
// NEXTTERM(p,d) is the rational equivalent of
// thisterm *= p
// d <d is usually an expansion of operations to get thisterm updated.>
// pret += thisterm
#define NEXTTERM(p,d,precision) mulrat(&thisterm,p,precision);d addrat( &pret, thisterm, precision )
//-----------------------------------------------------------------------------
//
// External variables used in the math package.
//
//-----------------------------------------------------------------------------
extern bool g_ftrueinfinite; // set to true to allow infinite precision
// don't use unless you know what you are doing
// used to help decide when to stop calculating.
extern long g_ratio; // Internally calculated ratio of internal radix
//-----------------------------------------------------------------------------
//
// External functions defined in the math package.
//
//-----------------------------------------------------------------------------
// Call whenever decimal separator character changes.
extern void SetDecimalSeparator(wchar_t decimalSeparator);
// Call whenever either radix or precision changes, is smarter about recalculating constants.
extern void ChangeConstants(uint32_t radix, int32_t precision);
extern bool equnum(_In_ PNUMBER a, _In_ PNUMBER b ); // returns true of a == b
extern bool lessnum(_In_ PNUMBER a, _In_ PNUMBER b ); // returns true of a < b
extern bool zernum(_In_ PNUMBER a ); // returns true of a == 0
extern bool zerrat(_In_ PRAT a ); // returns true if a == 0/q
extern std::wstring NumberToString(_Inout_ PNUMBER& pnum, int format, uint32_t radix, int32_t precision);
// returns a text representation of a PRAT
extern std::wstring RatToString(_Inout_ PRAT& prat, int format, uint32_t radix, int32_t precision);
extern long numtolong(_In_ PNUMBER pnum, uint32_t radix );
extern long rattolong(_In_ PRAT prat, uint32_t radix, int32_t precision);
ULONGLONG rattoUlonglong(_In_ PRAT prat, uint32_t radix, int32_t precision);
extern PNUMBER _createnum(_In_ ULONG size ); // returns an empty number structure with size digits
extern PNUMBER nRadixxtonum(_In_ PNUMBER a, uint32_t radix, int32_t precision);
extern PNUMBER gcd(_In_ PNUMBER a, _In_ PNUMBER b );
extern PNUMBER StringToNumber(std::wstring_view numberString, uint32_t radix, int32_t precision); // takes a text representation of a number and returns a number.
// takes a text representation of a number as a mantissa with sign and an exponent with sign.
extern PRAT StringToRat(bool mantissaIsNegative, std::wstring_view mantissa, bool exponentIsNegative, std::wstring_view exponent, uint32_t radix, int32_t precision);
extern PNUMBER longfactnum(long inlong, uint32_t radix);
extern PNUMBER longprodnum(long start, long stop, uint32_t radix);
extern PNUMBER longtonum(long inlong, uint32_t radix);
extern PNUMBER Ulongtonum(unsigned long inlong, uint32_t radix);
extern PNUMBER numtonRadixx(PNUMBER a, uint32_t radix);
// creates a empty/undefined rational representation (p/q)
extern PRAT _createrat( void );
// returns a new rat structure with the acos of x->p/x->q taking into account
// angle type
extern void acosanglerat( _Inout_ PRAT *px, ANGLE_TYPE angletype, uint32_t radix, int32_t precision);
// returns a new rat structure with the acosh of x->p/x->q
extern void acoshrat( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
// returns a new rat structure with the acos of x->p/x->q
extern void acosrat( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
// returns a new rat structure with the asin of x->p/x->q taking into account
// angle type
extern void asinanglerat( _Inout_ PRAT *px, ANGLE_TYPE angletype, uint32_t radix, int32_t precision);
extern void asinhrat( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
// returns a new rat structure with the asinh of x->p/x->q
// returns a new rat structure with the asin of x->p/x->q
extern void asinrat( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
// returns a new rat structure with the atan of x->p/x->q taking into account
// angle type
extern void atananglerat( _Inout_ PRAT *px, ANGLE_TYPE angletype, uint32_t radix, int32_t precision);
// returns a new rat structure with the atanh of x->p/x->q
extern void atanhrat( _Inout_ PRAT *px, int32_t precision);
// returns a new rat structure with the atan of x->p/x->q
extern void atanrat( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
// returns a new rat structure with the cosh of x->p/x->q
extern void coshrat( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
// returns a new rat structure with the cos of x->p/x->q
extern void cosrat( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
// returns a new rat structure with the cos of x->p/x->q taking into account
// angle type
extern void cosanglerat( _Inout_ PRAT *px, ANGLE_TYPE angletype, uint32_t radix, int32_t precision);
// returns a new rat structure with the exp of x->p/x->q this should not be called explicitly.
extern void _exprat( _Inout_ PRAT *px, int32_t precision);
// returns a new rat structure with the exp of x->p/x->q
extern void exprat( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
// returns a new rat structure with the log base 10 of x->p/x->q
extern void log10rat( _Inout_ PRAT *px, int32_t precision);
// returns a new rat structure with the natural log of x->p/x->q
extern void lograt( _Inout_ PRAT *px, int32_t precision);
extern PRAT longtorat( long inlong );
extern PRAT Ulongtorat( unsigned long inulong );
extern PRAT numtorat( _In_ PNUMBER pin, uint32_t radix);
extern void sinhrat( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
extern void sinrat( _Inout_ PRAT *px );
// returns a new rat structure with the sin of x->p/x->q taking into account
// angle type
extern void sinanglerat( _Inout_ PRAT *px, ANGLE_TYPE angletype, uint32_t radix, int32_t precision);
extern void tanhrat( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
extern void tanrat( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
// returns a new rat structure with the tan of x->p/x->q taking into account
// angle type
extern void tananglerat( _Inout_ PRAT *px, ANGLE_TYPE angletype, uint32_t radix, int32_t precision);
extern void _dupnum(_In_ PNUMBER dest, _In_ PNUMBER src);
extern void _destroynum( _In_ PNUMBER pnum );
extern void _destroyrat( _In_ PRAT prat );
extern void addnum( _Inout_ PNUMBER *pa, _In_ PNUMBER b, uint32_t radix);
extern void addrat( _Inout_ PRAT *pa, _In_ PRAT b, int32_t precision);
extern void andrat( _Inout_ PRAT *pa, _In_ PRAT b, uint32_t radix, int32_t precision);
extern void divnum( _Inout_ PNUMBER *pa, _In_ PNUMBER b, uint32_t radix, int32_t precision);
extern void divnumx( _Inout_ PNUMBER *pa, _In_ PNUMBER b, int32_t precision);
extern void divrat( _Inout_ PRAT *pa, _In_ PRAT b, int32_t precision);
extern void fracrat( _Inout_ PRAT *pa , uint32_t radix, int32_t precision);
extern void factrat( _Inout_ PRAT *pa, uint32_t radix, int32_t precision);
extern void modrat( _Inout_ PRAT *pa, _In_ PRAT b );
extern void gcdrat( _Inout_ PRAT *pa, uint32_t radix, int32_t precision);
extern void intrat( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
extern void mulnum( _Inout_ PNUMBER *pa, _In_ PNUMBER b, uint32_t radix);
extern void mulnumx( _Inout_ PNUMBER *pa, _In_ PNUMBER b );
extern void mulrat( _Inout_ PRAT *pa, _In_ PRAT b, int32_t precision);
extern void numpowlong( _Inout_ PNUMBER *proot, long power, uint32_t radix, int32_t precision);
extern void numpowlongx( _Inout_ PNUMBER *proot, long power );
extern void orrat( _Inout_ PRAT *pa, _In_ PRAT b, uint32_t radix, int32_t precision);
extern void powrat( _Inout_ PRAT *pa, _In_ PRAT b , uint32_t radix, int32_t precision);
extern void powratNumeratorDenominator(_Inout_ PRAT *pa, _In_ PRAT b, uint32_t radix, int32_t precision);
extern void powratcomp(_Inout_ PRAT *pa, _In_ PRAT b, uint32_t radix, int32_t precision);
extern void ratpowlong( _Inout_ PRAT *proot, long power, int32_t precision);
extern void remnum( _Inout_ PNUMBER *pa, _In_ PNUMBER b, uint32_t radix);
extern void rootrat( _Inout_ PRAT *pa, _In_ PRAT b , uint32_t radix, int32_t precision);
extern void scale2pi( _Inout_ PRAT *px, uint32_t radix, int32_t precision);
extern void scale( _Inout_ PRAT *px, _In_ PRAT scalefact, uint32_t radix, int32_t precision);
extern void subrat( _Inout_ PRAT *pa, _In_ PRAT b, int32_t precision);
extern void xorrat( _Inout_ PRAT *pa, _In_ PRAT b, uint32_t radix, int32_t precision);
extern void lshrat( _Inout_ PRAT *pa, _In_ PRAT b , uint32_t radix, int32_t precision);
extern void rshrat( _Inout_ PRAT *pa, _In_ PRAT b, uint32_t radix, int32_t precision);
extern bool rat_equ( _In_ PRAT a, _In_ PRAT b, int32_t precision);
extern bool rat_neq( _In_ PRAT a, _In_ PRAT b, int32_t precision);
extern bool rat_gt( _In_ PRAT a, _In_ PRAT b, int32_t precision);
extern bool rat_ge( _In_ PRAT a, _In_ PRAT b, int32_t precision);
extern bool rat_lt( _In_ PRAT a, _In_ PRAT b, int32_t precision);
extern bool rat_le( _In_ PRAT a, _In_ PRAT b, int32_t precision);
extern void inbetween( _In_ PRAT *px, _In_ PRAT range, int32_t precision);
extern void trimit( _Inout_ PRAT *px, int32_t precision);
extern void _dumprawrat(_In_ const wchar_t *varname, _In_ PRAT rat, std::wostream& out);
extern void _dumprawnum(_In_ const wchar_t *varname, _In_ PNUMBER num, std::wostream& out);

View File

@@ -0,0 +1,719 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
//----------------------------------------------------------------------------
// Package Title ratpak
// File support.c
// Copyright (C) 1995-96 Microsoft
// Date 10-21-96
//
//
// Description
//
// Contains support functions for rationals and numbers.
//
// Special Information
//
//
//
//----------------------------------------------------------------------------
#include "pch.h"
#include "ratpak.h"
using namespace std;
void _readconstants( void );
#if defined( GEN_CONST )
static int cbitsofprecision = 0;
#define READRAWRAT(v)
#define READRAWNUM(v)
#define DUMPRAWRAT(v) _dumprawrat(#v,v, wcout)
#define DUMPRAWNUM(v) fprintf( stderr, \
"// Autogenerated by _dumprawrat in support.c\n" ); \
fprintf( stderr, "NUMBER init_" #v "= {\n" ); \
_dumprawnum(v, wcout); \
fprintf( stderr, "};\n" )
#else
#define DUMPRAWRAT(v)
#define DUMPRAWNUM(v)
#define READRAWRAT(v) createrat(v); DUPNUM((v)->pp,(&(init_p_##v))); \
DUPNUM((v)->pq,(&(init_q_##v)));
#define READRAWNUM(v) DUPNUM(v,(&(init_##v)))
#define INIT_AND_DUMP_RAW_NUM_IF_NULL(r, v) if (r == nullptr) { r = longtonum(v, BASEX); DUMPRAWNUM(v); }
#define INIT_AND_DUMP_RAW_RAT_IF_NULL(r, v) if (r == nullptr) { r = longtorat(v); DUMPRAWRAT(v); }
static constexpr int RATIO_FOR_DECIMAL = 9;
static constexpr int DECIMAL = 10;
static constexpr int CALC_DECIMAL_DIGITS_DEFAULT = 32;
static int cbitsofprecision = RATIO_FOR_DECIMAL * DECIMAL * CALC_DECIMAL_DIGITS_DEFAULT;
#include "ratconst.h"
#endif
bool g_ftrueinfinite = false; // Set to true if you don't want
// chopping internally
// precision used internally
PNUMBER num_one= nullptr;
PNUMBER num_two= nullptr;
PNUMBER num_five= nullptr;
PNUMBER num_six= nullptr;
PNUMBER num_ten= nullptr;
PRAT ln_ten= nullptr;
PRAT ln_two= nullptr;
PRAT rat_zero= nullptr;
PRAT rat_one= nullptr;
PRAT rat_neg_one= nullptr;
PRAT rat_two= nullptr;
PRAT rat_six= nullptr;
PRAT rat_half= nullptr;
PRAT rat_ten= nullptr;
PRAT pt_eight_five= nullptr;
PRAT pi= nullptr;
PRAT pi_over_two= nullptr;
PRAT two_pi= nullptr;
PRAT one_pt_five_pi= nullptr;
PRAT e_to_one_half= nullptr;
PRAT rat_exp= nullptr;
PRAT rad_to_deg= nullptr;
PRAT rad_to_grad= nullptr;
PRAT rat_qword= nullptr;
PRAT rat_dword= nullptr; // unsigned max ulong
PRAT rat_word= nullptr;
PRAT rat_byte= nullptr;
PRAT rat_360= nullptr;
PRAT rat_400= nullptr;
PRAT rat_180= nullptr;
PRAT rat_200= nullptr;
PRAT rat_nRadix= nullptr;
PRAT rat_smallest= nullptr;
PRAT rat_negsmallest= nullptr;
PRAT rat_max_exp= nullptr;
PRAT rat_min_exp= nullptr;
PRAT rat_max_fact = nullptr;
PRAT rat_min_fact = nullptr;
PRAT rat_min_long= nullptr; // min signed long
PRAT rat_max_long= nullptr; // max signed long
//----------------------------------------------------------------------------
//
// FUNCTION: ChangeConstants
//
// ARGUMENTS: base changing to, and precision to use.
//
// RETURN: None
//
// SIDE EFFECTS: sets a mess of constants.
//
//
//----------------------------------------------------------------------------
void ChangeConstants(uint32_t radix, int32_t precision)
{
// ratio is set to the number of digits in the current radix, you can get
// in the internal BASEX radix, this is important for length calculations
// in translating from radix to BASEX and back.
uint64_t limit = static_cast<uint64_t>(BASEX) / static_cast<uint64_t>(radix);
g_ratio = 0;
for (uint32_t digit = 1; digit < limit; digit *= radix )
{
g_ratio++;
}
g_ratio += !g_ratio;
destroyrat(rat_nRadix);
rat_nRadix=longtorat( radix );
// Check to see what we have to recalculate and what we don't
if (cbitsofprecision < (g_ratio * static_cast<int32_t>(radix) * precision))
{
g_ftrueinfinite = false;
INIT_AND_DUMP_RAW_NUM_IF_NULL(num_one, 1L);
INIT_AND_DUMP_RAW_NUM_IF_NULL(num_two, 2L);
INIT_AND_DUMP_RAW_NUM_IF_NULL(num_five, 5L);
INIT_AND_DUMP_RAW_NUM_IF_NULL(num_six, 6L);
INIT_AND_DUMP_RAW_NUM_IF_NULL(num_ten, 10L);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_six, 6L);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_two, 2L);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_zero, 0L);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_one, 1L);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_neg_one, -1L);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_ten, 10L);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_word, 0xffff);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_word, 0xff);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_400, 400);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_360, 360);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_200, 200);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_180, 180);
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_max_exp, 100000);
// 3248, is the max number for which calc is able to compute factorial, after that it is unable to compute due to overflow.
// Hence restricted factorial range as at most 3248.Beyond that calc will throw overflow error immediately.
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_max_fact, 3249);
// -1000, is the min number for which calc is able to compute factorial, after that it takes too long to compute.
INIT_AND_DUMP_RAW_RAT_IF_NULL(rat_min_fact, -1000);
DUPRAT(rat_smallest, rat_nRadix);
ratpowlong(&rat_smallest, -precision, precision);
DUPRAT(rat_negsmallest, rat_smallest);
rat_negsmallest->pp->sign = -1;
DUMPRAWRAT(rat_smallest);
DUMPRAWRAT(rat_negsmallest);
if (rat_half == nullptr)
{
createrat(rat_half);
DUPNUM(rat_half->pp, num_one);
DUPNUM(rat_half->pq, num_two);
DUMPRAWRAT(rat_half);
}
if (pt_eight_five == nullptr)
{
createrat(pt_eight_five);
pt_eight_five->pp = longtonum(85L, BASEX);
pt_eight_five->pq = longtonum(100L, BASEX);
DUMPRAWRAT(pt_eight_five);
}
DUPRAT(rat_qword, rat_two);
numpowlong(&(rat_qword->pp), 64, BASEX, precision);
subrat(&rat_qword, rat_one, precision);
DUMPRAWRAT(rat_qword);
DUPRAT(rat_dword, rat_two);
numpowlong(&(rat_dword->pp), 32, BASEX, precision);
subrat(&rat_dword, rat_one, precision);
DUMPRAWRAT(rat_dword);
DUPRAT(rat_max_long, rat_two);
numpowlong(&(rat_max_long->pp), 31, BASEX, precision);
DUPRAT(rat_min_long, rat_max_long);
subrat(&rat_max_long, rat_one, precision); // rat_max_long = 2^31 -1
DUMPRAWRAT(rat_max_long);
rat_min_long->pp->sign *= -1; // rat_min_long = -2^31
DUMPRAWRAT(rat_min_long);
DUPRAT(rat_min_exp, rat_max_exp);
rat_min_exp->pp->sign *= -1;
DUMPRAWRAT(rat_min_exp);
cbitsofprecision = g_ratio * radix * precision;
// Apparently when dividing 180 by pi, another (internal) digit of
// precision is needed.
long extraPrecision = precision + g_ratio;
DUPRAT(pi, rat_half);
asinrat(&pi, radix, extraPrecision);
mulrat(&pi, rat_six, extraPrecision);
DUMPRAWRAT(pi);
DUPRAT(two_pi, pi);
DUPRAT(pi_over_two, pi);
DUPRAT(one_pt_five_pi, pi);
addrat(&two_pi, pi, extraPrecision);
DUMPRAWRAT(two_pi);
divrat(&pi_over_two, rat_two, extraPrecision);
DUMPRAWRAT(pi_over_two);
addrat(&one_pt_five_pi, pi_over_two, extraPrecision);
DUMPRAWRAT(one_pt_five_pi);
DUPRAT(e_to_one_half, rat_half);
_exprat(&e_to_one_half, extraPrecision);
DUMPRAWRAT(e_to_one_half);
DUPRAT(rat_exp, rat_one);
_exprat(&rat_exp, extraPrecision);
DUMPRAWRAT(rat_exp);
// WARNING: remember lograt uses exponent constants calculated above...
DUPRAT(ln_ten, rat_ten);
lograt(&ln_ten, extraPrecision);
DUMPRAWRAT(ln_ten);
DUPRAT(ln_two, rat_two);
lograt(&ln_two, extraPrecision);
DUMPRAWRAT(ln_two);
destroyrat(rad_to_deg);
rad_to_deg = longtorat(180L);
divrat(&rad_to_deg, pi, extraPrecision);
DUMPRAWRAT(rad_to_deg);
destroyrat(rad_to_grad);
rad_to_grad = longtorat(200L);
divrat(&rad_to_grad, pi, extraPrecision);
DUMPRAWRAT(rad_to_grad);
}
else
{
_readconstants();
DUPRAT(rat_smallest, rat_nRadix);
ratpowlong(&rat_smallest, -precision, precision);
DUPRAT(rat_negsmallest, rat_smallest);
rat_negsmallest->pp->sign = -1;
}
}
//----------------------------------------------------------------------------
//
// FUNCTION: intrat
//
// ARGUMENTS: pointer to x PRAT representation of number
//
// RETURN: no return value x PRAT is smashed with integral number
//
//
//----------------------------------------------------------------------------
void intrat( PRAT *px, uint32_t radix, int32_t precision)
{
// Only do the intrat operation if number is nonzero.
// and only if the bottom part is not one.
if ( !zernum( (*px)->pp ) && !equnum( (*px)->pq, num_one ) )
{
wstring ratStr = RatToString(*px, FMT_FLOAT, radix, precision);
PNUMBER pnum = StringToNumber(ratStr, radix, precision);
destroyrat( *px );
*px = numtorat( pnum, radix);
destroynum( pnum );
PRAT pret = nullptr;
DUPRAT(pret,*px);
modrat( &pret, rat_one );
subrat( px, pret, precision);
destroyrat( pret );
}
}
//---------------------------------------------------------------------------
//
// FUNCTION: rat_equ
//
// ARGUMENTS: PRAT a and PRAT b
//
// RETURN: true if equal false otherwise.
//
//
//---------------------------------------------------------------------------
bool rat_equ( PRAT a, PRAT b, int32_t precision)
{
PRAT rattmp= nullptr;
DUPRAT(rattmp,a);
rattmp->pp->sign *= -1;
addrat( &rattmp, b, precision);
bool bret = zernum( rattmp->pp );
destroyrat( rattmp );
return( bret );
}
//---------------------------------------------------------------------------
//
// FUNCTION: rat_ge
//
// ARGUMENTS: PRAT a, PRAT b and long precision
//
// RETURN: true if a is greater than or equal to b
//
//
//---------------------------------------------------------------------------
bool rat_ge( PRAT a, PRAT b, int32_t precision)
{
PRAT rattmp= nullptr;
DUPRAT(rattmp,a);
b->pp->sign *= -1;
addrat( &rattmp, b, precision);
b->pp->sign *= -1;
bool bret = ( zernum( rattmp->pp ) ||
rattmp->pp->sign * rattmp->pq->sign == 1 );
destroyrat( rattmp );
return( bret );
}
//---------------------------------------------------------------------------
//
// FUNCTION: rat_gt
//
// ARGUMENTS: PRAT a and PRAT b
//
// RETURN: true if a is greater than b
//
//
//---------------------------------------------------------------------------
bool rat_gt( PRAT a, PRAT b, int32_t precision)
{
PRAT rattmp= nullptr;
DUPRAT(rattmp,a);
b->pp->sign *= -1;
addrat( &rattmp, b, precision);
b->pp->sign *= -1;
bool bret = ( !zernum( rattmp->pp ) &&
rattmp->pp->sign * rattmp->pq->sign == 1 );
destroyrat( rattmp );
return( bret );
}
//---------------------------------------------------------------------------
//
// FUNCTION: rat_le
//
// ARGUMENTS: PRAT a, PRAT b and long precision
//
// RETURN: true if a is less than or equal to b
//
//
//---------------------------------------------------------------------------
bool rat_le( PRAT a, PRAT b, int32_t precision)
{
PRAT rattmp= nullptr;
DUPRAT(rattmp,a);
b->pp->sign *= -1;
addrat( &rattmp, b, precision);
b->pp->sign *= -1;
bool bret = ( zernum( rattmp->pp ) ||
rattmp->pp->sign * rattmp->pq->sign == -1 );
destroyrat( rattmp );
return( bret );
}
//---------------------------------------------------------------------------
//
// FUNCTION: rat_lt
//
// ARGUMENTS: PRAT a, PRAT b and long precision
//
// RETURN: true if a is less than b
//
//
//---------------------------------------------------------------------------
bool rat_lt( PRAT a, PRAT b, int32_t precision)
{
PRAT rattmp= nullptr;
DUPRAT(rattmp,a);
b->pp->sign *= -1;
addrat( &rattmp, b, precision);
b->pp->sign *= -1;
bool bret = ( !zernum( rattmp->pp ) &&
rattmp->pp->sign * rattmp->pq->sign == -1 );
destroyrat( rattmp );
return( bret );
}
//---------------------------------------------------------------------------
//
// FUNCTION: rat_neq
//
// ARGUMENTS: PRAT a and PRAT b
//
// RETURN: true if a is not equal to b
//
//
//---------------------------------------------------------------------------
bool rat_neq( PRAT a, PRAT b, int32_t precision)
{
PRAT rattmp= nullptr;
DUPRAT(rattmp,a);
rattmp->pp->sign *= -1;
addrat( &rattmp, b, precision);
bool bret = !( zernum( rattmp->pp ) );
destroyrat( rattmp );
return( bret );
}
//---------------------------------------------------------------------------
//
// function: scale
//
// ARGUMENTS: pointer to x PRAT representation of number, and scaling factor
//
// RETURN: no return, value x PRAT is smashed with a scaled number in the
// range of the scalefact.
//
//---------------------------------------------------------------------------
void scale( PRAT *px, PRAT scalefact, uint32_t radix, int32_t precision )
{
PRAT pret = nullptr;
DUPRAT(pret,*px);
// Logscale is a quick way to tell how much extra precision is needed for
// scaleing by scalefact.
long logscale = g_ratio * ( (pret->pp->cdigit+pret->pp->exp) -
(pret->pq->cdigit+pret->pq->exp) );
if ( logscale > 0 )
{
precision += logscale;
}
divrat( &pret, scalefact, precision);
intrat(&pret, radix, precision);
mulrat( &pret, scalefact, precision);
pret->pp->sign *= -1;
addrat( px, pret, precision);
destroyrat( pret );
}
//---------------------------------------------------------------------------
//
// function: scale2pi
//
// ARGUMENTS: pointer to x PRAT representation of number
//
// RETURN: no return, value x PRAT is smashed with a scaled number in the
// range of 0..2pi
//
//---------------------------------------------------------------------------
void scale2pi( PRAT *px, uint32_t radix, int32_t precision )
{
PRAT pret = nullptr;
PRAT my_two_pi = nullptr;
DUPRAT(pret,*px);
// Logscale is a quick way to tell how much extra precision is needed for
// scaleing by 2 pi.
long logscale = g_ratio * ( (pret->pp->cdigit+pret->pp->exp) -
(pret->pq->cdigit+pret->pq->exp) );
if ( logscale > 0 )
{
precision += logscale;
DUPRAT(my_two_pi,rat_half);
asinrat( &my_two_pi, radix, precision);
mulrat( &my_two_pi, rat_six, precision);
mulrat( &my_two_pi, rat_two, precision);
}
else
{
DUPRAT(my_two_pi,two_pi);
logscale = 0;
}
divrat( &pret, my_two_pi, precision);
intrat(&pret, radix, precision);
mulrat( &pret, my_two_pi, precision);
pret->pp->sign *= -1;
addrat( px, pret, precision);
destroyrat( my_two_pi );
destroyrat( pret );
}
//---------------------------------------------------------------------------
//
// FUNCTION: inbetween
//
// ARGUMENTS: PRAT *px, and PRAT range.
//
// RETURN: none, changes *px to -/+range, if px is outside -range..+range
//
//---------------------------------------------------------------------------
void inbetween( PRAT *px, PRAT range, int32_t precision)
{
if ( rat_gt(*px,range, precision) )
{
DUPRAT(*px,range);
}
else
{
range->pp->sign *= -1;
if ( rat_lt(*px, range, precision) )
{
DUPRAT(*px,range);
}
range->pp->sign *= -1;
}
}
//---------------------------------------------------------------------------
//
// FUNCTION: _dumprawrat
//
// ARGUMENTS: const wchar *name of variable, PRAT x, output stream out
//
// RETURN: none, prints the results of a dump of the internal structures
// of a PRAT, suitable for READRAWRAT to stderr.
//
//---------------------------------------------------------------------------
void _dumprawrat( const wchar_t *varname, PRAT rat, wostream& out)
{
_dumprawnum(varname, rat->pp, out );
_dumprawnum(varname, rat->pq, out );
}
//---------------------------------------------------------------------------
//
// FUNCTION: _dumprawnum
//
// ARGUMENTS: const wchar *name of variable, PNUMBER num, output stream out
//
// RETURN: none, prints the results of a dump of the internal structures
// of a PNUMBER, suitable for READRAWNUM to stderr.
//
//---------------------------------------------------------------------------
void _dumprawnum(const wchar_t *varname, PNUMBER num, wostream& out)
{
int i;
out << L"NUMBER " << varname << L" = {\n";
out << L"\t"<< num->sign << L",\n";
out << L"\t" << num->cdigit << L",\n";
out << L"\t" << num->exp << L",\n";
out << L"\t{ ";
for ( i = 0; i < num->cdigit; i++ )
{
out << L" "<< num->mant[i] << L",";
}
out << L"}\n";
out << L"};\n";
}
void _readconstants( void )
{
READRAWNUM(num_one);
READRAWNUM(num_two);
READRAWNUM(num_five);
READRAWNUM(num_six);
READRAWNUM(num_ten);
READRAWRAT(pt_eight_five);
READRAWRAT(rat_six);
READRAWRAT(rat_two);
READRAWRAT(rat_zero);
READRAWRAT(rat_one);
READRAWRAT(rat_neg_one);
READRAWRAT(rat_half);
READRAWRAT(rat_ten);
READRAWRAT(pi);
READRAWRAT(two_pi);
READRAWRAT(pi_over_two);
READRAWRAT(one_pt_five_pi);
READRAWRAT(e_to_one_half);
READRAWRAT(rat_exp);
READRAWRAT(ln_ten);
READRAWRAT(ln_two);
READRAWRAT(rad_to_deg);
READRAWRAT(rad_to_grad);
READRAWRAT(rat_qword);
READRAWRAT(rat_dword);
READRAWRAT(rat_word);
READRAWRAT(rat_byte);
READRAWRAT(rat_360);
READRAWRAT(rat_400);
READRAWRAT(rat_180);
READRAWRAT(rat_200);
READRAWRAT(rat_smallest);
READRAWRAT(rat_negsmallest);
READRAWRAT(rat_max_exp);
READRAWRAT(rat_min_exp);
READRAWRAT(rat_max_fact);
READRAWRAT(rat_min_fact);
READRAWRAT(rat_min_long);
READRAWRAT(rat_max_long);
}
//---------------------------------------------------------------------------
//
// FUNCTION: trimit
//
// ARGUMENTS: PRAT *px, long precision
//
//
// DESCRIPTION: Chops off digits from rational numbers to avoid time
// explosions in calculations of functions using series.
// It can be shown that it is enough to only keep the first n digits
// of the largest of p or q in the rational p over q form, and of course
// scale the smaller by the same number of digits. This will give you
// n-1 digits of accuracy. This dramatically speeds up calculations
// involving hundreds of digits or more.
// The last part of this trim dealing with exponents never affects accuracy
//
// RETURN: none, modifies the pointed to PRAT
//
//---------------------------------------------------------------------------
void trimit( PRAT *px, int32_t precision)
{
if ( !g_ftrueinfinite )
{
long trim;
PNUMBER pp=(*px)->pp;
PNUMBER pq=(*px)->pq;
trim = g_ratio * (min((pp->cdigit+pp->exp),(pq->cdigit+pq->exp))-1) - precision;
if ( trim > g_ratio )
{
trim /= g_ratio;
if ( trim <= pp->exp )
{
pp->exp -= trim;
}
else
{
memmove( pp->mant, &(pp->mant[trim-pp->exp]), sizeof(MANTTYPE)*(pp->cdigit-trim+pp->exp) );
pp->cdigit -= trim-pp->exp;
pp->exp = 0;
}
if ( trim <= pq->exp )
{
pq->exp -= trim;
}
else
{
memmove( pq->mant, &(pq->mant[trim-pq->exp]), sizeof(MANTTYPE)*(pq->cdigit-trim+pq->exp) );
pq->cdigit -= trim-pq->exp;
pq->exp = 0;
}
}
trim = min(pp->exp,pq->exp);
pp->exp -= trim;
pq->exp -= trim;
}
}

View File

@@ -0,0 +1,297 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
//----------------------------------------------------------------------------
// File trans.c
// Copyright (C) 1995-96 Microsoft
// Date 01-16-95
//
//
// Description
//
// Contains sin, cos and tan for rationals
//
//
//----------------------------------------------------------------------------
#include "pch.h"
#include "ratpak.h"
void scalerat( _Inout_ PRAT *pa, ANGLE_TYPE angletype, uint32_t radix, int32_t precision )
{
switch ( angletype )
{
case ANGLE_RAD:
scale2pi( pa, radix, precision);
break;
case ANGLE_DEG:
scale( pa, rat_360, radix, precision);
break;
case ANGLE_GRAD:
scale( pa, rat_400, radix, precision);
break;
}
}
//-----------------------------------------------------------------------------
//
// FUNCTION: sinrat, _sinrat
//
// ARGUMENTS: x PRAT representation of number to take the sine of
//
// RETURN: sin of x in PRAT form.
//
// EXPLANATION: This uses Taylor series
//
// n
// ___ 2j+1
// \ ] j X
// \ -1 * ---------
// / (2j+1)!
// /__]
// j=0
// or,
// n
// ___ 2
// \ ] -X
// \ thisterm ; where thisterm = thisterm * ---------
// / j j+1 j (2j)*(2j+1)
// /__]
// j=0
//
// thisterm = X ; and stop when thisterm < precision used.
// 0 n
//
//-----------------------------------------------------------------------------
void _sinrat( PRAT *px, int32_t precision)
{
CREATETAYLOR();
DUPRAT(pret,*px);
DUPRAT(thisterm,*px);
DUPNUM(n2,num_one);
xx->pp->sign *= -1;
do {
NEXTTERM(xx,INC(n2) DIVNUM(n2) INC(n2) DIVNUM(n2), precision);
} while ( !SMALL_ENOUGH_RAT( thisterm, precision) );
DESTROYTAYLOR();
// Since *px might be epsilon above 1 or below -1, due to TRIMIT we need
// this trick here.
inbetween(px, rat_one, precision);
// Since *px might be epsilon near zero we must set it to zero.
if ( rat_le(*px, rat_smallest, precision) && rat_ge(*px, rat_negsmallest, precision) )
{
DUPRAT(*px,rat_zero);
}
}
void sinrat( PRAT *px, uint32_t radix, int32_t precision)
{
scale2pi(px, radix, precision);
_sinrat(px, precision);
}
void sinanglerat( _Inout_ PRAT *pa, ANGLE_TYPE angletype, uint32_t radix, int32_t precision)
{
scalerat( pa, angletype, radix, precision);
switch ( angletype )
{
case ANGLE_DEG:
if ( rat_gt( *pa, rat_180, precision) )
{
subrat(pa, rat_360, precision);
}
divrat( pa, rat_180, precision);
mulrat( pa, pi, precision);
break;
case ANGLE_GRAD:
if ( rat_gt( *pa, rat_200, precision) )
{
subrat(pa,rat_400, precision);
}
divrat( pa, rat_200, precision);
mulrat( pa, pi, precision);
break;
}
_sinrat( pa, precision);
}
//-----------------------------------------------------------------------------
//
// FUNCTION: cosrat, _cosrat
//
// ARGUMENTS: x PRAT representation of number to take the cosine of
//
// RETURN: cosin of x in PRAT form.
//
// EXPLANATION: This uses Taylor series
//
// n
// ___ 2j j
// \ ] X -1
// \ ---------
// / (2j)!
// /__]
// j=0
// or,
// n
// ___ 2
// \ ] -X
// \ thisterm ; where thisterm = thisterm * ---------
// / j j+1 j (2j)*(2j+1)
// /__]
// j=0
//
// thisterm = 1 ; and stop when thisterm < precision used.
// 0 n
//
//-----------------------------------------------------------------------------
void _cosrat( PRAT *px, uint32_t radix, int32_t precision)
{
CREATETAYLOR();
destroynum(pret->pp);
destroynum(pret->pq);
pret->pp=longtonum( 1L, radix);
pret->pq=longtonum( 1L, radix);
DUPRAT(thisterm,pret)
n2=longtonum(0L, radix);
xx->pp->sign *= -1;
do {
NEXTTERM(xx,INC(n2) DIVNUM(n2) INC(n2) DIVNUM(n2), precision);
} while ( !SMALL_ENOUGH_RAT( thisterm, precision) );
DESTROYTAYLOR();
// Since *px might be epsilon above 1 or below -1, due to TRIMIT we need
// this trick here.
inbetween(px, rat_one, precision);
// Since *px might be epsilon near zero we must set it to zero.
if ( rat_le(*px, rat_smallest, precision) && rat_ge(*px, rat_negsmallest, precision) )
{
DUPRAT(*px,rat_zero);
}
}
void cosrat( PRAT *px, uint32_t radix, int32_t precision)
{
scale2pi(px, radix, precision);
_cosrat(px, radix, precision);
}
void cosanglerat( _Inout_ PRAT *pa, ANGLE_TYPE angletype, uint32_t radix, int32_t precision)
{
scalerat( pa, angletype, radix, precision);
switch ( angletype )
{
case ANGLE_DEG:
if ( rat_gt( *pa, rat_180, precision) )
{
PRAT ptmp= nullptr;
DUPRAT(ptmp,rat_360);
subrat(&ptmp, *pa, precision);
destroyrat(*pa);
*pa=ptmp;
}
divrat( pa, rat_180, precision);
mulrat( pa, pi, precision);
break;
case ANGLE_GRAD:
if ( rat_gt( *pa, rat_200, precision) )
{
PRAT ptmp= nullptr;
DUPRAT(ptmp,rat_400);
subrat(&ptmp, *pa, precision);
destroyrat(*pa);
*pa=ptmp;
}
divrat( pa, rat_200, precision);
mulrat( pa, pi, precision);
break;
}
_cosrat( pa, radix, precision);
}
//-----------------------------------------------------------------------------
//
// FUNCTION: tanrat, _tanrat
//
// ARGUMENTS: x PRAT representation of number to take the tangent of
//
// RETURN: tan of x in PRAT form.
//
// EXPLANATION: This uses sinrat and cosrat
//
//-----------------------------------------------------------------------------
void _tanrat( PRAT *px, uint32_t radix, int32_t precision)
{
PRAT ptmp= nullptr;
DUPRAT(ptmp,*px);
_sinrat(px, precision);
_cosrat(&ptmp, radix, precision);
if ( zerrat( ptmp ) )
{
destroyrat(ptmp);
throw( CALC_E_DOMAIN );
}
divrat(px, ptmp, precision);
destroyrat(ptmp);
}
void tanrat( PRAT *px, uint32_t radix, int32_t precision)
{
scale2pi(px, radix, precision);
_tanrat(px, radix, precision);
}
void tananglerat( _Inout_ PRAT *pa, ANGLE_TYPE angletype, uint32_t radix, int32_t precision)
{
scalerat( pa, angletype, radix, precision);
switch ( angletype )
{
case ANGLE_DEG:
if ( rat_gt( *pa, rat_180, precision) )
{
subrat(pa, rat_180, precision);
}
divrat( pa, rat_180, precision);
mulrat( pa, pi, precision);
break;
case ANGLE_GRAD:
if ( rat_gt( *pa, rat_200, precision) )
{
subrat(pa, rat_200, precision);
}
divrat( pa, rat_200, precision);
mulrat( pa, pi, precision);
break;
}
_tanrat( pa, radix, precision);
}

View File

@@ -0,0 +1,231 @@
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.
//-----------------------------------------------------------------------------
// Package Title ratpak
// File transh.c
// Copyright (C) 1995-96 Microsoft
// Date 01-16-95
//
//
// Description
//
// Contains hyperbolic sin, cos, and tan for rationals.
//
//
//-----------------------------------------------------------------------------
#include "pch.h"
#include "ratpak.h"
bool IsValidForHypFunc(PRAT px, int32_t precision)
{
PRAT ptmp = nullptr;
bool bRet = true;
DUPRAT(ptmp,rat_min_exp);
divrat(&ptmp, rat_ten, precision);
if ( rat_lt( px, ptmp, precision) )
{
bRet = false;
}
destroyrat( ptmp );
return bRet;
}
//-----------------------------------------------------------------------------
//
// FUNCTION: sinhrat, _sinhrat
//
// ARGUMENTS: x PRAT representation of number to take the sine hyperbolic
// of
// RETURN: sinh of x in PRAT form.
//
// EXPLANATION: This uses Taylor series
//
// n
// ___ 2j+1
// \ ] X
// \ ---------
// / (2j+1)!
// /__]
// j=0
// or,
// n
// ___ 2
// \ ] X
// \ thisterm ; where thisterm = thisterm * ---------
// / j j+1 j (2j)*(2j+1)
// /__]
// j=0
//
// thisterm = X ; and stop when thisterm < precision used.
// 0 n
//
// if x is bigger than 1.0 (e^x-e^-x)/2 is used.
//
//-----------------------------------------------------------------------------
void _sinhrat( PRAT *px, int32_t precision)
{
if ( !IsValidForHypFunc(*px, precision))
{
// Don't attempt exp of anything large or small
throw( CALC_E_DOMAIN );
}
CREATETAYLOR();
DUPRAT(pret,*px);
DUPRAT(thisterm,pret);
DUPNUM(n2,num_one);
do {
NEXTTERM(xx,INC(n2) DIVNUM(n2) INC(n2) DIVNUM(n2), precision);
} while ( !SMALL_ENOUGH_RAT( thisterm, precision) );
DESTROYTAYLOR();
}
void sinhrat( PRAT *px, uint32_t radix, int32_t precision)
{
PRAT tmpx= nullptr;
if ( rat_ge( *px, rat_one, precision) )
{
DUPRAT(tmpx,*px);
exprat(px, radix, precision);
tmpx->pp->sign *= -1;
exprat(&tmpx, radix, precision);
subrat( px, tmpx, precision);
divrat( px, rat_two, precision);
destroyrat( tmpx );
}
else
{
_sinhrat( px, precision);
}
}
//-----------------------------------------------------------------------------
//
// FUNCTION: coshrat
//
// ARGUMENTS: x PRAT representation of number to take the cosine
// hyperbolic of
//
// RETURN: cosh of x in PRAT form.
//
// EXPLANATION: This uses Taylor series
//
// n
// ___ 2j
// \ ] X
// \ ---------
// / (2j)!
// /__]
// j=0
// or,
// n
// ___ 2
// \ ] X
// \ thisterm ; where thisterm = thisterm * ---------
// / j j+1 j (2j)*(2j+1)
// /__]
// j=0
//
// thisterm = 1 ; and stop when thisterm < precision used.
// 0 n
//
// if x is bigger than 1.0 (e^x+e^-x)/2 is used.
//
//-----------------------------------------------------------------------------
void _coshrat( PRAT *px, uint32_t radix, int32_t precision)
{
if ( !IsValidForHypFunc(*px, precision))
{
// Don't attempt exp of anything large or small
throw( CALC_E_DOMAIN );
}
CREATETAYLOR();
pret->pp=longtonum( 1L, radix);
pret->pq=longtonum( 1L, radix);
DUPRAT(thisterm,pret)
n2=longtonum(0L, radix);
do {
NEXTTERM(xx,INC(n2) DIVNUM(n2) INC(n2) DIVNUM(n2), precision);
} while ( !SMALL_ENOUGH_RAT( thisterm, precision) );
DESTROYTAYLOR();
}
void coshrat( PRAT *px, uint32_t radix, int32_t precision)
{
PRAT tmpx= nullptr;
(*px)->pp->sign = 1;
(*px)->pq->sign = 1;
if ( rat_ge( *px, rat_one, precision) )
{
DUPRAT(tmpx,*px);
exprat(px, radix, precision);
tmpx->pp->sign *= -1;
exprat(&tmpx, radix, precision);
addrat( px, tmpx, precision);
divrat( px, rat_two, precision);
destroyrat( tmpx );
}
else
{
_coshrat( px, radix, precision);
}
// Since *px might be epsilon below 1 due to TRIMIT
// we need this trick here.
if ( rat_lt(*px, rat_one, precision) )
{
DUPRAT(*px,rat_one);
}
}
//-----------------------------------------------------------------------------
//
// FUNCTION: tanhrat
//
// ARGUMENTS: x PRAT representation of number to take the tangent
// hyperbolic of
//
// RETURN: tanh of x in PRAT form.
//
// EXPLANATION: This uses sinhrat and coshrat
//
//-----------------------------------------------------------------------------
void tanhrat( PRAT *px, uint32_t radix, int32_t precision)
{
PRAT ptmp= nullptr;
DUPRAT(ptmp,*px);
sinhrat(px, radix, precision);
coshrat(&ptmp, radix, precision);
mulnumx(&((*px)->pp),ptmp->pq);
mulnumx(&((*px)->pq),ptmp->pp);
destroyrat(ptmp);
}