2826d37056
Fixes #202 This PR fixes code style for the project files. The Problem Different files in the project use different code style. That is not consistent and leads to harder maintenance of the project. Description of the changes: Have investigated and determined the most used code style across the given codebase Have configured IDE and applied code style to all project files. Have crafted clang-formatter config. see https://clang.llvm.org/docs/ClangFormat.html https://clang.llvm.org/docs/ClangFormatStyleOptions.html Some cases were fixed manually How changes were validated: manual/ad-hoc testing, automated testing All tests pass as before because these are only code style changes. Additional Please review, and let me know if I have any mistake in the code style. In case of any mistake, I will change the configuration and re-apply it to the project.
253 lines
6.7 KiB
C++
253 lines
6.7 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;
|
|
PRAT ratprec = nullptr;
|
|
PRAT ratRadix = nullptr;
|
|
int32_t oldprec;
|
|
|
|
// Set up constants and initial conditions
|
|
oldprec = precision;
|
|
ratprec = i32torat(oldprec);
|
|
|
|
// 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);
|
|
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
|
|
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)) && (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);
|
|
}
|