Line data    Source code

       1             : /*----------------------------------------------------------------------------*/
       2             : /*  CP2K: A general program to perform molecular dynamics simulations         */
       3             : /*  Copyright 2000-2025 CP2K developers group <https://cp2k.org>              */
       4             : /*                                                                            */
       5             : /*  SPDX-License-Identifier: MIT                                              */
       6             : /*----------------------------------------------------------------------------*/
       7             : 
       8             : /*
       9             :  *  libgrpp - a library for the evaluation of integrals over
      10             :  *            generalized relativistic pseudopotentials.
      11             :  *
      12             :  *  Copyright (C) 2021-2023 Alexander Oleynichenko
      13             :  */
      14             : 
      15             : #include "grpp_factorial.h"
      16             : 
      17             : #include <assert.h>
      18             : #include <stdint.h>
      19             : 
      20             : static uint64_t pretabulated_factorials[] = {1,
      21             :                                              1,
      22             :                                              2,
      23             :                                              6,
      24             :                                              24,
      25             :                                              120,
      26             :                                              720,
      27             :                                              5040,
      28             :                                              40320,
      29             :                                              362880,
      30             :                                              3628800,
      31             :                                              39916800,
      32             :                                              479001600,
      33             :                                              6227020800,
      34             :                                              87178291200,
      35             :                                              1307674368000,
      36             :                                              20922789888000,
      37             :                                              355687428096000,
      38             :                                              6402373705728000,
      39             :                                              121645100408832000,
      40             :                                              2432902008176640000};
      41             : 
      42   123077355 : double libgrpp_factorial(int n) {
      43   123077355 :   if (n < 0) {
      44             :     return 1;
      45   123077355 :   } else if (n <= 20) {
      46     7948411 :     return (double)pretabulated_factorials[n];
      47             :   } else {
      48   115128944 :     return n * libgrpp_factorial(n - 1);
      49             :   }
      50             : }
      51             : 
      52             : /*
      53             :  * Calculates ratio of two factorials:
      54             :  *   n!
      55             :  *  ----
      56             :  *   m!
      57             :  */
      58    16155986 : double libgrpp_factorial_ratio(int n, int m) {
      59    16155986 :   if (n == m) {
      60             :     return 1.0;
      61             :   }
      62    16155972 :   if (n < m) {
      63           0 :     return 1.0 / libgrpp_factorial_ratio(m, n);
      64             :   } else { // n > m
      65    16155972 :     double prod = 1.0;
      66   829640644 :     for (int i = m + 1; i <= n; i++) {
      67   813484672 :       prod *= i;
      68             :     }
      69             :     return prod;
      70             :   }
      71             : }
      72             : 
      73             : static uint64_t pretabulated_double_factorials[] = {
      74             :     1, // 0!!
      75             :     1,
      76             :     2,
      77             :     3,
      78             :     8,
      79             :     15, // 5!!
      80             :     48,
      81             :     105,
      82             :     384,
      83             :     945,
      84             :     3840, // 10!!
      85             :     10395,
      86             :     46080,
      87             :     135135,
      88             :     645120,
      89             :     2027025, // 15!!
      90             :     10321920,
      91             :     34459425,
      92             :     185794560,
      93             :     654729075,
      94             :     3715891200, // 20!!
      95             :     13749310575,
      96             :     81749606400,
      97             :     316234143225,
      98             :     1961990553600,
      99             :     7905853580625, // 25!!
     100             :     51011754393600,
     101             :     213458046676875,
     102             :     1428329123020800,
     103             :     6190283353629375,
     104             :     42849873690624000 // 30!!
     105             : };
     106             : 
     107   385758263 : double libgrpp_double_factorial(int n) {
     108   385758263 :   assert(n >= -1 && n <= 30);
     109   385758263 :   if (n == -1) {
     110             :     return 1;
     111   318726549 :   } else if (n <= 30) {
     112   318726549 :     return (double)pretabulated_double_factorials[n];
     113             :   } else {
     114           0 :     return n * libgrpp_double_factorial(n - 2);
     115             :   }
     116             : }