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_binomial.h"
      16              : 
      17              : #include <stdint.h>
      18              : 
      19              : /* The code is borrowed from RosettaCode:
      20              :  * https://rosettacode.org/wiki/Evaluate_binomial_coefficients#C
      21              :  * We go to some effort to handle overflow situations.
      22              :  */
      23              : 
      24              : static uint64_t gcd_ui(uint64_t x, uint64_t y);
      25              : 
      26              : /*
      27              :  * returns binomial coefficient:
      28              :  * ( n )
      29              :  * ( k )
      30              :  */
      31    223426091 : uint64_t libgrpp_binomial(uint64_t n, uint64_t k) {
      32    223426091 :   uint64_t d, g, r = 1;
      33              : 
      34    223426091 :   if (k == 0) {
      35              :     return 1;
      36              :   }
      37    169724120 :   if (k == 1) {
      38              :     return n;
      39              :   }
      40    142204491 :   if (k >= n) {
      41     65667054 :     return (k == n);
      42              :   }
      43     76537437 :   if (k > n / 2) {
      44     37216833 :     k = n - k;
      45              :   }
      46    361244604 :   for (d = 1; d <= k; d++) {
      47    288523203 :     if (r >= UINT64_MAX / n) { /* Possible overflow */
      48      3816036 :       uint64_t nr, dr;         /* reduced numerator / denominator */
      49      3816036 :       g = gcd_ui(n, d);
      50      3816036 :       nr = n / g;
      51      3816036 :       dr = d / g;
      52      3816036 :       g = gcd_ui(r, dr);
      53      3816036 :       r = r / g;
      54      3816036 :       dr = dr / g;
      55      3816036 :       if (r >= UINT64_MAX / nr)
      56              :         return 0; /* Unavoidable overflow */
      57            0 :       r *= nr;
      58            0 :       r /= dr;
      59            0 :       n--;
      60              :     } else {
      61    284707167 :       r *= n--;
      62    284707167 :       r /= d;
      63              :     }
      64              :   }
      65              :   return r;
      66              : }
      67              : 
      68      7632072 : static uint64_t gcd_ui(uint64_t x, uint64_t y) {
      69      7632072 :   uint64_t t;
      70              : 
      71      7632072 :   if (y < x) {
      72      3816036 :     t = x;
      73      3816036 :     x = y;
      74      3816036 :     y = t;
      75              :   }
      76     19080180 :   while (y > 0) {
      77     11448108 :     t = y;
      78     11448108 :     y = x % y;
      79     11448108 :     x = t; /* y1 <- x0 % y0 ; x1 <- y0 */
      80              :   }
      81      7632072 :   return x;
      82              : }