248 lines
6.6 KiB
C++
248 lines
6.6 KiB
C++
// 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 "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;
|
|
|
|
// Set up constants and initial conditions
|
|
PRAT ratprec = i32torat(precision);
|
|
|
|
// Find the best 'A' for convergence to the required precision.
|
|
a = i32torat(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 = i32torat(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);
|
|
const auto ratRadix = i32torat(radix);
|
|
DUPRAT(tmp, ratRadix);
|
|
lograt(&tmp, precision);
|
|
subrat(&term, tmp, precision);
|
|
precision += rattoi32(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 = i32tonum(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
|
|
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(_Inout_ 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)) && (SIGN(*px) == -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);
|
|
}
|