// 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;
    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);
}