The Wayback Machine - https://web.archive.org/web/20240112044106/https://github.com/stdlib-js/stdlib/pull/72
Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add double factorial function #72

Closed
wants to merge 6 commits into from

Conversation

chrisdamba
Copy link
Contributor

@chrisdamba chrisdamba commented Dec 17, 2016

Resolves #44 .

Checklist

Please ensure the following tasks are completed before submitting this pull request.

  • Read, understood, and followed the contributing guidelines, including the relevant style guides.
  • Read and understand the Code of Conduct.
  • Read and understood the licensing terms.
  • Searched for existing issues and pull requests before submitting this pull request.
  • Filed an issue (or an issue already existed) prior to submitting this pull request.
  • Rebased onto latest develop.
  • Submitted against develop branch.

Description

What is the purpose of this pull request?

This pull request:

  • Provides an implementation for the double factorial function

Related Issues

Does this pull request have any related issues?

This pull request:

Questions

Any questions for reviewers of this pull request?

Follow-up questions are in the responses to the review suggestions.

Other

Any other information relevant to this pull request? This may include screenshots, references, and/or implementation notes.

Encountering errors when trying to use factorial2 method in Scipy for generating the fixtures.


@stdlib-js/reviewers

Copy link
Member

@Planeshifter Planeshifter left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you very much for this pull request! I left some comments.


The double factorial should not be confused with the factorial function iterated twice, written as (n!)! and not n!!

The [double factorial][double-factorial-function] is implemented for integer values in terms of the [Factorial][@stdlib/math/base/special/factorial] function using the relation:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing [@stdlib/math/base/special/factorial]: https://github.com/stdlib-js/stdlib in the link section at the end of the file.

'use strict';

var incrspace = require( '@stdlib/math/generics/utils/incrspace' );
var factorial = require( './../lib' );
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

factorial should be replaced by factorial2 in this file.

function factorial2( x ) {
var n;
if (
(isInteger( x ) && x < 0) ||
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use isNegativeInteger module instead.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Alternatively, keep isInteger, but cache the result so you can use the boolean later at line 66.

var isInt;

isInt = isInteger( x );
if ( isInt && x < 0 ) {
   return NaN;
}
...
if ( isInt ) {
  ...
}
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good call :)

@@ -0,0 +1,38 @@
{
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The contents of the package.json files have changed, take a look at the snippet. Fields that need to be updated: author, contributors, repository, homepage, bugs, and license.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not so sure what values that need to be updated for this.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

{
  "name": "@stdlib/math/base/special/factorial2",
  "version": "0.0.0",
  "description": "Evaluate a double factorial.",
  "author": {
    "name": "The Stdlib Authors",
    "url": "https://github.com/stdlib-js/stdlib/graphs/contributors"
  },
  "contributors": [
    {
      "name": "The Stdlib Authors",
      "url": "https://github.com/stdlib-js/stdlib/graphs/contributors"
    }
  ],
  "scripts": {},
  "main": "./lib",
  "repository": {
    "type": "git",
    "url": "git://github.com/stdlib-js/stdlib.git"
  },
  "homepage": "https://github.com/stdlib-js/stdlib",
  "keywords": [
    "stdlib",
    "stdmath",
    "mathematics",
    "math",
    "special functions",
    "special",
    "function",
    "factorial",
    "double factorial",
    "fact",
    "combinatorics",
    "gamma",
    "number"
  ],
  "bugs": {
    "url": "https://github.com/stdlib-js/stdlib/issues"
  },
  "dependencies": {},
  "devDependencies": {},
  "engines": {
    "node": ">=0.10.0"
  },
  "license": "Apache-2.0"
}

var integers = require( './fixtures/julia/integers.json' );
var decimals = require( './fixtures/julia/decimals.json' );

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add empty line before // TESTS //.

} else {
delta = abs( v - expected[ i ] );
tol = EPS * Math.max( 1.0, abs( v ), abs( expected[ i ] ) );
t.ok( delta <= tol, 'within tolerance. x: ' + x[ i ] + '. Value: ' + v + '. Expected: ' + expected[ i ] + '. Tolerance: ' + tol + '.' );
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you replace the last two lines with

tol = EPS * abs( expected[ i ] );
t.ok( delta <= tol, 'within tolerance. x: '+x[i]+'. Value: '+v+'. E: '+expected[i]+'. Δ: '+delta+'. tol: '+tol );

and check whether the tests still pass?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This passes.

v = factorial2( x[ i ] );
delta = abs( v - expected[ i ] );
tol = 2.0e-14 * Math.max( 1.0, abs( v ), abs( expected[ i ] ) );
t.ok( delta <= tol, 'within tolerance. x: ' + x[ i ] + '. Value: ' + v + '. Expected: ' + expected[ i ] + '. Tolerance: ' + tol + '.' );
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you replace the last two lines with

tol = 2.0e-14 * abs( expected[ i ] );
t.ok( delta <= tol, 'within tolerance. x: '+x[i]+'. Value: '+v+'. E: '+expected[i]+'. Δ: '+delta+'. tol: '+tol );

and check whether the tests still pass?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The tolerance should be scaled via epsilon.

tol = <scalar> * EPS * abs( expected[ i ] );
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I scale the tolerance with epsilon, the tests break.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you explain the rationale behind scaling the tolerance with epsilon? Also, why can't we apply the same for the factorial tests?

@kgryte kgryte added feature Feature requests. math Specific to mathematical functionality. labels Dec 17, 2016
@kgryte
Copy link
Member

kgryte commented Dec 17, 2016

Awesome. @chrisdamba Can you also complete the PR template checklist and update the description section.


<!-- </equation> -->

The double factorial should not be confused with the factorial function iterated twice, written as (n!)! and not n!!
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that the double factorial should not...

<br>
</div>

The [Gamma][@stdlib/math/base/special/gamma] function extends the [double factorial][double-factorial-function] function for non-integer values.
Copy link
Member

@kgryte kgryte Dec 17, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

...function to non-integer values.

var DOUBLE_FACTORIALS = require( './doublefactorials.json' );
var NINF = require( '@stdlib/math/constants/float64-ninf' );
var PINF = require( '@stdlib/math/constants/float64-pinf' );

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add another newline.

var isInteger = require( '@stdlib/math/base/utils/is-integer' );
var isNegativeZero = require( '@stdlib/math/base/utils/is-negative-zero' );
var sqrt = require( '@stdlib/math/base/special/sqrt' );
var DOUBLE_FACTORIALS = require( './doublefactorials.json' );
Copy link
Member

@kgryte kgryte Dec 17, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be the last require statement. Order of requires:

  • Builtins; e.g., path, Buffer, etc.
  • External packages; e.g., debug, object-keys, etc.
  • Internal packages; e.g., @stdlib/math/base/special/sqrt.
  • Relative modules; e.g., ./doublefactorials.json.
@@ -0,0 +1,97 @@
#!/usr/bin/env julia
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might be better to use Python and SciPy for generating test fixtures, as SciPy actually has a factorial2 implementation against which we can test.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think so too but I've discovered some errors whilst trying to generate the fixtures, as in the attached screenshot.
selection_001

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will try to reproduce locally.

Copy link
Member

@kgryte kgryte Dec 18, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did some investigating. Based on the SciPy and Boost implementations, the double factorial is simply not implemented for doubles. In the code above, if you cast x to integers, then the fixtures work.

As a result, unlikely to get either of those two reference implementations to provide test fixtures.

Further, I am willing to hazard a guess that the Gamma relation is not straightforward. Boost only applies a Gamma approximation for odd integers. SciPy applies two different approximations based on whether an input integer is odd or even. While the double factorial can be extended to reals, neither of those two implementations can serve as reference implementations for that particular use case.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So in this case do we leave the fixtures as they are?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Based on my assessment, looks like the package implementation needs refactoring. Notably, unless we know of some implementation which actually does extend the double factorial function to all reals, we should stick to Boost and Scipy and only support integer inputs. For large inputs, similar to Boost, if not covered by the table, we can then apply a Gamma approximation*.

Once that is updated, the fixtures can be updated, and we'll have a better idea of where we stand.

  • NOTE: I would be curious to know how the table lookup compares to integer multiplication in terms of performance as done in SciPy. Did you ever investigate this?
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Definitely concur, needs some refactoring. I haven't actually done the investigation but I believe it's worth doing the comparison. Will do so and update you the results.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Awesome. Thanks! Really appreciate your tending to this. 👍

@kgryte
Copy link
Member

kgryte commented Dec 17, 2016

@chrisdamba Once you make the requested changes, if you could squash the commits into a single commit, that would be awesome! This would make the commit history a bit cleaner for us, so that we can identify when a unit of functionality was initially added.

Thanks again for all your work!

@codecov-io
Copy link

Current coverage is 88.2077% (diff: 96.6666%)

Merging #72 into develop will increase coverage by 0.0086%

@@            develop        #72   diff @@
==========================================
  Files          2202       2204     +2   
  Lines         29201      29231    +30   
  Methods           0          0          
  Messages          0          0          
  Branches          0          0          
==========================================
+ Hits          25755      25784    +29   
- Misses         3446       3447     +1   
  Partials          0          0          

Sunburst

Diff Coverage File Path
••••••••• 96% new ...tdlib/math/base/special/factorial2/lib/factorial2.js
•••••••••• 100% new ...es/@stdlib/math/base/special/factorial2/lib/index.js

Powered by Codecov. Last update 875afa1...a7fb88f

@@ -0,0 +1,38 @@
{
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

{
  "name": "@stdlib/math/base/special/factorial2",
  "version": "0.0.0",
  "description": "Evaluate a double factorial.",
  "author": {
    "name": "The Stdlib Authors",
    "url": "https://github.com/stdlib-js/stdlib/graphs/contributors"
  },
  "contributors": [
    {
      "name": "The Stdlib Authors",
      "url": "https://github.com/stdlib-js/stdlib/graphs/contributors"
    }
  ],
  "scripts": {},
  "main": "./lib",
  "repository": {
    "type": "git",
    "url": "git://github.com/stdlib-js/stdlib.git"
  },
  "homepage": "https://github.com/stdlib-js/stdlib",
  "keywords": [
    "stdlib",
    "stdmath",
    "mathematics",
    "math",
    "special functions",
    "special",
    "function",
    "factorial",
    "double factorial",
    "fact",
    "combinatorics",
    "gamma",
    "number"
  ],
  "bugs": {
    "url": "https://github.com/stdlib-js/stdlib/issues"
  },
  "dependencies": {},
  "devDependencies": {},
  "engines": {
    "node": ">=0.10.0"
  },
  "license": "Apache-2.0"
}

Note that the double factorial should not be confused with the factorial function iterated twice, written as (n!)! and not n!!

The [double factorial][double-factorial-function] is implemented for integer values in terms of the [Factorial][@stdlib/math/base/special/factorial] function using the relation:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Lowercase "factorial":

...in terms of the factorial function using the relation:


The [double factorial][double-factorial-function] is implemented for integer values in terms of the [Factorial][@stdlib/math/base/special/factorial] function using the relation:

<!-- <equation class="equation" label="eq:double_factorial_function_and_factorial" align="center" raw="(2n)!! = 2^n * n!" alt="Factorial function extension via the Factorial function"> -->
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Lowercase "Factorial"

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

2^n \cdot n!


The [Gamma][@stdlib/math/base/special/gamma] function extends the [double factorial][double-factorial-function] function to non-integer values.

<!-- <equation class="equation" label="eq:double_factorial_function_and_gamma" align="center" raw="(2n-1)!! = \Gamma \left( \frac{(2n+1)}{2} \right) * \frac{2^n}{\sqrt{\pi}}" alt="Factorial function extension via the Gamma function"> -->
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do not believe you need the * (asterisk) in your equation. This applies here and immediately below.


<!-- <equation class="equation" label="eq:double_factorial_function_and_factorial" align="center" raw="(2n)!! = 2^n * n!" alt="Factorial function extension via the Factorial function"> -->

<div class="equation" align="center" data-raw-text="(2n)!! = 2^n * n!" data-equation="eq:double_factorial_function_and_factorial">
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

2^n \cdot n!

* // returns NaN
*/
function factorial2( x ) {
var n;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Variables should be listed in order of length.

var isInt;
var n;
@@ -0,0 +1,97 @@
#!/usr/bin/env julia
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will try to reproduce locally.

@@ -96,7 +97,7 @@ tape( 'the function evaluates the double factorial function (positive integers <
} else {
delta = abs( v - expected[ i ] );
tol = EPS * Math.max( 1.0, abs( v ), abs( expected[ i ] ) );
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line should be

tol = EPS * abs( expected[ i ] );
@@ -126,7 +127,7 @@ tape( 'the function evaluates the factorial function (decimal values)', function
v = factorial2( x[ i ] );
delta = abs( v - expected[ i ] );
tol = 2.0e-14 * Math.max( 1.0, abs( v ), abs( expected[ i ] ) );
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line should be

tol = <num> * EPS * abs( expected[ i ] );

For example,

tol = 100.0 * EPS * abs( expected[ i ] );

The reason why we scale epsilon is that doing so provides more context as to how much results are off by and lessens the presence of so-called "magic" numbers. For instance, if we know that

actual = <number>;
expected = 1.0;
tol = 1.0 * EPS * abs( expected );
delta = abs( actual - expected );
bool = ( delta <= tol );

we know that actual value differs from the expected by only 1 bit. Similarly for other values, although the correspondence between units epsilon and number of differing bits is not 1:1.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great explanation, thank you for this!

Copy link
Member

@kgryte kgryte Dec 19, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No problem! :)

// VARIABLES //

var MAX_FACTORIAL = 170; // TODO: consider extracting as a constant
var PIO2 = 1.57079632679489661923; // pi/2
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This can be replaced by using @stdlib/math/constants/float64-half-pi.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can see there are other areas in the library which need to be updated as well to use the constant e.g in acos, asin and atan.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed. We have an internal backlog of various clean-up tasks. :( Currently picking them off as we have time. Thanks for the reminder!

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't mind picking up this clean-up task if there's an issue not assigned.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can create an issue to this end, and, if you want to submit a separate PR addressing them, that would be great!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@chrisdamba See issue #74.

@chrisdamba chrisdamba changed the title Feature/factorial2 Add double factorial function Dec 20, 2016
@rreusser
Copy link
Member

rreusser commented May 6, 2017

I'd love to help get this merged. Am I to conclude that, at least for a start, restricting the implementation to integers is the right move? Is the function then strictly a lookup table?

@Planeshifter Planeshifter reopened this May 6, 2017
@Planeshifter
Copy link
Member

Yes, restricting to integer inputs should be fine.

A table lookup and for large inputs a Gamma approximation should work, as Athan pointed out. Maybe have a look at how Boost does it.

@kgryte
Copy link
Member

kgryte commented Sep 21, 2023

Closing this PR due to inactivity.

@kgryte kgryte closed this Sep 21, 2023
@kgryte kgryte mentioned this pull request Oct 24, 2023
1 task
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
feature Feature requests. math Specific to mathematical functionality.
5 participants