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   171874909 : uint64_t libgrpp_binomial(uint64_t n, uint64_t k) {
      32   171874909 :   uint64_t d, g, r = 1;
      33             : 
      34   171874909 :   if (k == 0) {
      35             :     return 1;
      36             :   }
      37   131251759 :   if (k == 1) {
      38             :     return n;
      39             :   }
      40   110357799 :   if (k >= n) {
      41    50880212 :     return (k == n);
      42             :   }
      43    59477587 :   if (k > n / 2) {
      44    28894895 :     k = n - k;
      45             :   }
      46   280864964 :   for (d = 1; d <= k; d++) {
      47   224355405 :     if (r >= UINT64_MAX / n) { /* Possible overflow */
      48     2968028 :       uint64_t nr, dr;         /* reduced numerator / denominator */
      49     2968028 :       g = gcd_ui(n, d);
      50     2968028 :       nr = n / g;
      51     2968028 :       dr = d / g;
      52     2968028 :       g = gcd_ui(r, dr);
      53     2968028 :       r = r / g;
      54     2968028 :       dr = dr / g;
      55     2968028 :       if (r >= UINT64_MAX / nr)
      56             :         return 0; /* Unavoidable overflow */
      57           0 :       r *= nr;
      58           0 :       r /= dr;
      59           0 :       n--;
      60             :     } else {
      61   221387377 :       r *= n--;
      62   221387377 :       r /= d;
      63             :     }
      64             :   }
      65             :   return r;
      66             : }
      67             : 
      68     5936056 : static uint64_t gcd_ui(uint64_t x, uint64_t y) {
      69     5936056 :   uint64_t t;
      70             : 
      71     5936056 :   if (y < x) {
      72     2968028 :     t = x;
      73     2968028 :     x = y;
      74     2968028 :     y = t;
      75             :   }
      76    14840140 :   while (y > 0) {
      77     8904084 :     t = y;
      78     8904084 :     y = x % y;
      79     8904084 :     x = t; /* y1 <- x0 % y0 ; x1 <- y0 */
      80             :   }
      81     5936056 :   return x;
      82             : }