Line data    Source code

       1              : /*----------------------------------------------------------------------------*/
       2              : /*  CP2K: A general program to perform molecular dynamics simulations         */
       3              : /*  Copyright 2000-2026 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    548431275 : uint64_t libgrpp_binomial(uint64_t n, uint64_t k) {
      32    548431275 :   uint64_t d, g, r = 1;
      33              : 
      34    548431275 :   if (k == 0) {
      35              :     return 1;
      36              :   }
      37    425440680 :   if (k == 1) {
      38              :     return n;
      39              :   }
      40    361084863 :   if (k >= n) {
      41    165895600 :     return (k == n);
      42              :   }
      43    195189263 :   if (k > n / 2) {
      44     94703275 :     k = n - k;
      45              :   }
      46    922367836 :   for (d = 1; d <= k; d++) {
      47    736930665 :     if (r >= UINT64_MAX / n) { /* Possible overflow */
      48      9752092 :       uint64_t nr, dr;         /* reduced numerator / denominator */
      49      9752092 :       g = gcd_ui(n, d);
      50      9752092 :       nr = n / g;
      51      9752092 :       dr = d / g;
      52      9752092 :       g = gcd_ui(r, dr);
      53      9752092 :       r = r / g;
      54      9752092 :       dr = dr / g;
      55      9752092 :       if (r >= UINT64_MAX / nr)
      56              :         return 0; /* Unavoidable overflow */
      57            0 :       r *= nr;
      58            0 :       r /= dr;
      59            0 :       n--;
      60              :     } else {
      61    727178573 :       r *= n--;
      62    727178573 :       r /= d;
      63              :     }
      64              :   }
      65              :   return r;
      66              : }
      67              : 
      68     19504184 : static uint64_t gcd_ui(uint64_t x, uint64_t y) {
      69     19504184 :   uint64_t t;
      70              : 
      71     19504184 :   if (y < x) {
      72      9752092 :     t = x;
      73      9752092 :     x = y;
      74      9752092 :     y = t;
      75              :   }
      76     48760460 :   while (y > 0) {
      77     29256276 :     t = y;
      78     29256276 :     y = x % y;
      79     29256276 :     x = t; /* y1 <- x0 % y0 ; x1 <- y0 */
      80              :   }
      81     19504184 :   return x;
      82              : }