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    158414399 : double libgrpp_factorial(int n) {
      43    158414399 :   if (n < 0) {
      44              :     return 1;
      45    158414399 :   } else if (n <= 20) {
      46     10391471 :     return (double)pretabulated_factorials[n];
      47              :   } else {
      48    148022928 :     return n * libgrpp_factorial(n - 1);
      49              :   }
      50              : }
      51              : 
      52              : /*
      53              :  * Calculates ratio of two factorials:
      54              :  *   n!
      55              :  *  ----
      56              :  *   m!
      57              :  */
      58     20771982 : double libgrpp_factorial_ratio(int n, int m) {
      59     20771982 :   if (n == m) {
      60              :     return 1.0;
      61              :   }
      62     20771964 :   if (n < m) {
      63            0 :     return 1.0 / libgrpp_factorial_ratio(m, n);
      64              :   } else { // n > m
      65     20771964 :     double prod = 1.0;
      66   1066680828 :     for (int i = m + 1; i <= n; i++) {
      67   1045908864 :       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    410659055 : double libgrpp_double_factorial(int n) {
     108    410659055 :   assert(n >= -1 && n <= 30);
     109    410659055 :   if (n == -1) {
     110              :     return 1;
     111    337063671 :   } else if (n <= 30) {
     112    337063671 :     return (double)pretabulated_double_factorials[n];
     113              :   } else {
     114            0 :     return n * libgrpp_double_factorial(n - 2);
     115              :   }
     116              : }