* Bugs 7562 + 15517 fixed: factorial() improved
[scilab.git] / scilab / modules / elementary_functions / tests / unit_tests / factorial.tst
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
3 // Copyright (C) 2018 - Samuel GOUGEON
4 //
5 // Copyright (C) 2012 - 2016 - Scilab Enterprises
6 //
7 // This file is hereby licensed under the terms of the GNU GPL v2.0,
8 // pursuant to article 5.3.4 of the CeCILL v.2.1.
9 // This file was originally licensed under the terms of the CeCILL v2.1,
10 // and continues to be available under such terms.
11 // For more information, see the COPYING file which you should have received
12 // along with this program.
13 // ===========================================================================
14 // <-- CLI SHELL MODE -->
15 // <-- NO CHECK REF -->
16
17 x = factorial ( 1 );
18 assert_checkequal ( x , 1 );
19 x = factorial ( [1 2 3 4] );
20 assert_checkequal ( x , [1 2 6 24] );
21 x = factorial ( [1 2 3 4]' );
22 assert_checkequal ( x , [1 2 6 24]' );
23 x = factorial ( 170 );
24 assert_checkalmostequal ( x , 7.25741561530799896739e306 , 10 * %eps );
25 // Test with a matrix
26 n = [
27 1 2 3
28 4 5 6
29 7 8 9
30 ];
31 x = factorial ( n );
32 e = [
33     1.       2.        6.
34     24.      120.      720.
35     5040.    40320.    362880.
36 ];
37 assert_checkequal ( x , e );
38
39 // Test with an hypermatrix
40 clear n expected
41 n(2,2,2) = 1;
42 n(1:8) = 1:8;
43 expected(2,2,2) = 1;
44 expected(1:8) = e'(1:8)
45 assert_checkequal(factorial(n), expected);
46
47 // -----------
48 // Accurate values from Wolfram Alpha
49 expected = [
50    1.
51    1.
52    2.
53    6.
54    24.
55    120.
56    720.
57    5040.
58    40320.
59    362880.
60    3628800.
61    39916800.
62    479001600.
63    6227020800.
64    87178291200.
65    1307674368000.
66    20922789888000.
67    355687428096000.
68    6402373705728000.
69    121645100408832000.
70    2432902008176640000.
71    51090942171709440000.
72    1124000727777607700000.
73    25852016738884978000000.
74  ];
75 assert_checkequal(factorial(0:23)', expected);
76
77 expected = [
78    2.652528598121910D+32     8.1591528324789774D+47   3.04140932017133780D+64 ..
79    1.19785716699698918D+100  9.3326215443944153D+157  5.7133839564458546D+262
80    ];
81 assert_checkalmostequal(factorial([30 40 50 70 100 150]), expected, 5*%eps);
82
83 [f, p, m] = factorial(171);
84 assert_checkequal(f, %inf);
85 assert_checkequal(p, 309);
86 assert_checkalmostequal(m, 1.2410180702176678, 1.5e-13); // 2e-14
87
88 [f, p, m] = factorial(300);
89 assert_checkequal(f, %inf);
90 assert_checkequal(p, 614);
91 assert_checkalmostequal(m, 3.060575122164406360, 2e-13); // 1e-13
92
93 [f, p, m] = factorial(200:100:900);
94 m_exp = [7.886578673647905036  3.06057512216440636  6.40345228466238953 ..
95          1.220136825991110068  1.26557231622543074  2.42204012475027218 ..
96          7.710530113353860041  6.75268022096458416
97         ];
98 p_exp = [ 374  614  868  1134  1408  1689  1976  2269  ];
99 assert_checkequal(p, p_exp);
100 assert_checkalmostequal(m, m_exp, 5e-13);
101
102 [f, p, m] = factorial(1000:1000:1e4);
103 m_exp = [
104    4.02387260077093774  ..
105    3.31627509245063324   4.149359603437854086   1.828801951514065013 ..
106    4.22857792660554352   2.683999765726739596   8.842007956963112247 ..
107    5.18418106048087694   8.099589986687190858   2.846259680917054519
108     ];
109 p_exp = [ 2567  5735  9130  12673  16325  20065  23877  27752  31681  35659 ];
110 assert_checkequal(p, p_exp);
111 assert_checkalmostequal(m, m_exp, 6e-12);  // 1.2e-11
112
113 [f, p, m] = factorial(2e4:1e4:1e5);
114 m_exp = [
115    1.81920632023034513   2.75953724621938460   2.09169242212132363  ..
116    3.34732050959714483   1.56413770802606692   1.17681241537969008 ..
117    3.09772225166224929   1.58011915409779537   2.82422940796034787
118     ];
119 p_exp = [ 77337  121287  166713  213236  260634  308759  357506  406798  456573 ];
120 assert_checkequal(p, p_exp);
121 assert_checkalmostequal(m, m_exp, 1.5e-10);  //3e-10
122
123 [f, p, m] = factorial(round(10^(6:14)));
124 m_exp = [
125    8.26393168833124006  1.202423340051590346  1.61720379492146239  ..
126    9.90462657922299373  2.325796205673083365  3.74892859910502696  ..
127    1.40366116037375609  2.403330084340115344  1.64560205598729788
128    ];
129 p_exp = [5565708  65657059  756570556  8565705522  95657055186 ..
130          1056570551815  11565705518103  125657055180974  1356570551809682
131         ];
132 //r_err = [3e-9 3e-8 3e-7 2e-6 5e-6 4e-4 5e-3 7e-2 2e-1]; // with 1st implementation
133 //r_err = [3e-9 4e-8 2e-7 2e-6 3e-5 2e-4 2e-3 5e-2 8e-2]; // with gammaln(n+1)
134 r_err  = [1e-10 6e-8 2e-8 2e-7 2e-6 7e-6 2e-4 5e-4 8e-3];
135 for i = 1:length(m_exp)
136     assert_checkequal(p(i), p_exp(i));
137     assert_checkalmostequal(m(i), m_exp(i), r_err(i));
138 end